#### Continental Extension Cookbook # This cookbook is based off numerous published studies, three of which are 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" #### Global parameters set Dimension = 3 set Start time = 0 set End time = 10e6 set Resume computation = auto set Use years in output instead of seconds = true set Nonlinear solver scheme = single Advection, iterated Newton Stokes set Nonlinear solver tolerance = 1e-4 set Max nonlinear iterations = 100 # Normally set to 50 or 100 set CFL number = 0.5 set Maximum time step = 50e3 # Normally set to 20 Kyr set Output directory = output-continental_extension set Timing output frequency = 1 set Pressure normalization = no # 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 Z repetitions = 5 set X extent = 200e3 set Y extent = 100e3 set Z extent = 100e3 end end # Adaptively refine the mesh to 2.5 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 global refinement = 2 set Initial adaptive refinement = 1 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,z set Function expression = if ( z>=50e3 && x>=60.e3 && x<=140.e3, 3, 2) end set Run postprocessors on initial refinement = false end # Quadratic elements for compostion, temperature, and Stokes velocity subsection Discretization set Composition polynomial degree = 2 set Stokes velocity polynomial degree = 2 set Temperature polynomial degree = 2 end # Solver parameters subsection Solver parameters subsection Stokes solver parameters set Stokes solver type = block AMG set Number of cheap Stokes solver steps = 0 end subsection Newton solver parameters set Max Newton line search iterations = 5 set Max pre-Newton nonlinear iterations = 10 set Maximum linear Stokes solver tolerance = 1e-1 set Nonlinear Newton solver switch tolerance = 1e-4 set SPD safety factor = 0.9 set Stabilization preconditioner = SPD set Stabilization velocity block = SPD set Use Newton failsafe = false set Use Newton residual scaling method = false set Use Eisenstat Walker method for Picard iterations = true end end # A normal free surface advection may become unstable when significant topography develops subsection Mesh deformation set Mesh deformation boundary indicators = top: free surface subsection Free surface set Surface velocity projection = normal end end # Velocity on boundaries characterized by functions # The outward velocity (x-direction) on the left and right walls is 0.65 cm/year. # The vertical velocity at the base is 0.25 cm/year (balances outflow on sides). #the obliquity angle set by the parameter "o" (o degrees initially fro 2 myrs then 30 degrees for until 4 myr) # 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 z:function set Tangential velocity boundary indicators = front, back subsection Function set Variable names = x,y,z set Function constants = v=0.0065,o=30., c=100.e3, w=200.e3, d=100.e3 set Function expression = if (x60.e3 && x<140.e3 && z>50.e3, 0.5 + rand_seed(1), 0); \ if(z>=80.e3, 1, 0); \ if(z<80.e3 && z>=60.e3, 1, 0); \ if(z<60.e3 && z>=0.e3, 1, 0); end end # Composition: fixed on bottom, 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 Model name = initial temperature 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.40e-6 W/m^3; # Heat Production - lower crust (A2) = 1.10e-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; # Initial temperature field subsection Initial temperature model set List of model names = function subsection Function set Variable names = x,y,z set Function constants = h=100.e3,ts1=273,ts2=633,ts3=893, \ A1=1.40e-6,A2=1.10e-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-z)<=20.e3, \ ts1 + (qs1/k1)*(h-z) - (A1*(h-z)*(h-z))/(2.0*k1), \ if( (h-y)>20.e3 && (h-z)<=40.e3, \ ts2 + (qs2/k2)*(h-z-20.e3) - (A2*(h-z-20.e3)*(h-z-20.e3))/(2.0*k2), \ ts3 + (qs3/k3)*(h-z-40.e3) - (A3*(h-z-40.e3)*(h-z-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 Compositional heating values = 0.0, 0.0, 0.0, 1.40e-6, 1.10e-6, 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 subsection Visco Plastic # Reference temperature and viscosity set Reference temperature = 273 set Reference viscosity = 1e22 # 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 set Define thermal conductivities = true set Thermal conductivities = 2.5 set Heat capacities = 875 set Densities = 3300, 1.0, 1.0, 2700, 2850, 3300 set Thermal expansivities = 3.1e-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) # Rutter & Brodie (2004), J. Struct. Geol., v.26, p.2011-2023. # "Experimental grain size-sensitive flow of hot-pressed Brazilian quartz aggregates" # 3. Lower crust and weak seed (wet anorthite) # Gleason and Tullis (1995), Tectonophysics, v.472, p.213–225. # " A flow law for dislocation creep of quartz aggregates determined with the molten salt cell" # Note that the viscous pre-factors below are scaled to plane strain from unixial strain experiments. set Prefactors for dislocation creep = 6.52e-16, 1.00e-50, 1.00e-50, 8.57e-28, 7.13e-18, 6.52e-16 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, 345.e3, 530.e3 set Activation volumes for dislocation creep = 18.e-6, 0.0, 0.0, 0., 0., 18.e-6 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, 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 Time between graphical output = 100.e3 set Interpolate output = false set Output format = vtu end end subsection Checkpointing set Steps between checkpoint = 50 end