Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[WIP] update continental extension cookbook #6117

Open
wants to merge 1 commit into
base: main
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
202 changes: 143 additions & 59 deletions cookbooks/continental_extension/continental_extension.prm
Original file line number Diff line number Diff line change
Expand Up @@ -2,33 +2,35 @@
# This cookbook is based off numerous published studies, five 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
# Rift migration explains continental margin asymmetry and crustal hyperextension.
# https://doi.org/10.1038/ncomms5014.
# 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
# Depth-dependent extension, two-stage breakup and cratonic underplating at rifted margins.
# https://doi.org/10.1038/nature09988.
# 3. Naliboff, J., and Buiter, S.H. (2015), Earth Planet. Sci. Lett., v.421, p.58-67,
# "Rift Reactivation and migration during multiphase extension"
# "Rift Reactivation and migration during multiphase extension".
# https://doi.org/10.1016/j.epsl.2015.03.050.
# 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"
# network evolution". https://doi.org/10.1029/2019GL086611.
# 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"
# detachment faults: Solid-block rotation and apparent unbending".
# https://doi.org/10.1029/2021GC009681

#### Global parameters
set Dimension = 2
set Start time = 0
set End time = 5e6
set Use years in output instead of seconds = true
set Nonlinear solver scheme = single Advection, iterated defect correction Stokes
set Nonlinear solver scheme = iterated Advection and Stokes
set Nonlinear solver tolerance = 1e-4
set Max nonlinear iterations = 100
set Max nonlinear iterations = 50
set CFL number = 0.5
set Maximum time step = 20e3
set Maximum time step = 10e3
set Output directory = output-continental_extension
set Pressure normalization = no

#### Parameters describing the model

# Governing equations
subsection Formulation
set Formulation = Boussinesq approximation
Expand Down Expand Up @@ -64,20 +66,19 @@
end
end

# Use the Eisenstat Walker method to automatically determine the
# linear solver tolerance during the defect Picard iterations.
# Adjusting the Maximum linear solver tolerance will affect how long
# it takes to reach a solution, but not the actual value of the
# solution.
# For models with a nonlinear visco-plastic rheology
# and a minimum-maximum viscosity contrast greater
# 4 orders of magnitude, the AMG solver combined
# with a large number of cheap Stokes solver steps
# produces the lowest Stokes solver times.
subsection Solver parameters
subsection Stokes solver parameters
set Number of cheap Stokes solver steps = 0
set Stokes solver type = block AMG
end
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

subsection Newton solver parameters
set Maximum linear Stokes solver tolerance = 1e-2
set Use Eisenstat Walker method for Picard iterations = true
end
end

Expand All @@ -86,7 +87,11 @@
# 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 Mesh deformation boundary indicators = top: free surface, top: diffusion

# Allow the mesh on the left and right boundaries to deform in order
# order to prevent significant mesh distortion at the top corners.
set Additional tangential mesh velocity boundary indicators = left, right

subsection Free surface
set Surface velocity projection = normal
Expand Down Expand Up @@ -115,7 +120,7 @@
end
end

# Number and names of compositional fields
# Number and names of compositional fields, which are tracked with particles.
# The five compositional fields represent:
# 1. The plastic strain that accumulates over time, with the initial plastic strain removed
# 2. The plastic strain that accumulated over time, including the initial plastic strain values
Expand All @@ -124,22 +129,37 @@
# 5. The mantle lithosphere
subsection Compositional fields
set Number of fields = 5
set Names of fields = noninitial_plastic_strain, plastic_strain, crust_upper, crust_lower, mantle_lithosphere
set Names of fields = plastic_strain, \
noninitial_plastic_strain, \
crust_upper, \
crust_lower, \
mantle_lithosphere
set Types of fields = strain, \
strain, \
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, \
crust_upper: initial crust_upper, \
crust_lower: initial crust_lower, \
mantle_lithosphere: initial mantle_lithosphere
end

# Initial values of different compositional fields
# The upper crust (20 km thick), lower crust (20 km thick)
# and mantle (60 km thick) are continuous horizontal layers
# of constant thickness. The non initial plastic strain is set
# to 0 and the initial plastic strain is randomized between
# 0.5 and 1.5.
# 0.5 and 1.5 in a repeatable manner using the rand_seed() function.
subsection Initial composition model
set Model name = function

subsection Function
set Variable names = x,y
set Function expression = 0; \
if(x>50.e3 && x<150.e3 && y>50.e3, 0.5 + rand_seed(1), 0); \
set Function expression = if(x>50.e3 && x<150.e3 && y>50.e3, 0.5 + rand_seed(1), 0); \
0; \
if(y>=80.e3, 1, 0); \
if(y<80.e3 && y>=60.e3, 1, 0); \
if(y<60.e3, 1, 0);
Expand All @@ -149,7 +169,7 @@
# 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
set List of model names = initial composition
end

