addressing comments

Author: PrajaktaPMohite

cookbooks/continental_compression_with_imposed_faults/continental_compression_with_imposed_faults.prm

@@ -44,9 +44,9 @@ subsection Formulation
   set Temperature equation = reference density profile
 end
 
-# constraining adiabatic temperature profile in compositional fields. 
+# constraining adiabatic temperature profile in compositional fields.
 # each component referring to compositional fields. First three
-# components are referring to strain and then chemical composition fields. 
+# components are referring to strain and then chemical composition fields.
 # see subsection Compositional fields
 subsection Adiabatic conditions model
   subsection Compute profile
@@ -63,20 +63,30 @@ subsection Adiabatic conditions model
   end
 end
 
+# Specify the reference depth-dependent compositional profile used when
+# computing the adiabatic conditions. A function is required for each compositional
+# field, and the order of the functions should correspond to the same order of in
+# which the fields were listed under "set Names of fields".  As such, the first
+# three fields correspond to various types of "strain", and should be set to
+# 0, while the remaining fields correspond to fields representing chemical
+# compositions. Note that in the function expression "x" refers to depth.
+
 subsection Discretization
   set Composition polynomial degree           = 2
   set Stokes velocity polynomial degree       = 2
   set Temperature polynomial degree           = 2
-  # The gradient at quadrature points will not be a precise approximation of the true gradient 
-  # since it excludes the contribution of jumps in the compositional field between cells.
-  # Even though the finite element solution varies from cell to cell, the gradient of the 
-  # solution at quadrature locations inside each cell will always be precisely zero.
   set Use discontinuous composition discretization = true
+  # Discontinuous composition discretization is required when using viscoelastic or
+  # other history-dependent rheologies in ASPECT. This ensures that compositional fields
+  # remain cell-local and consistent with the assumptions of the constitutive law.
+  # While it may increase the number of degrees of freedom, it is essential for maintaining
+  # the accuracy and stability of stress evolution in models involving localized rheological heterogeneity.
 end
 
 # Model geometry (400x400 km, 20 km spacing)
 subsection Geometry model
   set Model name = box
+
   subsection Box
     set X repetitions = 20
     set Y repetitions = 20
@@ -95,6 +105,7 @@ subsection Mesh refinement
   set Initial global refinement          = 3
   set Time steps between mesh refinement = 0
   set Strategy = minimum refinement function
+
   subsection Minimum refinement function
     set Coordinate system   = cartesian
     set Variable names      = x,y
@@ -109,9 +120,11 @@ end
 subsection Mesh deformation
   set Mesh deformation boundary indicators        = top: free surface, top: diffusion
   set Additional tangential mesh velocity boundary indicators = left, right
+
   subsection Free surface
     set Surface velocity projection = normal
   end
+
   subsection Diffusion
     # Diffusivity term. Increasing this value will result
     # in a smoother free surface and lower topography
@@ -121,25 +134,25 @@ subsection Mesh deformation
 end
 
 # Velocity on boundaries characterized by functions
-# A horizontal velocity (x-direction) of 0.1 cm/year is applied to the left and right walls 
-# resulting in inflow, which is balanced by the 0.2 cm/yr vertical (y-direction) outflow on the bottom boundary.
+# A horizontal velocity (x-direction) of 0.1 cm/year (~0.0001 m/s)  is applied to the left and right walls
+# resulting in inflow, which is balanced by the 0.2 cm/yr  (~0.0002 m/s) vertical (y-direction) outflow on the bottom boundary.
 # Velocity components parallel to the base (x-velocity) and side walls (y-velocity)
 # are unconstrained (i.e. 'free').
+
 subsection Boundary velocity model
-  set Prescribed velocity boundary indicators = left x: function, right x:function, bottom y:function
+  set Prescribed velocity boundary indicators = left x:function, right x:function, bottom y:function
 
   subsection Function
     set Variable names      = x,y
-    set Function constants  = v=-0.0001, w=400.e3, d=400.e3    # 1mmyr
-    # negative sign indicates compression   (v=-0.0025)
+    set Function constants  = v=0.0001, w=400.e3, d=400.e3    # 0.1 cm/yr
     set Function expression = if (x < w/2 , -v, v) ; v*2*d/w
   end
 end
 
 # Number and name of compositional fields, which are tracked with particles
 # The seven compositional fields represent:
 # 1. Plastic strain that accumulated over time, including the initial plastic strain values
