#### ModiContinental Extension Cookbook # This PRM is a modified version of the continetnal extension cookbook, which is based off numerous published studies # and includes the five listed below. # For additional information, see these publications and references therein. # 1. Brune, S., Heine, C., Perez-Gussinye, M., and Sobolev, S.V. (2014), Nat. Comm., v.5, n.4014, # Rift migration explains continental margin asymmetry and crustal hyperextension # 2. Huismans, R., and Beaumont, C. (2011), Nature, v.473, p.71-75. # Depth-dependent extension, two-stage breakup and cratonic underplating at rifted margins # 3. Naliboff, J., and Buiter, S.H. (2015), Earth Planet. Sci. Lett., v.421, p.58-67, # "Rift Reactivation and migration during multiphase extension" # 4. Naliboff, J., Glerum, A., Sascha, S., Peron-Pinvidic, G., and Wrona, T. (2020), Geophys. # Res. Lett., 47, e2019GL086611, "Development of 3‐D rift heterogeneity through fault # network evolution" # 5. Sandiford, D., Brune, S., Glerum, A., Naliboff, J., and Whittaker, J.M. (2021), Geophys. # Geochem. Geosys., 22, e2021GC009681, "Kinematics of footwall exhumation at oceanic # detachment faults: Solid-block rotation and apparent unbending" #### Global parameters set Dimension = 2 set Start time = 0 set End time = 10e6 set Use years in output instead of seconds = true set Nonlinear solver scheme = single Advection, iterated defect correction Stokes set Nonlinear solver tolerance = 1e-5 set Max nonlinear iterations = 200 set CFL number = 0.25 set Maximum time step = 10e3 set Output directory = output-continental_extension_2D_particles set Pressure normalization = no set Use operator splitting = false #### Parameters describing the model # Governing equations subsection Formulation set Formulation = Boussinesq approximation end # Q2Q1 Element for Stokes, Q2 for temperature, Q2 for composition subsection Discretization set Composition polynomial degree = 2 set Stokes velocity polynomial degree = 2 set Temperature polynomial degree = 2 end # Model geometry (200x100 km, 20 km spacing) subsection Geometry model set Model name = box subsection Box set X repetitions = 10 set Y repetitions = 5 set X extent = 200e3 set Y extent = 100e3 end end # Globally refine the mesh to 1.25 km spacing, and then # adaptively refine the mesh to 0.625 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 = 4 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>=50e3 && x>=40.e3 && x<=160.e3, 5, 4) end end # Solver parameters subsection Solver parameters subsection Stokes solver parameters set Stokes solver type = block GMG set Number of cheap Stokes solver steps = 1000 set Linear solver tolerance = 1e-7 set GMRES solver restart length = 100 set Use full A block as preconditioner = true end subsection Newton solver parameters set Maximum linear Stokes solver tolerance = 1e-3 set Use Eisenstat Walker method for Picard iterations = true end subsection Advection solver parameters set GMRES solver restart length = 50 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 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 # The outward velocity (x-direction) on the left and right walls is 0.25 cm/year # The vertical velocity at the base is 0.25 cm/year (balances outflow on sides) # 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.0025, w=200.e3, d=100.e3 set Function expression = if (x < w/2 , -v, v) ; v*2*d/w end end # Number and name of compositional fields # The eleven compositional fields represent: # 1. The plastic strain that accumualates over time, with the initial plastic strain removed # 2. The plastic strain that accumulated over time, including the initial plastic strain values # 3. The upper crust # 4. The lower crust # 5. The mantle lithosphere subsection Compositional fields set Number of fields = 5 set Names of fields = plastic_strain, \ noninitial_plastic_strain, \ crust_upper, \ crust_lower, \ mantle_lithosphere set Compositional field methods = particles set Mapped particle properties = plastic_strain: plastic_strain, \ noninitial_plastic_strain: noninitial_plastic_strain, \ crust_upper: initial crust_upper, \ crust_lower: initial crust_lower, \ mantle_lithosphere: initial mantle_lithosphere end # Spatial domain of different compositional fields subsection Initial composition model set List of model names = ascii data subsection Ascii data model set Data directory = set Data file name = composition.txt end end # Composition: fixed on bottom (inflow boundary), free on sides and top subsection Boundary composition model set Fixed composition boundary indicators = bottom set List of model names = initial composition 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). subsection Boundary temperature model set Fixed temperature boundary indicators = bottom, top set List of model names = box subsection Box set Bottom temperature = 1613 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.055 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; # Thermal Conductivity - all layers = 2.5 (W/(m K)); # To satisfy these constraints, the following values are required: # Surface Temperature - lower crust (ts2) = 633 K # - mantle (ts3) = 893 K # Surface Heat Flow - lower crust (qs2) = 0.035 W/m^2; # - mantle (qs3) = 0.030 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=100e3, ts1=273, ts2=633, ts3=893, \ A1=1.e-6, A2=0.25e-6, A3=0.0, \ k1=2.5, k2=2.5, k3=2.5, \ qs1=0.055, qs2=0.035, qs3=0.030 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), \ ts3 + (qs3/k3)*(h-y-40.e3) - (A3*(h-y-40.e3)*(h-y-40.e3))/(2.0*k3) ) ); 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 subsection Compositional heating set Use compositional field for heat production averaging = 1, 0, 0, 1, 1, 1 set Compositional heating values = 0., 0., 0., 1.0e-6, 0.25e-6, 0.0 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). Table 1 of Naliboff and Buiter (2015), # Earth Planet. Sci. Lett., v.421, p. 58-67 contains values for each of these flow laws. 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. # Density values of 1 are assigned to "strain" fields, which are not taken into # account when computing material properties. set Densities = 3300, 1.0, 1.0, 2700, 2900, 3300 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 # 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 (wet anorthite) # Rybacki and Dresen (2000), J. Geophys. Res., 105( B11), p.26017– 26036. # "Dislocation and diffusion creep of synthetic 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 1.0 set Prefactors for dislocation creep = 7.37e-15, 1.0, 1.0, 1.37e-26, 1.43e-14, 7.37e-15 set Stress exponents for dislocation creep = 3.5, 1.0, 1.0, 4.0, 3.0, 3.5 set Activation energies for dislocation creep = 530.e3, 0.0, 0.0, 223.e3, 356.e3, 530.e3 set Activation volumes for dislocation creep = 18.e-6, 0.0, 0.0, 0.0, 0.0, 18.e-6 # 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. set Strain weakening mechanism = plastic weakening with plastic strain only 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 end end # Gravity model subsection Gravity model set Model name = vertical subsection Vertical set Magnitude = 9.81 end end # Post processing subsection Postprocess set List of postprocessors = basic statistics, composition statistics, heat flux densities, heat flux statistics, mass flux statistics, matrix statistics, particles, pressure statistics, temperature statistics, topography, velocity statistics, visualization subsection Visualization set List of output variables = density, heat flux map, named additional outputs, strain rate, viscosity set Output format = vtu set Time between graphical output = 10.e3 set Interpolate output = true set Write higher order output = true set Output format = vtu end subsection Particles set Minimum particles per cell = 40 set Maximum particles per cell = 60 set Load balancing strategy = remove and add particles set List of particle properties = initial composition, viscoplastic strain invariants, pT path, position set Interpolation scheme = quadratic 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 end end subsection Interpolator subsection Quadratic least squares set Use boundary extrapolation = false set Use quadratic least squares limiter = true end end set Time between data output = 10.e3 set Data output format = vtu end end # Checkpointing subsection Checkpointing set Steps between checkpoint = 10 end