diff --git a/cookbooks/glacial_cycling_melt_transport/doc/glacial_cycling_melt_transport.md b/cookbooks/glacial_cycling_melt_transport/doc/glacial_cycling_melt_transport.md new file mode 100644 index 00000000000..cc3ff70bc9d --- /dev/null +++ b/cookbooks/glacial_cycling_melt_transport/doc/glacial_cycling_melt_transport.md @@ -0,0 +1,33 @@ +(sec:cookbooks:glacial-cycling-melt-trasport)= +# Glacial Unloading and Melt Migration with Viscoplastic Rheologies + +*This section was contributed by Prajakta Mohite & John Naliboff.* + +This cookbook demonstrates interactions between time and spatially-dependent surface loading, lithospheric deformation, and reactive melt transport. It is motivated by observations and prior modeling investigations highlighting the role of glacial loading and unloading on volcanic activity, and designed to serve as a potential template for investigating the effects of glacial cycles on tectonic and magmatic processes, such as in Antarctica and Iceland. +For example, the model in this cookbook can be used or extended to investigate the following topics: + +1. Quantify how changes in surface loading drive stress evolution in the crust and upper mantle. + +2. Examine the onset, distribution, and temporal evolution of melt during and after rapid unloading. + +3. Compare the relative effects of different rheological laws and material properties on crustal and mantle dynamics. + +## The input file + +The model in this cookbook is designed to investigate the mechanical and magmatic response of the continental lithosphere and asthenosphere to rapid deglaciation, using a physically motivated and time-dependent parabolic ice load. The model domain is 400 (X) x 200 (Y) km wide, and contains layers representing the continental crust, mantle lithosphere, and asthenosphere. The surface load is imposed as a 200 km wide parabolic ice sheet with a maximum height of 1 km, which is centered along the horizontal extent of the model and applied through a traction boundary condition. The load remains constant for a prescribed period of 1 Myr, after which it decreases linearly to zero over an interval of 20 kyr, thereby simulating rapid deglaciation. + +This is defined in the parameter file as: +```{literalinclude} traction_boundary.part.prm +``` +where ```{math}/rho``` , ``` /g```, and thickness represent ice density, gravitational acceleration, and ice thickness, respectively, and the spatial function defines the parabolic load geometry. Respectively, the model side and bottom boundaries contain no-slip and free-slip velocity boundary conditions. + +The highly nonlinear system of equations is solved with the nonlinear solver scheme 'iterated Advection and Stokes', following prior investigations and cookbooks {ref}`sec:cookbooks:global-melt` modeling two-phase reactive melt transport in ASPECT {cite:t}`dannberg:heister:2016`. Respectively, the linear and nonlinear solver tolerance are set to, although we note that the model dynamics may slightly or moderately differ when using stricter values. Here, the values are selected as a compromise between accuracy and model run times. We note that the introduction of a free surface, as compared to a free-slip top boundary, introduces significant challenges for both the linear and nonlinear solvers. Likewise, we note that using an open traction boundary at the model base with the traction magnitude equal to the adiabatic pressure may also introduce significant nonlinear solver issues. + +The initial temperature field is defined using a depth-dependent conductive temperature profile through the lithosphere in a similar fashion to the continental extension cookbook, and an adiabatic gradient of 0.5 C/km from the base of the lithosphere (80 km depth) to the model base (200 km depth). Although temperatures are specified for the model sides, these values are not used, as the boundaries are insulating (zero net heat flux). Temperature evolves in the model domain through flow induced by spatiotemporal variations in surface loading, and includes the effects of adiabatic and shear heating introduced through the use of the extended Boussinesq approximation [cookbooks/continental_extension/](https://github.com/geodynamics/aspect/tree/main/cookbooks/continental_extension). +The model rheology follows nonlinear dislocation creep similar to the setup in the continental extension cookbook, with distinct flow parameters for the upper crust, lower crust, and mantle. Melt generation is parameterized using the {cite:t}`katz:etal:2003` model, with a prescribed melt freezing rate of 0.5 and maximum melt migration depth of 30 km. + +The configuration below allows us to combine two-phase flow transport (Katz 2003) with viscoplastic rheologies[reaction-model:katz2003_mantle_melting.cc](https://github.com/geodynamics/aspect/blob/main/source/material_model/reaction_model/katz2003_mantle_melting.cc). + +This is defined in the parameter file as: +```{literalinclude} reactive_fluid_transport_model.part.prm +``` diff --git a/cookbooks/glacial_cycling_melt_transport/doc/reactive_fluid_transport_model.part.prm b/cookbooks/glacial_cycling_melt_transport/doc/reactive_fluid_transport_model.part.prm new file mode 100644 index 00000000000..32669717051 --- /dev/null +++ b/cookbooks/glacial_cycling_melt_transport/doc/reactive_fluid_transport_model.part.prm @@ -0,0 +1,16 @@ +subsection Material model + set Model name = reactive fluid transport + + subsection Reactive Fluid Transport Model + set Base model = visco plastic + set Fluid-solid reaction scheme = katz2003 + + # Katz 2003 melting model + subsection Katz 2003 model + set Reference permeability = 1e-7 + set Melt extraction depth = 30.0e3 + set Freezing rate = 0.5 + set Melting time scale for operator splitting = 10e2 + set Exponential melt weakening factor = 20 + end + end \ No newline at end of file diff --git a/cookbooks/glacial_cycling_melt_transport/doc/traction_boundary.part.prm b/cookbooks/glacial_cycling_melt_transport/doc/traction_boundary.part.prm new file mode 100644 index 00000000000..ee97d1d977e --- /dev/null +++ b/cookbooks/glacial_cycling_melt_transport/doc/traction_boundary.part.prm @@ -0,0 +1,13 @@ +# Prescribe a fixed vertical traction on the top boundary +# Adding parabolic ice load on top boundary at the middle of the model. The ice load is linearly +# decreasing over time. +subsection Boundary traction model + set Prescribed traction boundary indicators = top: function + + subsection Function + set Variable names = x,y,t + set Function constants = rho=917, g=9.81, thickness=1e3, x_center=200e3, half_width=100e3, t1=1e6, t2=1.02e6 + set Function expression = 0; - ( if(t<=t1, rho*g*thickness, if(t20.e3 && x<=40.e3, 1, 0); \ + if( x>40.e3 && x<=80.e3, 1, 0); \ + if( x>80.e3, 1, 0) + end +end + +subsection Discretization + set Composition polynomial degree = 2 + set Stokes velocity polynomial degree = 2 + set Temperature polynomial degree = 1 + set Use discontinuous composition discretization = false +end + +# Model geometry (200x200 km, 20 km spacing) +subsection Geometry model + set Model name = box + subsection Box + set X repetitions = 20 + set Y repetitions = 10 + set X extent = 400e3 + set Y extent = 200e3 + end +end + +# Mesh refinement specifications (refined to 1.25 km adaptively) +subsection Mesh refinement + set Initial adaptive refinement = 0 + set Initial global refinement = 4 + set Time steps between mesh refinement = 0 +end + +# Mesh deformation and free surface +# 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 + + subsection Free surface + set Surface velocity projection = normal + end +end + +subsection Boundary velocity model + set Tangential velocity boundary indicators = bottom + set Zero velocity boundary indicators = left, right +end + +# Prescribe a fixed vertical traction on the top boundary +# Adding parabolic ice load on top boundary at the middle of the model. The ice load is linearly +# decreasing over time. +subsection Boundary traction model + set Prescribed traction boundary indicators = top: function + + subsection Function + set Variable names = x,y,t + set Function constants = rho=917, g=9.81, thickness=1e3, x_center=200e3, half_width=100e3, t1=1e6, t2=1.02e6 + set Function expression = 0; - ( if(t<=t1, rho*g*thickness, if(t=180.e3, 1, 0); \ + if(y<180.e3 && y>=160.e3, 1, 0); \ + if(y<160.e3 && y>=120.e3, 1, 0); \ + if(y<120.e3 && y>=0.e3, 1, 0); + end +end + +# Composition boundary conditions +subsection Boundary composition model + set Fixed composition boundary indicators = bottom, left, right + set Allow fixed composition on outflow boundaries = true + 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 = 1753 + 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.065 mW/m^22 +# 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; +# Heat Production - asthenosphere (A4)= 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) = 713 K +# - mantle (ts3) = 1053 K +# - asthenosphere (ts4)= 1693 K +# Surface Heat Flow - lower crust (qs2) = 0.045 W/m^2; +# - mantle (qs3) = 0.04 W/m^2; +# - asthenosphere (qs4)= 0.04 W/m^2; +subsection Initial temperature model + set Model name = function + + subsection Function + set Variable names = x,y + set Function constants = h=200e3, ts1=273, ts2=713, ts3=1053, ts4=1693,\ + A1=1.e-6, A2=0.25e-6, A3=0.0, A4=0.0,\ + k1=2.5, k2=2.5, k3=2.5, k4=2.5, \ + qs1=0.065, qs2= 0.045, qs3=0.04, qs4=0.04 + 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), \ + if( (h-y)>40.e3 && (h-y)<=80.e3, \ + ts3 + (qs3/k3)*(h-y-40.e3) - (A3*(h-y-40.e3)*(h-y-40.e3))/(2.0*k3), \ + ts4 + ((h-y-80.e3)/1.e3)*0.5 ) ) ); + 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, adiabatic heating, shear heating + subsection Compositional heating + set Use compositional field for heat production averaging = 1, 0, 0, 1, 1, 1, 1 + set Compositional heating values = 0., 0., 0., 1.e-6, 0.25e-6, 0.0, 0.0 + end + subsection Adiabatic heating + set Use simplified adiabatic heating = true + end +end + +subsection Material model + set Model name = reactive fluid transport + + subsection Reactive Fluid Transport Model + set Base model = visco plastic + set Fluid-solid reaction scheme = katz2003 + + # Katz 2003 melting model + subsection Katz 2003 model + set Reference permeability = 1e-7 + set Melt extraction depth = 30.0e3 + set Freezing rate = 0.5 + set Melting time scale for operator splitting = 10e2 + set Exponential melt weakening factor = 20 + end + end + + 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 = 1e20 + set Maximum viscosity = 1e24 + + # 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. + + set Densities = background: 3300, \ + crust_upper: 2700, \ + crust_lower: 2900, \ + mantle_lithosphere: 3300, \ + asthenosphere: 3300 + + set Thermal expansivities = 2e-5 + + # Harmonic viscosity averaging + set Viscosity averaging scheme = geometric + set Viscous flow law = composite + + # 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 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 = background: 7.37e-15, \ + crust_upper: 1.37e-26, \ + crust_lower: 5.71e-23, \ + mantle_lithosphere: 7.37e-15, \ + asthenosphere: 7.37e-15 + + set Stress exponents for dislocation creep = background: 3.5, \ + crust_upper: 4.0, \ + crust_lower: 3.0, \ + mantle_lithosphere: 3.5, \ + asthenosphere: 3.5 + + set Activation energies for dislocation creep = background: 530.e3, \ + crust_upper: 223.e3, \ + crust_lower: 345.e3, \ + mantle_lithosphere: 530.e3, \ + asthenosphere: 530.e3 + + set Activation volumes for dislocation creep = background: 18.e-6, \ + crust_upper: 0, \ + crust_lower: 0, \ + mantle_lithosphere: 18.e-6, \ + asthenosphere: 18.e-6 + + set Prefactors for diffusion creep = background: 2.37e-15, \ + crust_upper: 1.00e-50, \ + crust_lower: 1.00e-50, \ + mantle_lithosphere: 1.00e-50, \ + asthenosphere: 1.00e-50 + + set Grain size exponents for diffusion creep = background: 3.0, \ + crust_upper: 0., \ + crust_lower: 0., \ + mantle_lithosphere: 0., \ + asthenosphere: 0. + + set Activation energies for diffusion creep = background: 375.e3, \ + crust_upper: 0., \ + crust_lower: 0., \ + mantle_lithosphere: 0., \ + asthenosphere: 0. + + set Activation volumes for diffusion creep = background: 2.e-6, \ + crust_upper: 0, \ + crust_lower: 0, \ + mantle_lithosphere: 0, \ + asthenosphere: 0 + # Plasticity parameters + set Cohesions = 1.e50 + + end +end + +# Gravity model +subsection Gravity model + set Model name = vertical + + subsection Vertical + set Magnitude = 9.815 + 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, 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 + +# postprocessing +subsection Postprocess + set List of postprocessors = visualization, composition statistics, melt statistics, velocity statistics, basic statistics, heat flux statistics, heat flux densities, pressure statistics, temperature statistics, particles + + # We mainly want to look at material properties of the solid and the melt. + subsection Visualization + set List of output variables = material properties, melt material properties, melt fraction, strain rate, principal stress, stress second invariant, strain rate tensor, shear stress + set Time between graphical output = 0 + set Output format = vtu + set Time between graphical output = 100.e3 + set Interpolate output = true + set Write higher order output = true + set Output format = vtu + + subsection Material properties + set List of material properties = density, viscosity + end + + subsection Melt material properties + set List of properties = compaction viscosity, permeability, fluid density, is melt cell + end + end + + subsection Particles + set Time between data output = 500.e3 + set Data output format = vtu + end +end + +# Checkpointing +subsection Checkpointing + set Steps between checkpoint = 10 +end + +# Termination criteria +subsection Termination criteria + set Termination criteria = end time +end diff --git a/doc/sphinx/user/cookbooks/geophysical-setups.md b/doc/sphinx/user/cookbooks/geophysical-setups.md index e0594a059fd..a501a3cd48c 100644 --- a/doc/sphinx/user/cookbooks/geophysical-setups.md +++ b/doc/sphinx/user/cookbooks/geophysical-setups.md @@ -106,6 +106,7 @@ cookbooks/continental_extension/doc/continental_extension.md cookbooks/inner_core_convection/doc/inner_core_convection.md cookbooks/lower_crustal_flow/doc/lower_crustal_flow.md cookbooks/global_melt/doc/global_melt.md +cookbooks/glacial_cycling_melt_transport/doc/glacial_cycling_melt_transport.md cookbooks/mid_ocean_ridge/doc/mid_ocean_ridge.md cookbooks/grain_size_ridge/doc/grain_size_ridge.md cookbooks/transform_fault_behn_2007/doc/transform_fault_behn_2007.md