-# 2. Plastic strain that accumualates over time, with the initial plastic strain removed
+# 2. Plastic strain that accumualated over time, with the initial plastic strain removed
 # 3. Viscous strain that accumulated over time, including the initial viscous strain values
 # 4. Upper crust
 # 5. Lower crust
@@ -174,7 +187,7 @@ end
 
 # Initial values of different compositional fields
 # The upper crust (20 km thick), lower crust (20 km thick)
-# and mantle lithosphere(40 km thick) and asthenosphere (320 km thick )are continuous 
+# and mantle lithosphere (40 km thick) and asthenosphere (320 km thick) are continuous
 # horizontal layers of constant thickness. The non initial plastic strain is
 # set to 0, while the initial plastic and viscous strain is
 # randomized between 0.5 and 1.5.
@@ -203,16 +216,23 @@ subsection Boundary composition model
 end
 
 # Temperature boundary conditions
-# Top and bottom (fixed) temperatures are consistent with the initial temperature field
-# Note that while temperatures are specified for the model sides, these values are
-# not used as the sides are not specified "Fixed temperature boundaries".  Rather,
-# these boundaries are insulating (zero net heat flux).
-# Initial temperature model can be computed using plugin from continental_extension cookbook
+# Top and bottom (fixed) temperatures are consistent with the initial temperature field.
+# Side boundaries are set to insulating, i.e. a net zero heat flux.
+# These boundaries are insulating (i.e., zero net heat flux), so the temperature
+# field evolves due to internal heat production and conduction alone.
+# The initial geotherm defined here is consistent with these constraints and can
+# also be generated using the 'continental geotherm' plugin in ASPECT.
+# For an example of how to implement this plugin, refer to:
+# tests/continental_extension.prm
+
 subsection Boundary temperature model
   set Fixed temperature boundary indicators = bottom, top
   set List of model names = box
+
   subsection Box
     set Bottom temperature = 1669
+    # (calculated at y = 0 using the asthenosphere expression)
+    # Continuity in temperature and heat flux at each layer boundary
     set Top temperature    =  273
   end
 end
@@ -236,11 +256,13 @@ end
 #   Surface Heat Flow    - lower crust (qs2) = 0.06125 W/m^2;
 #                        - mantle      (qs3) = 0.04125 W/m^2;
 #       - asthenosphere (qs4)= 0.03625 W/m^2;
+
 # Note: The continental geotherm initial temperature model
 #       plugin can be used to compute an identical geotherm
 #       for the lithosphere. An example of how to use this
 #       plugin is illustrated in the test for this cookbook
 #       (tests/continental_extension.prm).
+
 subsection Initial temperature model
   set Model name = function
   subsection Function
@@ -323,6 +345,10 @@ subsection Material model
     # so no need to specify diffusion creep parameters below, which are
     # only used if "diffusion" or "composite" option is selected.
     set Viscous flow law = dislocation
+    # The model currently uses a dislocation creep flow law, so the diffusion creep parameters
+    # listed below are NOT active. They are included here for reference and future use —
+    # if the viscous flow law is switched to "composite", the model will incorporate both
+    # dislocation and diffusion creep mechanisms.
 
     # Dislocation creep parameters for
     # 1. Background material/mantle (dry olivine)
@@ -422,8 +448,8 @@ subsection Particles
   set Load balancing strategy     = remove and add particles
   set List of particle properties = initial composition, viscoplastic strain invariants, pT path, position
   set Interpolation scheme        = bilinear least squares
-  set Update ghost particles      = true
   set Particle generator name     = reference cell
+
   subsection Generator
     subsection Reference cell
       set Number of particles per cell per direction = 7
@@ -435,17 +461,21 @@ end
 # Post processing
 subsection Postprocess
   set List of postprocessors = basic statistics, composition statistics, heat flux densities, heat flux statistics, mass flux statistics, particles, pressure statistics, temperature statistics, topography, velocity statistics, visualization