# Temperature boundary conditions
Expand All @@ -159,12 +179,7 @@
# 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
set List of model names = initial temperature
end

# Initial temperature field
Expand Down Expand Up @@ -220,25 +235,29 @@
# 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.
# creep flow-laws, with distinct values for the upper crust (quartzite), lower
# crust (wet anorthite) and mantle (dry olivine).
subsection Material model
set Model name = visco plastic

# Average the viscosity over the element to reduce complexity
# and improve solver performance.
set Material averaging = none

subsection Visco Plastic
# Reference temperature and viscosity

# Reference temperature used for calculating densities
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 Minimum strain rate = 1.e-22
set Reference strain rate = 1.e-16

# Limit the viscosity with minimum and maximum values
set Minimum viscosity = 1e18
set Minimum viscosity = 1e20
set Maximum viscosity = 1e26

# Thermal diffusivity is adjusted to match thermal conductivities
Expand All @@ -247,10 +266,18 @@
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
# Below, the density of each field classified as a chemical
# compositon (i.e., lithology) in the Compositional fields

Check warning on line 270 in cookbooks/continental_extension/continental_extension.prm

View workflow job for this annotation

GitHub Actions / Check for new typos

perhaps "compositon" should be "composition".
# subsection is assigned using the name of the field. This
# approach allows excluding values for fields not classified
# as chemical compositions, which are automatically excluded
# during volume fraction calculations of material properties.
set Densities = background: 3300, \
crust_upper: 2700, \
crust_lower: 2900, \
mantle_lithosphere: 3300

set Thermal expansivities = 2e-5

# Harmonic viscosity averaging
set Viscosity averaging scheme = harmonic
Expand All @@ -264,20 +291,38 @@
# 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
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, 0.0, 18.e-6
# "Rheology of the upper mantle and the mantle wedge:a view from the experimentalists".
# https://doi.org/10.1029/138GM06
# 2. Upper crust (quartzite)
# Gleason and Tullis (1995), Tectonophysics, v.247 p.1-23.
# "A flow law for dislocation creep of quartz aggregates determined with the molten salt cell"
# https://doi.org/10.1016/0040-1951(95)00011-B
# 3. Lower crust (wet anorthite)
# Rybacki and Dresen (2000), J. Geophys. Res., v.111(B3).
# "Dislocation and diffusion creep of synthetic anorthite aggregates".
# https://doi.org/10.1029/2000JB900223
# Note that the viscous pre-factors below are scaled to plane strain from unixial strain experiments,
# which is performed in the script viscous_prefactor_scaling_continental_extension_cookbook.py located
# within this cookbook folder.
set Prefactors for dislocation creep = background: 7.3704e-15, \
crust_upper: 1.3718e-26, \
crust_lower: 1.4332e-14, \
mantle_lithosphere: 7.3704e-15

set Stress exponents for dislocation creep = background: 3.5, \
crust_upper: 4.0, \
crust_lower: 3.0, \
mantle_lithosphere: 3.5

set Activation energies for dislocation creep = background: 530.e3, \
crust_upper: 223.e3, \
crust_lower: 345.e3, \
mantle_lithosphere: 530.e3

set Activation volumes for dislocation creep = background: 18.e-6, \
crust_upper: 0, \
crust_lower: 0, \
mantle_lithosphere: 18.e-6

# Plasticity parameters
set Angles of internal friction = 30
Expand All @@ -290,6 +335,17 @@
set End plasticity strain weakening intervals = 1.5
set Cohesion strain weakening factors = 0.25
set Friction strain weakening factors = 0.25

# Following Duretz et al. (2021), a plastic damper is used to regularize
# the Drucker Prager plasticity formulation, which provides stable shear
# band widths for a given plastic damper viscosity, cell size, and strain
# rate. Note that the cell size of the current setup (2.5-1.25 km) is
# larger than the shear band width given by a plastic damper viscosity
# of 1e21, which is on the order of hundreds of meters. This can be
# verified by increasing the number of global refinements from 3 to 5 or 6.
set Use plastic damper = true
set Plastic damper viscosity = 1e21

end
end

Expand All @@ -302,9 +358,26 @@
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 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
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
set List of postprocessors = basic statistics, composition statistics, heat flux statistics, \
mass flux statistics, particles, pressure statistics, \
temperature statistics, topography, velocity statistics, visualization

subsection Visualization
set List of output variables = material properties, heat flux map, named additional outputs, strain rate
Expand All @@ -316,4 +389,15 @@
set List of material properties = density, viscosity
end
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
Loading
Loading