diff --git a/cookbooks/continental_compression_with_imposed_faults/continental_compression_with_imposed_faults.prm b/cookbooks/continental_compression_with_imposed_faults/continental_compression_with_imposed_faults.prm new file mode 100644 index 00000000000..64c9beb289b --- /dev/null +++ b/cookbooks/continental_compression_with_imposed_faults/continental_compression_with_imposed_faults.prm @@ -0,0 +1,488 @@ +# Continental Compression with Imposed Faults Cookbook +# This cookbook builds on the continental extension cookbook +# and a recent paper investigating rift inversion (Vasey et al., 2024, +# https://doi.org/10.1130/G51489.1) to demonstrate how to design +# models of continental compression that extend into the asthenosphere +# and include faults as zones of weakness in the initial conditions. +# The setup is motivated by observed deformation patterns in the Mongolian +# Altai, which has experienced low magnitudes of compressional deformation +# since the onset of the India-Asia collision and contains evidence for +# reactivation of normal faults formed during prior phases of deformation. +# (Howard et al., 2003, https://doi.org/10.1046/j.1365-2117.2003.00198.x). +# The faults locations and properties are defined using the Geodynamic +# World Builder. + + +#### Global parameters +set World builder file = continental_compression_with_imposed_faults.wb +set Dimension = 2 +set Start time = 0 +set End time = 3.5e7 +set Use years in output instead of seconds = true +set Nonlinear solver scheme = iterated Advection and Stokes +set Nonlinear solver tolerance = 1e-6 +set Max nonlinear iterations = 200 +set CFL number = 0.25 +set Maximum time step = 20e3 +set Output directory = output-continental_compression_with_imposed_faults +set Pressure normalization = no + +subsection Solver parameters + subsection Stokes solver parameters + set Stokes solver type = block AMG + set Number of cheap Stokes solver steps = 4000 + set Linear solver tolerance = 1e-8 + set GMRES solver restart length = 100 + set Use full A block as preconditioner = true + end +end + +# Governing equations +subsection Formulation + set Formulation = custom + set Mass conservation = incompressible + set Temperature equation = reference density profile +end + +# 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. +# see subsection Compositional fields +subsection Adiabatic conditions model + subsection Compute profile + set Composition reference profile = function + set Function expression = \ + 0; \ + 0; \ + 0; \ + if( x<=20.e3, 1, 0); \ + if( x>20.e3 && x<=40.e3, 1, 0); \ + if( x>40.e3 && x<=80.e3, 1, 0); \ + if( x>80.e3, 1, 0); \ + + 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 + 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 + set X extent = 400e3 + set Y extent = 400e3 + end +end + +# Globally refine the mesh to 2.5 km spacing, and then +# adaptively refine the mesh to 1.25 km spacing above y=50 km +# and between x=40 and x=160 km. These values ensure areas +# undergoing brittle deformation are in the high-resolution +# region. +subsection Mesh refinement + set Initial adaptive refinement = 1 + 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 + set Function expression = if ( y>=300e3 && x>=40.e3 && x<=280.e3, 4, 3) + end +end + +# Advecting the free surface using a normal, rather than vertical, +# projection. To reduce mesh instabilities and associated solver +# issues when deformation becomes large, diffusion is applied to +# the free surface at each time step. +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 + # amplitudes. + set Hillslope transport coefficient = 1.e-8 + end +end + +# Velocity on boundaries characterized by functions +# 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 + + subsection Function + set Variable names = x,y + 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 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 +# 6. Mantle lithosphere +# 7. Asthenosphere + +subsection Compositional fields + set Number of fields = 7 + set Names of fields = plastic_strain, \ + noninitial_plastic_strain, \ + viscous_strain, \ + crust_upper, \ + crust_lower, \ + mantle_lithosphere, \ + asthenosphere + set Types of fields = strain, \ + strain, \ + strain, \ + chemical composition, \ + chemical composition, \ + chemical composition, \ + chemical composition + set Compositional field methods = particles + set Mapped particle properties = plastic_strain: plastic_strain, \ + noninitial_plastic_strain: noninitial_plastic_strain, \ + viscous_strain: viscous_strain, \ + crust_upper: initial crust_upper, \ + crust_lower: initial crust_lower, \ + mantle_lithosphere: initial mantle_lithosphere, \ + asthenosphere: initial asthenosphere +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 +# 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. +subsection Initial composition model + set List of model names = function, world builder + set List of model operators = add, add + + subsection Function + set Variable names = x,y + set Function expression = if(x>40.e3 && x<360.e3 && y>360.e3, 0.5 + rand_seed(1), 0); \ + 0; \ + if(x>40.e3 && x<360.e3 && y>300.e3, 0.5 + rand_seed(1), 0); \ + if(y>=380.e3, 1, 0); \ + if(y<380.e3 && y>=360.e3, 1, 0); \ + if(y<360.e3 && y>=320.e3, 1, 0); \ + if(y<320.e3 && y>=0.e3, 1, 0); + + end +end + +# Composition: fixed on the bottom (outflow boundary) and side boundary, and free on top +subsection Boundary composition model + set Fixed composition boundary indicators = bottom, left, right + set List of model names = initial composition + set Allow fixed composition on outflow boundaries = true +end + +# Temperature boundary conditions +# 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 + +# Initial temperature field +# Typical continental geotherm based on equations 4-6 from: +# D.S. Chapman (1986), "Thermal gradients in the continental crust", +# Geological Society of London Special Publications, v.24, p.63-70. +# The initial constraints are: +# Surface Temperature - upper crust (ts1) = 273 K +# Surface Heat Flow - upper crust (qs1) = 0.06125 mW/m^2 +# Heat Production - upper crust (A1) = 1.00e-6 W/m^3; +# Heat Production - lower crust (A2) = 0.25e-6 W/m^3; +# Heat Production - mantle (A3) = 0.00e-6 W/m^3; +# Heat Production - asthenosphere (A4)= 0.00e-6 W/m^3; +# Thermal Conductivity - all layers = 2.5 (W/(m K)); +# To satisfy these constraints, the following values are required: +# Surface Temperature - lower crust (ts2) = 683 K +# - mantle (ts3) = 993 K +# - asthenosphere (ts4)= 1573 K +# 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 + set Variable names = x,y + set Function constants = h=400e3, ts1=273, ts2=683, ts3=993, ts4=1573,\ + A1=1.e-6, A2=0.25e-6, A3=0.0, A4=0.0,\ + k1=2.5, k2=2.5, k3=2.5, k4=2.5, \ + qs1=0.06125, qs2=0.04125, qs3=0.03625, qs4=0.03625 + set Function expression = if( (h-y)<=20.e3, \ + ts1 + (qs1/k1)*(h-y) - (A1*(h-y)*(h-y))/(2.0*k1), \ + if( (h-y)>20.e3 && (h-y)<=40.e3, \ + ts2 + (qs2/k2)*(h-y-20.e3) - (A2*(h-y-20.e3)*(h-y-20.e3))/(2.0*k2), \ + if( (h-y)>40.e3 && (h-y)<=80.e3, \ + ts3 + (qs3/k3)*(h-y-40.e3) - (A3*(h-y-40.e3)*(h-y-40.e3))/(2.0*k3), \ + ts4 + ((h-y-80.e3)/1.e3)*0.3 ) ) ); + end +end + +# Constant internal heat production values (W/m^3) for background material +# and compositional fields. +subsection Heating model + set List of model names = compositional heating, adiabatic heating, shear heating + subsection Compositional heating + set Use compositional field for heat production averaging = 1, 0, 0, 0, 1, 1, 1,1 + set Compositional heating values = 0., 0., 0., 0., 1.0e-6, 0.25e-6, 0.0, 0.0 + end + subsection Adiabatic heating + set Use simplified adiabatic heating = true + end +end + +# Material model +# Rheology: Non-linear viscous flow and Drucker Prager Plasticity +# Values for most rheological parameters are specified for a background material and +# each compositional field. Values for viscous deformation are based on dislocation +# creep flow-laws, with distinct values for the upper crust (wet quartzite), lower +# crust (wet anorthite) and mantle (dry olivine). +subsection Material model + set Model name = visco plastic + + set Material averaging = harmonic average only viscosity + + subsection Visco Plastic + + # Reference temperature and viscosity + set Reference temperature = 273 + + # The minimum strain-rate helps limit large viscosities values that arise + # as the strain-rate approaches zero. + # The reference strain-rate is used on the first non-linear iteration + # of the first time step when the velocity has not been determined yet. + set Minimum strain rate = 1.e-20 + set Reference strain rate = 1.e-16 + + # Limit the viscosity with minimum and maximum values + set Minimum viscosity = 1e18 + set Maximum viscosity = 1e26 + + # Thermal diffusivity is adjusted to match thermal conductivities + # assumed in assigning the initial geotherm + set Define thermal conductivities = true + set Thermal conductivities = 2.5 + set Heat capacities = 750. + + # The densities below are set to achieve the following densities at the reference temperature (273 K): + # upper_crust - 2700, lower_crust - 2900, background/mantle_lithosphere/asthenosphere - 3300. + set Densities = background: 3216.374269005848, \ + crust_upper: 2729.044834307992, \ + crust_lower: 2826.510721247563, \ + mantle_lithosphere: 3216.374269005848, \ + asthenosphere: 3216.374269005848 + + set Thermal expansivities = 2e-5 + + # Harmonic viscosity averaging + set Viscosity averaging scheme = harmonic + + # Choose to have the viscosity (pre-yield) follow a dislocation + # diffusion or composite flow law. Here, dislocation is selected + # 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) + # Hirth & Kohlstedt (2004), Geophys. Monogr. Am. Geophys. Soc., v.138, p.83-105. + # "Rheology of the upper mantle and the mantle wedge:a view from the experimentalists" + # 2. Upper crust (wet quartzite) + # Gleason and Tullis (1995) + # 3. Lower crust and weak seed (wet anorthite) + # Rybacki et al. (2006), J. Geophys. Res., v.111(B3). + # "Influence of water fugacity and activation volume on the flow properties of fine-grained + # anorthite aggregates" + # Note that the viscous pre-factors below are scaled to plane strain from unixial strain experiments. + # For ease of identification, fields tracking strain are assigned prefactors of 1e-50 + set Prefactors for dislocation creep = background: 7.37e-15, \ + crust_upper: 1.37e-26, \ + crust_lower: 5.71e-23, \ + mantle_lithosphere: 7.37e-15, \ + asthenosphere: 7.37e-15 + + set Stress exponents for dislocation creep = background: 3.5, \ + crust_upper: 4.0, \ + crust_lower: 3.0, \ + mantle_lithosphere: 3.5, \ + asthenosphere: 3.5 + + set Activation energies for dislocation creep = background: 530.e3, \ + crust_upper: 223.e3, \ + crust_lower: 345.e3, \ + mantle_lithosphere: 530.e3, \ + asthenosphere: 530.e3 + + set Activation volumes for dislocation creep = background: 18.e-6, \ + crust_upper: 0, \ + crust_lower: 0, \ + mantle_lithosphere: 18.e-6, \ + asthenosphere: 18.e-6 + + set Prefactors for diffusion creep = background: 2.37e-15, \ + crust_upper: 1.00e-50, \ + crust_lower: 1.00e-50, \ + mantle_lithosphere: 1.00e-50, \ + asthenosphere: 1.00e-50 + + set Grain size exponents for diffusion creep = background: 3.0, \ + crust_upper: 0., \ + crust_lower: 0., \ + mantle_lithosphere: 0., \ + asthenosphere: 0. + + set Activation energies for diffusion creep = background: 375.e3, \ + crust_upper: 0., \ + crust_lower: 0., \ + mantle_lithosphere: 0., \ + asthenosphere: 0. + + set Activation volumes for diffusion creep = background: 2.e-6, \ + crust_upper: 0, \ + crust_lower: 0, \ + mantle_lithosphere: 0, \ + asthenosphere: 0 + + # Plasticity parameters + set Angles of internal friction = 30 + set Cohesions = 20.e6 + + # The parameters below weaken the friction and cohesion by a + # a factor of 4 between plastic strain values of 0.5 and 1.5, + # and the pre-yield viscosity by a factor of 10 over the same + # strain interval. + set Strain weakening mechanism = plastic weakening with plastic strain and viscous weakening with viscous strain + set Start plasticity strain weakening intervals = 0.5 + set End plasticity strain weakening intervals = 1.5 + set Cohesion strain weakening factors = 0.25 + set Friction strain weakening factors = 0.25 + set Start prefactor strain weakening intervals = 0.5 + set End prefactor strain weakening intervals = 1.5 + set Prefactor strain weakening factors = 0.1 + + set Use plastic damper = true + set Plastic damper viscosity = 1e21 + + end +end + +# Gravity model +subsection Gravity model + set Model name = vertical + + subsection Vertical + set Magnitude = 9.81 + end +end + +subsection Particles + set Minimum particles per cell = 25 + set Maximum particles per cell = 100 + 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 Particle generator name = reference cell + + subsection Generator + subsection Reference cell + set Number of particles per cell per direction = 7 + end + end +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 + end +end + +# Checkpointing +subsection Checkpointing + set Steps between checkpoint = 10 +end diff --git a/cookbooks/continental_compression_with_imposed_faults/continental_compression_with_imposed_faults.wb b/cookbooks/continental_compression_with_imposed_faults/continental_compression_with_imposed_faults.wb new file mode 100755 index 00000000000..9b990ed5fe3 --- /dev/null +++ b/cookbooks/continental_compression_with_imposed_faults/continental_compression_with_imposed_faults.wb @@ -0,0 +1,66 @@ +{ + "version":"1.0", + "cross section":[[0,4e5],[4.2e5,4e5]], + "features": + [ + { + "model":"fault","name":"fault1","dip point":[0e5,2e5], + "min depth":0, "max depth":1e5, + "coordinates":[[2.5e5,8e5],[2.5e5,0]], + "segments": + [ + {"length":2e5,"thickness":[5000],"angle":[80]} + ], + "temperature models":[{"model":"uniform","temperature":273}], + "composition models": + [ + {"model":"smooth", "compositions":[0,2], "operation":"add", "side distance fault center":2500, "center fractions":[1.5,1.5],"side fractions":[0.5,0.5]} + ] + + }, { + "model":"fault","name":"fault2","dip point":[0e5,2e5], + "min depth":0, "max depth":1e5, + "coordinates":[[2.25e5,8e5],[2.25e5,0]], + "segments": + [ + {"length":2e5,"thickness":[5000],"angle":[70]} + ], + "temperature models":[{"model":"uniform","temperature":273}], + "composition models": + [ + {"model":"smooth", "compositions":[0,2], "operation":"add", "side distance fault center":2500, "center fractions":[1.5,1.5],"side fractions":[0.5,0.5]} + ] + + }, { + "model":"fault","name":"fault3","dip point":[0e5,2e5], + "min depth":0, "max depth":1e5, + "coordinates":[[1.4e5,8e5],[1.4e5,0]], + "segments": + [ + {"length":2e5,"thickness":[5000],"angle":[110]} + ], + "temperature models":[{"model":"uniform","temperature":273}], + "composition models": + [ + {"model":"smooth", "compositions":[0,2], "operation":"add", "side distance fault center":2500, "center fractions":[1.5,1.5],"side fractions":[0.5,0.5]} + ] + + }, { + "model":"fault","name":"fault4","dip point":[0e5,2e5], + "min depth":0, "max depth":1e5, + "coordinates":[[1.85e5,8e5],[1.85e5,0]], + "segments": + [ + {"length":2e5,"thickness":[5000],"angle":[90]} + ], + "temperature models":[{"model":"uniform","temperature":273}], + "composition models": + [ + {"model":"smooth", "compositions":[0,2], "operation":"add", "side distance fault center":2500, "center fractions":[1.5,1.5],"side fractions":[0.5,0.5]} + ] + + } + + ] +} + diff --git a/cookbooks/continental_compression_with_imposed_faults/doc/continental_compression_with_imposed_faults.md b/cookbooks/continental_compression_with_imposed_faults/doc/continental_compression_with_imposed_faults.md new file mode 100644 index 00000000000..67f633ed1b0 --- /dev/null +++ b/cookbooks/continental_compression_with_imposed_faults/doc/continental_compression_with_imposed_faults.md @@ -0,0 +1,56 @@ +(sec:cookbooks:Continental-Compression-with-imposed-faults)= +# Continental Compression with Pre-Existing Fault inheritance + +*This section was contributed by Prajakta Mohite and John Naliboff.* + +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: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 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 (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 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. + +```{figure-md} fig:initial_plastic_and_viscous_strain_and_density + + + Initial plastic strain (top left), viscous strain (top right) and density (bottom) highlighting the location of defined fault zones and randomized strain across a broader region in the upper 100 km. + The density plot contains temperature contours at intervals of 200 K, beginning at 373 K and ending at 1573 K (LAB temperature). +``` + +The position, dip angle and direction, thickness (5 km), and composition (plastic and viscous strain) of each fault are defined with the [Geodynamic World Builder fault feature](https://geodynamicworldbuilder.github.io), which is illustrated below for a single fault. + +```{literalinclude} single_fault_imposed.part.wb +``` + +# Model Evolution and Potential Expansion. + +```{figure-md} fig:strain_rate_and_density_0_myr + + + Strain rate (left) and density distribution (right) at 0 Myr in the upper 100 km from x = 100 to 300 km. Strain localizes in the upper crust along pre-defined fault zones, while density increases with depth, reflecting the compositional and thermal stratification of the lithosphere. +``` + +```{figure-md} fig:strain_rate_and_density_35_myr + + + 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 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 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. diff --git a/cookbooks/continental_compression_with_imposed_faults/doc/initial_plastic_and_viscous_strain_and_density.svg b/cookbooks/continental_compression_with_imposed_faults/doc/initial_plastic_and_viscous_strain_and_density.svg new file mode 100644 index 00000000000..b468189a3b6 --- /dev/null +++ b/cookbooks/continental_compression_with_imposed_faults/doc/initial_plastic_and_viscous_strain_and_density.svg @@ -0,0 +1 @@ + \ No newline at end of file diff --git a/cookbooks/continental_compression_with_imposed_faults/doc/single_fault_imposed.part.wb b/cookbooks/continental_compression_with_imposed_faults/doc/single_fault_imposed.part.wb new file mode 100644 index 00000000000..8a6de9e7779 --- /dev/null +++ b/cookbooks/continental_compression_with_imposed_faults/doc/single_fault_imposed.part.wb @@ -0,0 +1,23 @@ +{ + "version":"1.0", + "cross section":[[0,4e5],[4.2e5,4e5]], + "features": + [ + { + "model":"fault","name":"fault1","dip point":[0e5,2e5], + "min depth":0, "max depth":1e5, + "coordinates":[[2.0e5,8e5],[2.0e5,0]], + "segments": + [ + {"length":2e5,"thickness":[5000],"angle":[60]} + ], + "temperature models":[{"model":"uniform","temperature":273}], + "composition models": + [ + {"model":"smooth", "compositions":[0,2], "operation":"add", "side distance fault center":2500, "center fractions":[1.5,1.5],"side fractions":[0.5,0.5]} + ] + + } + ] + +} \ No newline at end of file diff --git a/cookbooks/continental_compression_with_imposed_faults/doc/strain_rate_and_density_0_myr.svg b/cookbooks/continental_compression_with_imposed_faults/doc/strain_rate_and_density_0_myr.svg new file mode 100644 index 00000000000..75c3089e1c8 --- /dev/null +++ b/cookbooks/continental_compression_with_imposed_faults/doc/strain_rate_and_density_0_myr.svg @@ -0,0 +1 @@ + \ No newline at end of file diff --git a/cookbooks/continental_compression_with_imposed_faults/doc/strain_rate_and_density_35_myr.svg b/cookbooks/continental_compression_with_imposed_faults/doc/strain_rate_and_density_35_myr.svg new file mode 100644 index 00000000000..6a01065517a --- /dev/null +++ b/cookbooks/continental_compression_with_imposed_faults/doc/strain_rate_and_density_35_myr.svg @@ -0,0 +1 @@ + \ No newline at end of file diff --git a/cookbooks/continental_compression_with_imposed_faults/doc/strain_weakening_mechanism.part.prm b/cookbooks/continental_compression_with_imposed_faults/doc/strain_weakening_mechanism.part.prm new file mode 100644 index 00000000000..970a49c3716 --- /dev/null +++ b/cookbooks/continental_compression_with_imposed_faults/doc/strain_weakening_mechanism.part.prm @@ -0,0 +1,12 @@ +subsection Material model + subsection Visco Plastic + set Strain weakening mechanism = plastic weakening with plastic strain and viscous weakening with viscous strain + set Start plasticity strain weakening intervals = 0.5 + set End plasticity strain weakening intervals = 1.5 + set Cohesion strain weakening factors = 0.25 + set Friction strain weakening factors = 0.25 + set Start prefactor strain weakening intervals = 0.5 + set End prefactor strain weakening intervals = 1.5 + set Prefactor strain weakening factors = 0.1 + end +end diff --git a/doc/modules/changes/20250619_PrajaktaPMohite b/doc/modules/changes/20250619_PrajaktaPMohite new file mode 100644 index 00000000000..37ccdb5e589 --- /dev/null +++ b/doc/modules/changes/20250619_PrajaktaPMohite @@ -0,0 +1,4 @@ +New: A cookbook for continental compression that shows how to integrate faults into the initial conditions using the Geodynamic World Builder. + +
+(Prajakta Mohite, John Naliboff, 2025/06/19) diff --git a/doc/sphinx/references.bib b/doc/sphinx/references.bib index 350547732cd..2310dcf6ae2 100644 --- a/doc/sphinx/references.bib +++ b/doc/sphinx/references.bib @@ -12619,4 +12619,38 @@ @article{gleason:tullis:1995 pages={1--23}, year={1995}, publisher={Elsevier B.V.} +}@Article{Vasey:etal:2024, +AUTHOR = {Vasey, Dylan A. and Naliboff, John B. and Cowgill, Eric and Brune, Sascha and Glerum, Anne and Zwaan, Frank}, +TITLE = {Impact of rift history on the structural style of intracontinental rift-inversion orogens}, +JOURNAL = {Geology}, +VOLUME = {52}, +YEAR = {2024}, +NUMBER = {6}, +PAGES = {429--434}, +URL = {https://doi.org/10.1130/G51489.1}, +DOI = {10.1130/G51489.1} +} + +@Article{Zwaan:etal:2025, +AUTHOR = {Frank Zwaan and Sascha Brune and Anne C. Glerum and Dylan A. Vasey and John B. Naliboff and Gianreto Manatschal and Eric C. Gaucher}, +TITLE = {Rift-inversion orogens are potential hot spots for natural H2 generation }, +JOURNAL = {Science Advances}, +VOLUME = {11}, +YEAR = {2025}, +NUMBER = {8}, +PAGES = {eadr3418}, +URL = {https://www.science.org/doi/abs/10.1126/sciadv.adr3418}, +DOI = {10.1126/sciadv.adr3418} +} + +@Article{Howard:etal:2003, +AUTHOR = {Howard, J. and Cunningham, W. and Davies, S. and Dijkstra, Arjan and Badarch, G.}, +TITLE = {The stratigraphic and structural evolution of the Dzereg Basin, western Mongolia: Clastic sedimentation, transpressional faulting and basin destruction in an intraplate, intracontinental setting}, +JOURNAL = {Basin Research}, +VOLUME = {15}, +YEAR = {2003}, +NUMBER = {}, +PAGES = {45 - 72}, +URL = {}, +DOI = {10.1046/j.1365-2117.2003.00198.x} } \ No newline at end of file diff --git a/doc/sphinx/user/cookbooks/geophysical-setups.md b/doc/sphinx/user/cookbooks/geophysical-setups.md index 20c994f2385..c4c0fcdf525 100644 --- a/doc/sphinx/user/cookbooks/geophysical-setups.md +++ b/doc/sphinx/user/cookbooks/geophysical-setups.md @@ -103,6 +103,7 @@ cookbooks/multicomponent_steinberger/doc/multicomponent_steinberger.md cookbooks/morency_doin_2004/doc/morency_doin_2004.md cookbooks/crustal_deformation/doc/crustal_deformation.md cookbooks/continental_extension/doc/continental_extension.md +cookbooks/continental_compression_with_imposed_faults/doc/continental_compression_with_imposed_faults.md cookbooks/inner_core_convection/doc/inner_core_convection.md cookbooks/lower_crustal_flow/doc/lower_crustal_flow.md cookbooks/global_melt/doc/global_melt.md diff --git a/tests/continental_compression_with_imposed_faults.prm b/tests/continental_compression_with_imposed_faults.prm new file mode 100644 index 00000000000..8587b4091ff --- /dev/null +++ b/tests/continental_compression_with_imposed_faults.prm @@ -0,0 +1,27 @@ +# A test to check whether the continental compression with imposed faults cookbook +# at low resolutions for a single time step with one nonlinear iteration. + +include $ASPECT_SOURCE_DIR/cookbooks/continental_compression_with_imposed_faults/continental_compression_with_imposed_faults.prm + +#### Global parameters +set World builder file = $ASPECT_SOURCE_DIR/cookbooks/continental_compression_with_imposed_faults/continental_compression_with_imposed_faults.wb +set Start time = 0 +set End time = 0 +set Max nonlinear iterations = 1 +set Nonlinear solver tolerance = 1e-3 +set Output directory = continental_compression_with_imposed_faults + +subsection Solver parameters + subsection Stokes solver parameters + set Linear solver tolerance = 1e-3 + end +end +subsection Mesh refinement + set Initial adaptive refinement = 0 + set Initial global refinement = 0 +end + +# Post processing +subsection Postprocess + set List of postprocessors = composition statistics, basic statistics, velocity statistics, particles +end diff --git a/tests/continental_compression_with_imposed_faults/screen-output b/tests/continental_compression_with_imposed_faults/screen-output new file mode 100644 index 00000000000..8b1698f38aa --- /dev/null +++ b/tests/continental_compression_with_imposed_faults/screen-output @@ -0,0 +1,31 @@ + +Number of active cells: 400 (on 1 levels) +Number of degrees of freedom: 30,684 (3,362+441+1,681+3,600+3,600+3,600+3,600+3,600+3,600+3,600) + +Number of mesh deformation degrees of freedom: 882 + Solving mesh displacement system... 0 iterations. +*** Timestep 0: t=0 years, dt=0 years + Solving mesh displacement system... 0 iterations. + Solving temperature system... 0 iterations. + Advecting particles... done. + Rebuilding Stokes preconditioner... + Solving Stokes system (AMG)... 23+0 iterations. + Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.60695e-16, 0, 0, 0, 0, 0, 0, 0, 1 + Relative nonlinear residual (total system) after nonlinear iteration 1: 1 + + + WARNING: The nonlinear solver in the current timestep failed to converge. + Acting according to the parameter 'Nonlinear solver failure strategy': + Continuing to the next timestep even though solution is not fully converged. + + Postprocessing: + Compositions min/max/mass: 0/1.534/1.379e+10 // 0/0/0 // 0/1.493/3.306e+10 // 0/1/8e+09 // 0/1/8e+09 // 0/1/1.6e+10 // 0/1/1.28e+11 + RMS, max velocity: 0.000131 m/year, 0.000224 m/year + Writing particle output: output-continental_compression_with_imposed_faults/particles/particles-00000 + +Termination requested by criterion: end time + + + + +WARNING: During this computation 1 nonlinear solver failures occurred! diff --git a/tests/continental_compression_with_imposed_faults/statistics b/tests/continental_compression_with_imposed_faults/statistics new file mode 100644 index 00000000000..715cb281fe3 --- /dev/null +++ b/tests/continental_compression_with_imposed_faults/statistics @@ -0,0 +1,38 @@ +# 1: Time step number +# 2: Time (years) +# 3: Time step size (years) +# 4: Number of mesh cells +# 5: Number of Stokes degrees of freedom +# 6: Number of temperature degrees of freedom +# 7: Number of degrees of freedom for all compositions +# 8: Number of nonlinear iterations +# 9: Iterations for temperature solver +# 10: Iterations for Stokes solver +# 11: Velocity iterations in Stokes preconditioner +# 12: Schur complement iterations in Stokes preconditioner +# 13: Minimal value for composition plastic_strain +# 14: Maximal value for composition plastic_strain +# 15: Global mass for composition plastic_strain +# 16: Minimal value for composition noninitial_plastic_strain +# 17: Maximal value for composition noninitial_plastic_strain +# 18: Global mass for composition noninitial_plastic_strain +# 19: Minimal value for composition viscous_strain +# 20: Maximal value for composition viscous_strain +# 21: Global mass for composition viscous_strain +# 22: Minimal value for composition crust_upper +# 23: Maximal value for composition crust_upper +# 24: Global mass for composition crust_upper +# 25: Minimal value for composition crust_lower +# 26: Maximal value for composition crust_lower +# 27: Global mass for composition crust_lower +# 28: Minimal value for composition mantle_lithosphere +# 29: Maximal value for composition mantle_lithosphere +# 30: Global mass for composition mantle_lithosphere +# 31: Minimal value for composition asthenosphere +# 32: Maximal value for composition asthenosphere +# 33: Global mass for composition asthenosphere +# 34: RMS velocity (m/year) +# 35: Max. velocity (m/year) +# 36: Number of advected particles +# 37: Particle file name +0 0.000000000000e+00 0.000000000000e+00 400 3803 1681 25200 1 0 23 24 66 0.00000000e+00 1.53409020e+00 1.37942383e+10 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.49349236e+00 3.30589775e+10 0.00000000e+00 1.00000000e+00 8.00000000e+09 0.00000000e+00 1.00000000e+00 8.00000000e+09 0.00000000e+00 1.00000000e+00 1.60000000e+10 0.00000000e+00 1.00000000e+00 1.28000000e+11 1.31136528e-04 2.23606798e-04 19600 output-continental_compression_with_imposed_faults/particles/particles-00000