+
   subsection Visualization
     set List of output variables = heat flux map, material properties, named additional outputs, strain rate
+
     subsection Material properties
       set List of material properties = density, viscosity
     end
+
     set Output format = vtu
     set Time between graphical output = 100.e3
     set Interpolate output            = true
     set Write higher order output     = true
     set Output format                 = vtu
   end
+
   subsection Particles
     set Time between data output    = 100.e3
     set Data output format          = vtu

cookbooks/continental_compression_with_imposed_faults/doc/continental_compression_with_imposed_faults.md

@@ -5,21 +5,21 @@
 
 Recent numerical modeling investigations have demonstrated the influence of extensional phases of deformation on the subsequent evolution of continental collision zones (eg., {cite}`Vasey:etal:2024`,{cite}`Zwaan:etal:2025`). Motivated by these studies and a wide range of observations that highlight fault reactivation processes during the Wilson cycle, this cookbook implements a 2D visco-plastic model of continental compression that incorporates pre-existing fault zones into the initial conditions.
 
-The cookbook builds directly on components of the {ref}`sec:cookbooks:crustal-deformation` cookbook, and the primary goal here is to:
+The cookbook builds directly on components of the {ref}`sec:cookbooks:continental-extension` cookbook, and the primary goal here is to:
 1. Highlight techniques for imposing fault zones in the initial conditions of a lithospheric deformation simulation and
 2. Demonstrate their effect on the evolution of deformation patterns.
 
 # Model Design
-The model domain spans 400 km x 400 km and uses adaptive refinement to resolve deformation patterns in the regions where faults are imposed at the onset of deformation. The initial thermal structure follows a conductive, continental-style geothermal through the lithosphere and an initial adiabatic profile in the asthenosphere. The governing equations follow the extended Boussinesq approximation, which includes both adiabatic and shear heating.
+The model domain spans 400 km x 400 km and uses adaptive refinement to resolve deformation patterns in the regions where faults are imposed at the onset of deformation. The initial thermal structure follows a conductive, continental-style geotherm through the lithosphere and an initial adiabatic profile in the asthenosphere. The governing equations follow the extended Boussinesq approximation, which includes both adiabatic and shear heating.
 
-Deformation is driven by horizontal velocity applied at the model sides (1 mm/yr), which are balanced by outflow at the model base (2 mm/yr). A free surface allows topography to develop through time, which is diffused at each time step to approximate landscape evolution processes and stabilize both linear and nonlinear solver behavior.
+Deformation is driven by horizontal velocity applied at the model sides (0.1 cm/yr), which are balanced by outflow at the model base (0.2 cm/yr). A free surface allows topography to develop through time, which is diffused at each time step to approximate landscape evolution processes and stabilize both linear and nonlinear solver behavior.
 
 The initial lithological structure includes distinct layers representing the upper crust (20 km), lower crust (20 km), mantle lithosphere (40 km), and asthenosphere (320 km). Each layer is represented by a nonlinear viscoplastic rheology combining dislocation creep and pressure-dependent plasticity. Respectively, the pre-yield viscosity and brittle material parameters (cohesion, friction) weaken by a factor of 10 and 4 as a function of accumulated viscous and brittle strain over strain intervals of 0.5-1.5.
 
 ```{literalinclude} strain_weakening_mechanism.part.prm
 ```
 
