############################################################## ################# Oblique Extension Models ################### ############################################################## #################### Global parameters ####################### set Dimension = 3 set Start time = 0 set End time = 10e6 set Use years in output instead of seconds = true set Nonlinear solver scheme = single Advection, single Stokes set Nonlinear solver tolerance = 1e-7 set Max nonlinear iterations = 10 set CFL number = 0.3 set Output directory = output-frontback_test_03 set Timing output frequency = 1 set Pressure normalization = no set Resume computation = auto ##################### Checkpointing ########################## subsection Checkpointing set Steps between checkpoint = 5 set Time between checkpoint = 0 end ############# Parameters describing the model ################# # Model geometry (800x1200x100 km, 3.125 km spacing) subsection Geometry model set Model name = box subsection Box set Y periodic = true set X repetitions = 3 set Y repetitions = 4 set Z repetitions = 1 set X extent = 300e3 set Y extent = 400e3 set Z extent = 100e3 end end # Mesh refinement specifications subsection Mesh refinement set Initial adaptive refinement = 0 set Initial global refinement = 5 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>=0, 5, 5) 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 1 cm/year # The vertical velocity at the base is balaced by the Boundary traction model (balances outflow on sides) subsection Boundary velocity model set Prescribed velocity boundary indicators = left:function, right:function #, bottom:function #, front yz:function, back yz:function #, bottom z:function # set Tangential velocity boundary indicators = front, back subsection Function set Variable names = x,y, z set Function constants = cm=0.01, year=1 set Function expression = if (x<100e3, -0.7071067812*cm/year, 0.7071067812*cm/year); if (x<100e3, -0.7071067812*cm/year, 0.7071067812*cm/year); 0 end end subsection Boundary traction model set Prescribed traction boundary indicators = bottom:initial lithostatic pressure subsection Initial lithostatic pressure set Representative point = 200e3, 0e3, 0e3 end end # Number and names of compositional fields # The four compositional fields represent the upper crust, lower crust, mantle # and a 'seed' placed to help localize deformation. subsection Compositional fields set Number of fields = 4 set Names of fields = upper, lower, mantle, plastic_strain end # Spatial domain of different compositional fields # The upper crust, lower crust and mantle are continuous horizontal layers # of varying thickness. The seed is implaced into the top 30 km. subsection Initial composition model set Model name = function subsection Function set Variable names = x,y, z set Function constants = scar_wid=20e3, ylim1=50e3, ylim2=210.e3, ylim3=350e3, xlim1=100e3, xlim3=180e3 set Function expression = if(z>=80.e3, 1, 0); \ if(z<80.e3 && z>=70.e3, 1, 0); \ if(z<70.e3 && z>-100.e3,1, 0); \ if((z>60.e3 && (y>=ylim1 && y=(xlim3)))) || \ (z>60.e3 && x>=xlim1 && x=(ylim2-scar_wid))) || \ (z>60.e3 && y>=ylim2 && y=(xlim1))), 0.5 + rand_seed(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 List of model names = box subsection Box set Bottom temperature = 1573 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: # Layer Surface Temperature - upper crust (ts1) = 273 K; # mantle (ts3) = 823 K; # Model Base Temperature - (tb) = 1573 K; # Heat Production - upper crust (A) = 1.5e-6 W/m^3; # Thermal Conductivity - upper crust (k1) = 2.5 (W/(m K)); # lower crust (k2) = 2.5 (W/(m K)); # mantle (k3) = 3.3 (W/(m K)); # To satisfy these constraints, the following values are required: # Layer Surface Heat Flow - upper crust (qs1) = 0.065357 W/m^2; # lower crust (qs2) = 0.035357 W/m^2; # mantle (qs3) = 0.035357 W/m^2; # Temperature - base of upper crust (ts2) = 681.5714 subsection Initial temperature model set Model name = function subsection Function set Variable names = x,y, z set Function constants = h=100e3,ts1=273,ts2=681.5714,ts3=823., \ k1=2.5,k2=2.5,k3=3.3,A=1.5e-6, \ qs1=0.0653571,qs2=0.035357,qs3=0.035357,qb3=0.035357 set Function expression = if( (h-z)<=20.e3, \ ts1 + (qs1/k1)*(h-z) - (A*(h-z)*(h-z))/(2.0*k1), \ if((h-z)>20.e3 && (h-z)<=30.e3, ts2 + (qs2/k2)*(h-z-20.e3),\ ts3 + (qs3/k3)*(h-z-30.e3))); 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 # background, upper, lower, mantle, seed # set Use compositional field for heat production averaging = 1, 1, 1, 1, 0.0 set Compositional heating values = 0.0, 1.0e-6, 0.25e-6, 0.0, 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 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. # background, upper, lower, mantle, seed set Densities = 3300, 2700, 2900, 3300, 1.0 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) # 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) # 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 # background, upper, lower, mantle, seed set Prefactors for dislocation creep = 6.52e-16, 8.57e-28, 7.13e-18, 6.52e-16, 1.00e-50 set Stress exponents for dislocation creep = 3.5, 4.0, 3.0, 3.5, 1.0 set Activation energies for dislocation creep = 530.e3, 223.e3, 345.e3, 530.e3, 0.0 set Activation volumes for dislocation creep = 18.e-6, 0.0, 0.0, 18.e-6, 0.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. 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 # Remove nullspace in y direction subsection Nullspace removal set Remove nullspace = net y translation end ######################## Postprocessing ######################## # Outputs every 500,000 years subsection Postprocess set List of postprocessors = velocity statistics, basic statistics, temperature statistics, visualization subsection Visualization set List of output variables = density, viscosity, strain rate, stress, shear stress set Output format = vtu set Time between graphical output = 0.5e6 end end subsection Solver parameters subsection Stokes solver parameters set Number of cheap Stokes solver steps = 0 end end