From c114e4330290e560bb28b4ce6476d329e47ec3eb Mon Sep 17 00:00:00 2001 From: naliboff Date: Sat, 19 Oct 2024 15:11:35 -0600 Subject: [PATCH] update continental extension cookbook --- .../continental_extension.prm | 202 +++++++++++++----- ..._scaling_continental_extension_cookbook.py | 106 +++++++++ 2 files changed, 249 insertions(+), 59 deletions(-) create mode 100644 cookbooks/continental_extension/viscous_prefactor_scaling_continental_extension_cookbook.py diff --git a/cookbooks/continental_extension/continental_extension.prm b/cookbooks/continental_extension/continental_extension.prm index 7f57d5f1d33..8a3c620e9ad 100644 --- a/cookbooks/continental_extension/continental_extension.prm +++ b/cookbooks/continental_extension/continental_extension.prm @@ -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 @@ -64,20 +66,19 @@ subsection Mesh refinement 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 @@ -86,7 +87,11 @@ end # 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 @@ -115,7 +120,7 @@ subsection Boundary velocity model 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 @@ -124,7 +129,22 @@ end # 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 @@ -132,14 +152,14 @@ end # 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); @@ -149,7 +169,7 @@ 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 + set List of model names = initial composition end # Temperature boundary conditions @@ -159,12 +179,7 @@ end # 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 @@ -220,25 +235,29 @@ end # 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 @@ -247,10 +266,18 @@ subsection Material model 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 + # 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 @@ -264,20 +291,38 @@ subsection Material model # 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 @@ -290,6 +335,17 @@ subsection Material model 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 @@ -302,9 +358,26 @@ subsection Gravity model 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 @@ -316,4 +389,15 @@ subsection Postprocess 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 diff --git a/cookbooks/continental_extension/viscous_prefactor_scaling_continental_extension_cookbook.py b/cookbooks/continental_extension/viscous_prefactor_scaling_continental_extension_cookbook.py new file mode 100644 index 00000000000..d28edf76ae0 --- /dev/null +++ b/cookbooks/continental_extension/viscous_prefactor_scaling_continental_extension_cookbook.py @@ -0,0 +1,106 @@ +# This script scales the prefactors terms from viscous flow laws commonly used +# in geodynamic modeling in a manner that is internally consistent with ASPECT's +# conventions following Dannberg et al. (2017, https://doi.org/10.1002/2017GC006944). + +import numpy as np + + +# 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 + +pfct_gt95_dqtz_disl_rprt = 1.1e-4 # Original reported prefactor in units of MPa^-n micrometers^m + +sexp_gt95_dqtz_disl_rprt = 4.0 # Original reported stress exponent + +gexp_gt95_dqtz_disl_rprt = 0 # Original reported grain size exponent + +# Convert prefactor to units of Pa^-n meters^m s^-1 + +pfct_gt95_dqtz_disl_unit = (pfct_gt95_dqtz_disl_rprt) * \ + (1.e-6**sexp_gt95_dqtz_disl_rprt) * \ + (1.e-6**gexp_gt95_dqtz_disl_rprt) + +# Calculate prefactor scaling term + +pfct_gt95_dqtz_disl_sfac = 2.**(sexp_gt95_dqtz_disl_rprt-1.) * 3.**((sexp_gt95_dqtz_disl_rprt+1.)/2.) + +# Calculated modified prefactor term for ASPECT + +pfct_gt95_dqtz_disl_aspt = pfct_gt95_dqtz_disl_unit * pfct_gt95_dqtz_disl_sfac + +print('') + +print('Gleason and Tullis (1995) Quarzite Dislocation Prefactor Scaling = ', pfct_gt95_dqtz_disl_sfac) + +print('Gleason and Tullis (1995) Quarzite Dislocation Prefactor ASPECT = ', pfct_gt95_dqtz_disl_aspt) + +print('') + + +# Rybacki and Dresen (2000), J. Geophys. Res., v.111(B3). +# "Dislocation and diffusion creep of synthetic anorthite aggregates". +# https://doi.org/10.1029/2000JB900223 + +pfct_ry00_want_disl_rprt = 10.**(2.6) # Original reported log(prefactor) in units of MPa^(-n-r) micrometers^m + +sexp_ry00_want_disl_rprt = 3.0 # Original reported stress exponent + +gexp_ry00_want_disl_rprt = 0 # Original reported grain size exponent + +wfug_ry00_want_disl_rprt = 0 # Original reported water fugacity + +# Convert prefactor to units of Pa^(-n-r) meters^m s^-1 + +pfct_ry00_want_disl_unit = (pfct_ry00_want_disl_rprt) * \ + (1.e-6**(sexp_ry00_want_disl_rprt + wfug_ry00_want_disl_rprt)) * \ + (1.e-6**gexp_ry00_want_disl_rprt) + +# Calculate prefactor scaling term + +pfct_ry00_want_disl_sfac = 2.**(sexp_ry00_want_disl_rprt-1.) * 3.**((sexp_ry00_want_disl_rprt+1.)/2.) + +# Calculated modified prefactor term for ASPECT + +pfct_ry00_want_disl_aspt = pfct_ry00_want_disl_unit * pfct_ry00_want_disl_sfac + +print('') + +print('Rybacki et al. (2000) Wet Anorthite Dislocation Prefactor Scaling = ', pfct_ry00_want_disl_sfac) + +print('Rybacki et al. (2000) Wet Anorthite Dislocation Prefactor ASPECT = ', pfct_ry00_want_disl_aspt) + +print('') + + +# 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". +# https://doi.org/10.1029/138GM06 + +pfct_hk04_doli_disl_rprt = 1.10e5 # Original reported prefactor in units of MPa^-n micrometers^m + +sexp_hk04_doli_disl_rprt = 3.5 # Original reported stress exponent + +gexp_hk04_doli_disl_rprt = 0 # Original reported grain size exponent + +# Convert prefactor to units of Pa^-n meters^m s^-1 + +pfct_hk04_doli_disl_unit = (pfct_hk04_doli_disl_rprt) * \ + (1.e-6**sexp_hk04_doli_disl_rprt) * \ + (1.e-6**gexp_hk04_doli_disl_rprt) + +# Calculate prefactor scaling term + +pfct_hk04_doli_disl_sfac = 2.**(sexp_hk04_doli_disl_rprt-1.) * 3.**((sexp_hk04_doli_disl_rprt+1.)/2.) + +# Calculated modified prefactor term for ASPECT + +pfct_hk04_doli_disl_aspt = pfct_hk04_doli_disl_unit * pfct_hk04_doli_disl_sfac + +print('') + +print('Hirth Kohlsted (2004) Dry Olivine Dislocation Prefactor Scaling = ', pfct_hk04_doli_disl_sfac) + +print('Hirth Kohlsted (2004) Dry Olivine Dislocation Prefactor ASPECT = ', pfct_hk04_doli_disl_aspt) + +print('')