-Faults are correspondingly integrated into the model initial conditions as zones of initial strain defined in the Geodynamic World Builder. Four faults extending to the base of the lithosphere with constant dip angles are included, two of which dip at 60 degrees toward each other in a style representing pre-existing normal faults. A third vertical fault (i.e., approximating a strike-slip fault) is centrally located between the two dipping faults, while the fourth fault is located to the right of two normal faults and dips at a higher angle. The faults maintain a constant width of 5 km and extend to a depth of 100 km. In addition to defined fault locations with constant strain values, randomized zones of plastic (0-40 km depth) and viscous strain (0-100 km depth) are imposed across a 300 km wide zone in the upper to approximate pervasive off-fault damage observed in many regions that have undergone significant tectonic deformation. The faults are setting in the middle of the randomized strain zones, fault1 in wb file is located at x = 250km. The plastic and viscous strain compositions value within the faults is the sum of the value defined in the world builder and the randomized values defined in these zones. 
+Faults are integrated into the model initial conditions as zones of initial strain defined in the Geodynamic World Builder. Four faults extending to the base of the lithosphere with constant dip angles are included, two of which dip at 60 degrees toward each other in a style representing pre-existing normal faults. A third vertical fault (i.e., approximating a strike-slip fault) is centrally located between the two dipping faults, while the fourth fault is located to the right of two normal faults and dips at a higher angle. The faults maintain a constant width of 5 km and extend to a depth of 100 km. In addition to defined fault locations with constant strain values, randomized zones of plastic (0-40 km depth) and viscous strain (0-100 km depth) are imposed across a 300 km wide zone in the upper to approximate pervasive off-fault damage observed in many regions that have undergone significant tectonic deformation. The faults are prescribed in the middle of the randomized strain zones, fault1 in wb file is located at x = 250km. The plastic and viscous strain composition values within the faults is the sum of the value defined in the world builder and the randomized values defined in these zones.
 
 Following Howard et al. (2003), the configuration of these faults is motivated by the inferred tectonic history of the Dzereg basin in the Mongolia Altai, which has been undergoing relatively slow compression since the onset of the India-Asia collision following a period of extensional deformation.
 
@@ -49,8 +49,8 @@ The position, dip angle and direction, thickness (5 km), and composition (plasti
  Strain rate (left) and density (right) distributions with temperature contours after 35 Myr of deformation.
 ```
 
-Deformation preferentially localizes along the two pre-existing shear zones forming a grabben-like geometry ({numref}`fig:strain_rate_and_density_0_myr`), which reflects their optimal orientation and position within the model domain to accommodate the imposed convergence. After 35 Myr of convergence ({numref}`fig:strain_rate_and_density_35_myr`) deformation remains strongly localized along these faults,  with crustal shortening and thickening occurring between them forming an approximately symmetric orogenic wedge. Some degree of asymmetry develops on the right side of the wedge, where the fault farthest to the right remains active in the upper crust and connects to the fault bounding the wedge near the brittle-ductile transition zone.
+Deformation preferentially localizes along the two pre-existing shear zones forming a graben-like geometry ({numref}`fig:strain_rate_and_density_0_myr`), which reflects their optimal orientation and position within the model domain to accommodate the imposed convergence. After 35 Myr of convergence ({numref}`fig:strain_rate_and_density_35_myr`) deformation remains strongly localized along these faults, with crustal shortening and thickening occurring between them forming an approximately symmetric orogenic wedge. Some degree of asymmetry develops on the right side of the wedge, where the fault farthest to the right remains active in the upper crust and connects to the fault bounding the wedge near the brittle-ductile transition zone.
 
 Given the nonlinearity of the rheology and governing equations, minor variations in fault strength, geometry, lithospheric structure, and boundary velocities may lead to significant variations in the spatiotemporal evolution of deformation. This cookbook provides a flexible framework for exploring the effects of these parameters, and application to hypothesis-driven questions such as fault reactivation, inversion of rift basins, or the partitioning of strain in complex orogens.
 
-Similarly, changing the numerical resolution is likely to also affect the results, as the brittle shear band width by introduced by the plastic damper (1e21 Pa s) is still likely not fully resolved at the maximum resolution of 2.5 km. Furthermore, varying degrees of minor variations in the model results will likely occur if a stricter nonlinear solver tolerance is selected. This cookbook provides a flexible framework for exploring the effects of these parameters, and application to hypothesis-driven questions such as fault reactivation, inversion of rift basins, or the partitioning of strain in complex orogens.
+Similarly, changing the numerical resolution is likely to also affect the results, as the brittle shear band width introduced by the plastic damper (1e21 Pa s) is still likely not fully resolved at the maximum resolution of 2.5 km. Furthermore, minor variations in the model results will likely occur if a stricter nonlinear solver tolerance is selected. This cookbook provides a flexible framework for exploring the effects of these parameters, and application to questions such as fault reactivation, inversion of rift basins, or the partitioning of strain in complex orogens.