diff --git a/cookbooks/dynamic_subduction_with_two_phases/2D_tian_slab.wb b/cookbooks/dynamic_subduction_with_two_phases/2D_tian_slab.wb new file mode 100644 index 00000000000..44c8c08ede5 --- /dev/null +++ b/cookbooks/dynamic_subduction_with_two_phases/2D_tian_slab.wb @@ -0,0 +1,80 @@ +// Composition 1 is bound fluid +// Composition 2 is peridotite +// Composition 3 is gabbro +// Composition 4 is MORB +// Composition 5 is sediment +// Composition 6 is weak zone +{ + "version": "0.6", + "gravity model":{"model":"uniform", "magnitude":10}, + "cross section":[[0,0],[50e3,0]], + "surface temperature":273, "potential mantle temperature":1573, + "thermal expansion coefficient":3.1e-5, "specific heat":1000, "thermal diffusivity":1.0e-6, + "features": + [ + {"model": "mantle layer", "name": "peridotite mantle", "coordinates": [[-500e3, 100e3], [-500e3, -100e3], [2000e3, -100e3], [2000e3, 100e3]], + "min depth": 0, "max depth":10000e3, + "composition models": + [{"model": "uniform", "min depth":0.0, "max depth":10000e3, "compositions": [2]}]}, + + {"model": "oceanic plate", "name": "Overriding Plate", + "coordinates": [[300e3, -100e3], [300e3, 100e3], [2000e3, 100e3], [2000e3, -100e3]], + "min depth": 0, + "max depth": 100e3, + "composition models": [{"model": "uniform", "compositions": [2], "min depth":0, "max depth": 100e3}], + "temperature models": [{"model": "half space model", "ridge coordinates": [[[1500e3,-100e3],[1500e3,100e3]]], + "spreading velocity": 0.05, "min depth": 0, "max depth":100e3}] + }, + + {"model": "oceanic plate", "name": "Subducting Plate", + "coordinates": [[300e3, -100e3], [300e3, 100e3], [-500e3, 100e3], [-500e3, -100e3]], + "min depth": 0, + "max depth": 150e3, + "composition models": + [ + {"model": "uniform", "compositions": [1], "min depth":0.0, "max depth": 10e3, "fractions": [0.02]}, + {"model": "uniform", "compositions": [1], "min depth":10e3, "max depth": 20e3, "fractions": [0.01]}, + {"model": "uniform", "compositions": [1], "min depth":20e3, "max depth": 30e3, "fractions": [0.005]}, + {"model": "uniform", "compositions": [5], "min depth":0.0, "max depth": 10e3, "operation":"add"}, + {"model": "uniform", "compositions": [4], "min depth":10e3, "max depth": 20e3, "operation":"add"}, + {"model": "uniform", "compositions": [3], "min depth":20e3, "max depth": 30e3, "operation":"add"}], + + "temperature models": [{"model": "half space model", "ridge coordinates": [[[-500e3,-100e3],[-500e3,100e3]]], + "spreading velocity": 0.05, "min depth": 0, "max depth": 150e3, "bottom temperature": -1}] + }, + {"model":"subducting plate", "name":"Interface", + + "coordinates":[[310e3,-100e3],[310e3,100e3]], + "dip point":[1e7,0],"max depth":50e3, + "segments":[{"length":150e3,"thickness":[150e3],"top truncation":[-50e3],"angle":[0,45]}], + + "composition models":[ + {"model":"uniform", "compositions":[6], "min distance slab top":0, "max distance slab top":20e3}] + }, + + {"model":"subducting plate", "name":"Slab", + + "coordinates":[[300e3,-100e3],[300e3,100e3]], + "dip point":[1e7,0],"max depth":1000e3, + "segments":[{"length":150e3,"thickness":[150e3],"top truncation":[-50e3],"angle":[0,45]}], + + "composition models":[ + {"model":"uniform", "compositions":[1], "min distance slab top":0, "max distance slab top":10e3, "fractions": [0.02]}, + {"model":"uniform", "compositions":[1], "min distance slab top":10e3, "max distance slab top":20e3, "fractions": [0.01]}, + {"model":"uniform", "compositions":[1], "min distance slab top":20e3, "max distance slab top":30e3, "fractions": [0.005]}, + {"model":"uniform", "compositions":[5], "min distance slab top":0, "max distance slab top":10e3, "operation":"add"}, + {"model":"uniform", "compositions":[4], "min distance slab top":10e3, "max distance slab top":20e3, "operation":"add"}, + {"model":"uniform", "compositions":[3], "min distance slab top":20e3, "max distance slab top":30e3, "operation":"add"}], + + "temperature models":[{"model":"mass conserving", + "reference model name": "half space model", + "density":3300, "thermal conductivity":3.3, "adiabatic heating":true, + "plate velocity":0.05, + "ridge coordinates":[[[-500e3,-100e3],[-500e3,100e3]]], + "coupling depth":80e3, + "forearc cooling factor":5.0, + "taper distance":0, + "min distance slab top":-200e3, "max distance slab top":300e3}] + } + ] +} diff --git a/cookbooks/dynamic_subduction_with_two_phases/doc/dynamic_subduction_with_two_phases.md b/cookbooks/dynamic_subduction_with_two_phases/doc/dynamic_subduction_with_two_phases.md new file mode 100644 index 00000000000..f7d55495609 --- /dev/null +++ b/cookbooks/dynamic_subduction_with_two_phases/doc/dynamic_subduction_with_two_phases.md @@ -0,0 +1,12 @@ +(sec:cookbooks:dynamic_subduction_with_two_phases)= +# Dynamic 2D Subduction Model with Reactive Two Phase Fluid Flow + +*This section was contributed by Daniel Douglas.* + +In this cookbook we extend the simplified kinematic subduction model built in INSERT_LINK_HERE to a more dynamic model of subduction in a 2D Cartesian box. The model begins with a slab that has just started to subduct with a leading slab dip of 45 degrees, and a weak plate interface to facilitate subduction. The initial conditions for this model were constructed using the Geodynamic World Builder, see this [tutorial](https://gwb.readthedocs.io/en/latest/user_manual/cookbooks/simple_subduction_2d_cartesian/doc/README.html) for an in depth explanation on how to setup your own World Builder file, and the file used to generate this model can be found [here](https://www.github.com/geodynamics/aspect/blob/main/cookbooks/dynamic_subduction_with_two_phases/WORLDBUILDER_FILE.WB). While we still kinematically impose a convergence velocity of 5cm/yr, this is only prescribed on the upper left boundary. We use the 'box with lithosphere boundary indicators' geometry model to only impose this velocity on the shallowest 100 km of the left boundary, and for depths larger than 100 km the boundary in left open to allow inflow/outflow. On the right boundary, we impose a free slip boundary condition on the shallowest 100 km, and an open boundary for depths greater than 100 km. For the top and bottom boundaries, a free slip boundary condition is imposed. With the exception of the entire left boundary, the composition is not held fixed. We fix the temperature on the boundaries for the top, bottom, and the shallowest 100 km on both the left and right boundaires, and allow it to dynamically evolve for depths greater than 100 km on the left and right boundaries. The subducting plate (slab and unsubducted oceanic lithosphere) has a layered composition comprised of (from shallowest to deepest): 10 km of sediment containing 2 wt% water, 10 km of mid-ocean ridge basalt (MORB) hydrated to 1 wt% water, and 10 km of gabbro hydrated to 0.5 wt% water. The rest of the model is composed of peridotite, and there is no free fluid at the beginning of the model. + +```{figure-md} fig:initial-temperature-hydration + + +Schematic diagram showing the model region of a subduction zone, and the actual schematic diagram of the model. +``` diff --git a/cookbooks/dynamic_subduction_with_two_phases/doc/fixed_slab_temperature_fluid.png b/cookbooks/dynamic_subduction_with_two_phases/doc/fixed_slab_temperature_fluid.png new file mode 100644 index 00000000000..670bc29bf0e Binary files /dev/null and b/cookbooks/dynamic_subduction_with_two_phases/doc/fixed_slab_temperature_fluid.png differ diff --git a/cookbooks/dynamic_subduction_with_two_phases/uncoupled-fixed-slab.prm b/cookbooks/dynamic_subduction_with_two_phases/uncoupled-fixed-slab.prm new file mode 100644 index 00000000000..173ef6d45c9 --- /dev/null +++ b/cookbooks/dynamic_subduction_with_two_phases/uncoupled-fixed-slab.prm @@ -0,0 +1,252 @@ +# This is a cookbook which simulates a simple 2D subduction system with +# a hydrated sediment layer, MORB layer, gabbro layer, and peridotite +# layer. The layers dehydrate as they are subducted, with the solid-fluid +# reactions governed by an approximation published by Tian et al., 2018. +set World builder file = 2D_tian_slab.wb +set Adiabatic surface temperature = 1600 +set Nonlinear solver scheme = iterated Advection and Stokes +set Output directory = tian-slab-uncoupled-no-smoothin +set Max nonlinear iterations = 50 +set Nonlinear solver tolerance = 1e-5 +set Dimension = 2 +set End time = 1e6 +set Surface pressure = 0 +# If this is larger than the melt generation then that means that the fluid will instantly dehydrate +# Perhaps this is unstable, especially if the solver tolerance is 1? Make this 10x less than the +# reaction time step in the volatile material model. +set Maximum time step = 1e3 +set Use operator splitting = true +set Resume computation = auto + +subsection Solver parameters + subsection Stokes solver parameters + set GMRES solver restart length = 150 + set Number of cheap Stokes solver steps = 2000 + set Use full A block as preconditioner = true + set Linear solver tolerance = 1e-8 + set Maximum number of expensive Stokes solver steps = 100 + end + subsection Operator splitting parameters + set Reaction time step = 1e2 + set Reaction time steps per advection step = 10 + end +end + +subsection Discretization + subsection Stabilization parameters + set beta = 0.052 # 0.5 + set cR = 0.11 # 1.0 + end +end + +subsection Checkpointing + set Steps between checkpoint = 100 +end + +subsection Compositional fields + set Number of fields = 7 + set Names of fields = porosity, bound_fluid, peridotite, gabbro, MORB, sediment, weak_zone1 + set Compositional field methods = darcy field, field, field, field, field, field, field +end + +subsection Temperature field + set Temperature method = field +end + +subsection Initial composition model + set List of model names = world builder + subsection World builder + set List of relevant compositions = porosity, bound_fluid, peridotite, gabbro, MORB, sediment, weak_zone1 + end +end + +subsection Boundary composition model + set Fixed composition boundary indicators = left, left lithosphere + set List of model names = initial composition +end + +subsection Geometry model + set Model name = box with lithosphere boundary indicators + subsection Box with lithosphere boundary indicators + set Lithospheric thickness = 100e3 + set X extent = 800e3 + set Y extent = 400e3 + + set X repetitions = 8 + set Y repetitions = 3 + set Y repetitions lithosphere = 1 + + set Box origin X coordinate = 100e3 + set Box origin Y coordinate = 0 + end +end + + +subsection Gravity model + set Model name = function + subsection Function + set Variable names = x,y + set Function expression = 0; -9.81 + end +end + + +subsection Initial temperature model + set Model name = world builder +end + +subsection Boundary temperature model + set Fixed temperature boundary indicators = top, bottom, left lithosphere, right lithosphere + set List of model names = initial temperature +end + + +subsection Material model + + set Model name = reactive fluid transport + set Material averaging = harmonic average only viscosity + + subsection Reactive Fluid Transport Model + set Base model = visco plastic + set Reference fluid density = 1000 + set Shear to bulk viscosity ratio = 0.1 + set Reference fluid viscosity = 1 + set Reference permeability = 1e-6 + set Exponential fluid weakening factor = 30 + set Fluid compressibility = 0 + set Fluid reaction time scale for operator splitting = 5e4 + set Fluid-solid reaction scheme = tian approximation + set Maximum weight percent water in peridotite = 2 + set Maximum weight percent water in gabbro = 1 + set Maximum weight percent water in MORB = 2 + set Maximum weight percent water in sediment = 3 + end + + subsection Visco Plastic + set Viscosity averaging scheme = harmonic + set Viscous flow law = composite + set Prefactors for diffusion creep = background:4.5e-15, \ + porosity:4.5e-15, \ + bound_fluid:4.5e-15, \ + peridotite:4.5e-15, \ + gabbro:4.5e-15, \ + MORB:4.5e-15, \ + sediment:4.5e-15, \ + weak_zone1:4.5e-15 + + set Stress exponents for diffusion creep = 1.0 + + set Activation volumes for diffusion creep = background:8.2e-6, \ + porosity:8.2e-6, \ + bound_fluid:8.2e-6, \ + peridotite:8.2e-6, \ + gabbro:8.2e-6, \ + MORB:8.2e-6, \ + sediment:8.2e-6, \ + weak_zone1:8.2e-6 + + set Activation energies for diffusion creep = background:375e3, \ + porosity:375e3, \ + bound_fluid:375e3, \ + peridotite:375e3, \ + gabbro:375e3, \ + MORB:375e3, \ + sediment:375e3, \ + weak_zone1:375e3 + + set Prefactors for dislocation creep = background:7.4e-15, \ + porosity:7.4e-15, \ + bound_fluid:7.4e-15, \ + peridotite:7.4e-15, \ + gabbro:7.4e-15, \ + MORB:7.4e-15, \ + sediment:7.4e-15, \ + weak_zone1:1e-50 + + set Stress exponents for dislocation creep = 3.5 + + set Activation volumes for dislocation creep = background:14e-6, \ + porosity:14e-6, \ + bound_fluid:14e-6, \ + peridotite:14e-6, \ + gabbro:14e-6, \ + MORB:14e-6, \ + sediment:14e-6, \ + weak_zone1:14e-6 + + set Activation energies for dislocation creep = background:530e3, \ + porosity:530e3, \ + bound_fluid:530e3, \ + peridotite:530e3, \ + gabbro:530e3, \ + MORB:530e3, \ + sediment:530e3, \ + weak_zone1:530e3 + + set Angles of internal friction = 0 + + set Cohesions = background:1e50, \ + porosity:1e50, \ + bound_fluid:1e50, \ + peridotite:1e50, \ + gabbro:1e50, \ + MORB:1e50, \ + sediment:1e50, \ + weak_zone1:1e6 + + set Minimum viscosity = 5e19 + set Maximum viscosity = 5e23 + end +end + +subsection Mesh refinement + set Coarsening fraction = 0.05 + set Refinement fraction = 0.8 + set Initial adaptive refinement = 2 + set Initial global refinement = 3 + set Strategy = isosurfaces, composition threshold, minimum refinement function + set Time steps between mesh refinement = 2 + subsection Isosurfaces + set Isosurfaces = max, max, bound_fluid: 0.005|0.04 + end + + # minimum of 4 global refinements + subsection Minimum refinement function + set Function expression = 3 + end + + # refine where the porosity is bigger than 1e-6 + subsection Composition threshold + set Compositional field thresholds = 1e-6, 10.0, 10.0, 10.0, 10.0, 10.0, 10.0 + end +end + +subsection Boundary velocity model + set Tangential velocity boundary indicators = top, bottom, right lithosphere + set Prescribed velocity boundary indicators = left lithosphere: function + subsection Function + set Function expression = 0.05; 0.0 + end +end + +subsection Boundary traction model + set Prescribed traction boundary indicators = left:initial lithostatic pressure, right:initial lithostatic pressure + subsection Initial lithostatic pressure + set Representative point = 100e3, 0 + end +end + +subsection Melt settings + set Include melt transport = false +end + +subsection Postprocess + set List of postprocessors = visualization, composition statistics, velocity statistics + subsection Visualization + set List of output variables = density, viscosity, thermal expansivity, named additional outputs + set Output format = vtu + set Time between graphical output = 1e4 + set Interpolate output = true + set Number of grouped files = 0 + end +end diff --git a/cookbooks/fixed_tian_slab/coupled-fixed-slab.prm b/cookbooks/fixed_tian_slab/coupled-fixed-slab.prm new file mode 100644 index 00000000000..1602f2420c3 --- /dev/null +++ b/cookbooks/fixed_tian_slab/coupled-fixed-slab.prm @@ -0,0 +1,262 @@ +# This is a cookbook which simulates a simple 2D subduction system with +# a hydrated sediment layer, MORB layer, gabbro layer, and peridotite +# layer. The layers dehydrate as they are subducted, with the solid-fluid +# reactions governed by an approximation published by Tian et al., 2018. +set World builder file = /home/900358141/2024_aspect_hack/fixed_slab_cookbook/2D_tian_slab.wb +set Adiabatic surface temperature = 1600 +set Nonlinear solver scheme = iterated Advection and Stokes +set Output directory = tian_fixed_slab +set Max nonlinear iterations = 50 +set Nonlinear solver tolerance = 1e-5 +set Dimension = 2 +set End time = 1e6 +set Surface pressure = 0 +# If this is larger than the melt generation then that means that the fluid will instantly dehydrate +# Perhaps this is unstable, especially if the solver tolerance is 1? Make this 10x less than the +# reaction time step in the volatile material model. +set Maximum time step = 1e3 +set Use operator splitting = true +set Resume computation = auto + +subsection Solver parameters + subsection Stokes solver parameters + # set Stokes solver type = block GMG + set GMRES solver restart length = 150 + set Number of cheap Stokes solver steps = 2000 + set Use full A block as preconditioner = true + set Linear solver tolerance = 1e-8 + set Maximum number of expensive Stokes solver steps = 100 + end + subsection Operator splitting parameters + set Reaction time step = 1e2 + set Reaction time steps per advection step = 10 + end +end + +subsection Discretization + # We choose relatively large values for the stabilization parameters: + # # However, note that in an application model with a higher resolution, + # # we would choose much smaller values for the stabilization parameters. + subsection Stabilization parameters + set beta = 0.5 + set cR = 1 + end +end + +subsection Checkpointing + set Steps between checkpoint = 100 + # set Steps between checkpoint = 0 +end + +subsection Compositional fields + set Number of fields = 7 + set Names of fields = porosity, bound_fluid, peridotite, gabbro, MORB, sediment, weak_zone + set Compositional field methods = field, field, static, static, static, static, static +end + +subsection Temperature field + # set Temperature method = field + set Temperature method = static +end + +subsection Initial composition model + set List of model names = world builder + subsection World builder + set List of relevant compositions = porosity, bound_fluid, peridotite, gabbro, MORB, sediment, weak_zone + end +end + +subsection Boundary composition model + set Fixed composition boundary indicators = left, top + set List of model names = initial composition +end + +subsection Geometry model + set Model name = box + subsection Box + set X extent = 1000e3 + set Y extent = 500e3 + set Y repetitions = 1 + set X repetitions = 2 + end +end + + +subsection Gravity model + set Model name = function + subsection Function + set Variable names = x,y + set Function expression = 0; -9.81 + end +end + + +subsection Initial temperature model + set Model name = world builder +end + +subsection Boundary temperature model + set Fixed temperature boundary indicators = top, bottom, left, right + set List of model names = initial temperature +end + + +subsection Material model + + set Model name = reactive fluid transport + # set Model name = visco plastic + set Material averaging = harmonic average only viscosity + + subsection Reactive Fluid Transport Model + set Base model = visco plastic + set Reference fluid density = 1000 + set Shear to bulk viscosity ratio = 0.1 + set Reference fluid viscosity = 1 + set Reference permeability = 1e-6 + set Exponential fluid weakening factor = 27 + set Fluid compressibility = 0 + set Fluid reaction time scale for operator splitting = 5e4 + set Fluid-solid reaction scheme = tian approximation + set Maximum weight percent water in peridotite = 2 + set Maximum weight percent water in gabbro = 1 + set Maximum weight percent water in MORB = 2 + set Maximum weight percent water in sediment = 3 + end + + subsection Visco Plastic + set Viscosity averaging scheme = harmonic + set Viscous flow law = composite + set Prefactors for diffusion creep = background:4.5e-15, \ + porosity:4.5e-15, \ + bound_fluid:4.5e-15, \ + peridotite:4.5e-15, \ + gabbro:4.5e-15, \ + MORB:4.5e-15, \ + sediment:4.5e-15, \ + weak_zone:5e-19 + + set Stress exponents for diffusion creep = 1.0 + + set Activation volumes for diffusion creep = background:8.2e-6, \ + porosity:8.2e-6, \ + bound_fluid:8.2e-6, \ + peridotite:8.2e-6, \ + gabbro:8.2e-6, \ + MORB:8.2e-6, \ + sediment:8.2e-6, \ + weak_zone:8.2e-6 + + set Activation energies for diffusion creep = background:375e3, \ + porosity:375e3, \ + bound_fluid:375e3, \ + peridotite:375e3, \ + gabbro:375e3, \ + MORB:375e3, \ + sediment:375e3, \ + weak_zone:375e3 + + set Prefactors for dislocation creep = background:7.4e-15, \ + porosity:7.4e-15, \ + bound_fluid:7.4e-15, \ + peridotite:7.4e-15, \ + gabbro:7.4e-15, \ + MORB:7.4e-15, \ + sediment:7.4e-15, \ + weak_zone:1e-50 + + set Stress exponents for dislocation creep = 3.5 + + set Activation volumes for dislocation creep = background:14e-6, \ + porosity:14e-6, \ + bound_fluid:14e-6, \ + peridotite:14e-6, \ + gabbro:14e-6, \ + MORB:14e-6, \ + sediment:14e-6, \ + weak_zone:14e-6 + + set Activation energies for dislocation creep = background:530e3, \ + porosity:530e3, \ + bound_fluid:530e3, \ + peridotite:530e3, \ + gabbro:530e3, \ + MORB:530e3, \ + sediment:530e3, \ + weak_zone:530e3 + + set Angles of internal friction = 0 + + set Cohesions = background:1e50, \ + porosity:1e50, \ + bound_fluid:1e50, \ + peridotite:1e50, \ + gabbro:1e50, \ + MORB:1e50, \ + sediment:1e50, \ + weak_zone:1e6 + + set Minimum viscosity = 5e19 + set Maximum viscosity = 5e23 + # set Use plastic damper = true + # set Plastic damper viscosity = 1e21 + + # set Strain weakening mechanism = plastic weakening with plastic strain only + # set Friction strain weakening factors = 0.5 + # set Cohesion strain weakening factors = 0.5 + # set Start plasticity strain weakening intervals = 0 + # set End plasticity strain weakening intervals = 1 + # set Thermal expansivities = 0 + end +end + +subsection Mesh refinement + set Coarsening fraction = 0.05 + set Refinement fraction = 0.8 + set Initial adaptive refinement = 2 + set Initial global refinement = 3 + set Strategy = isosurfaces, composition threshold, minimum refinement function + set Time steps between mesh refinement = 2 + subsection Isosurfaces + set Isosurfaces = max, max, bound_fluid: 0.005|0.04 + end + + # minimum of 5 global refinements + subsection Minimum refinement function + set Function expression = 3 + end + + # refine where the porosity is bigger than 1e-6 + subsection Composition threshold + set Compositional field thresholds = 1e-6, 10.0, 10.0, 10.0, 10.0, 10.0, 10.0 + end +end + +subsection Boundary velocity model + set Tangential velocity boundary indicators = top, bottom, right + # set Zero velocity boundary indicators = bottom +end + +subsection Boundary traction model + set Prescribed traction boundary indicators = left: initial lithostatic pressure + subsection Initial lithostatic pressure + set Representative point = 0, 0 + end +end + +subsection Melt settings + set Include melt transport = true +end + +subsection Postprocess + set List of postprocessors = visualization, composition statistics, velocity statistics + subsection Visualization + set List of output variables = density, viscosity, thermal expansivity, named additional outputs, melt material properties, melt fraction + subsection Melt material properties + set List of properties = fluid density, permeability, fluid viscosity, compaction viscosity + end + set Output format = vtu + set Time between graphical output = 10e3 + set Interpolate output = true + set Number of grouped files = 0 + end +end diff --git a/cookbooks/fixed_tian_slab/uncoupled-fixed-slab.prm b/cookbooks/fixed_tian_slab/uncoupled-fixed-slab.prm new file mode 100644 index 00000000000..748393f6339 --- /dev/null +++ b/cookbooks/fixed_tian_slab/uncoupled-fixed-slab.prm @@ -0,0 +1,260 @@ +# This is a cookbook which simulates a simple 2D subduction system with +# a hydrated sediment layer, MORB layer, gabbro layer, and peridotite +# layer. The layers dehydrate as they are subducted, with the solid-fluid +# reactions governed by an approximation published by Tian et al., 2018. +set World builder file = /home/900358141/2024_aspect_hack/fixed_slab_cookbook/2D_tian_slab.wb +set Adiabatic surface temperature = 1600 +set Nonlinear solver scheme = iterated Advection and Stokes +set Output directory = tian-slab-uncoupled +set Max nonlinear iterations = 50 +set Nonlinear solver tolerance = 1e-5 +set Dimension = 2 +set End time = 1e6 +set Surface pressure = 0 +# If this is larger than the melt generation then that means that the fluid will instantly dehydrate +# Perhaps this is unstable, especially if the solver tolerance is 1? Make this 10x less than the +# reaction time step in the volatile material model. +set Maximum time step = 1e3 +set Use operator splitting = true +set Resume computation = auto + +subsection Solver parameters + subsection Stokes solver parameters + # set Stokes solver type = block GMG + set GMRES solver restart length = 150 + set Number of cheap Stokes solver steps = 2000 + set Use full A block as preconditioner = true + set Linear solver tolerance = 1e-8 + set Maximum number of expensive Stokes solver steps = 100 + end + subsection Operator splitting parameters + set Reaction time step = 1e2 + set Reaction time steps per advection step = 10 + end +end + +subsection Discretization + # We choose relatively large values for the stabilization parameters: + # # However, note that in an application model with a higher resolution, + # # we would choose much smaller values for the stabilization parameters. + subsection Stabilization parameters + set beta = 0.5 + set cR = 1 + end +end + +subsection Checkpointing + set Steps between checkpoint = 100 + # set Steps between checkpoint = 0 +end + +subsection Compositional fields + set Number of fields = 8 + set Names of fields = porosity, bound_fluid, peridotite, gabbro, MORB, sediment, weak_zone1, weak_zone2 + set Compositional field methods = darcy field, static, static, static, static, static, static, static + # set Compositional field methods = darcy field, field, field, field, field, field, field +end + +subsection Temperature field + set Temperature method = static +end + +subsection Initial composition model + set List of model names = world builder, function + subsection World builder + set List of relevant compositions = porosity, bound_fluid, peridotite, gabbro, MORB, sediment, weak_zone1 + end + subsection Function + set Coordinate system = cartesian + set Variable names = x, y + set Function expression = 0;0;0;0;0;0;0;if(x<=50e3 & y>=450e3, 1, 0) +end + +subsection Boundary composition model + set Fixed composition boundary indicators = left, top + set List of model names = initial composition +end + +subsection Geometry model + set Model name = box + subsection Box + set X extent = 1000e3 + set Y extent = 500e3 + set Y repetitions = 10 + set X repetitions = 5 + end +end + + +subsection Gravity model + set Model name = function + subsection Function + set Variable names = x,y + set Function expression = 0; -9.81 + end +end + + +subsection Initial temperature model + set Model name = world builder +end + +subsection Boundary temperature model + set Fixed temperature boundary indicators = top, bottom, left, right + set List of model names = initial temperature +end + + +subsection Material model + + set Model name = reactive fluid transport + set Material averaging = harmonic average only viscosity + + subsection Reactive Fluid Transport Model + set Base model = visco plastic + set Reference fluid density = 1000 + set Shear to bulk viscosity ratio = 0.1 + set Reference fluid viscosity = 1 + set Reference permeability = 1e-6 + set Exponential fluid weakening factor = 27 + set Fluid compressibility = 0 + set Fluid reaction time scale for operator splitting = 5e4 + set Fluid-solid reaction scheme = tian approximation + set Maximum weight percent water in peridotite = 2 + set Maximum weight percent water in gabbro = 1 + set Maximum weight percent water in MORB = 2 + set Maximum weight percent water in sediment = 3 + end + + subsection Visco Plastic + set Viscosity averaging scheme = harmonic + set Viscous flow law = composite + set Prefactors for diffusion creep = background:4.5e-15, \ + porosity:4.5e-15, \ + bound_fluid:4.5e-15, \ + peridotite:4.5e-15, \ + gabbro:4.5e-15, \ + MORB:4.5e-15, \ + sediment:4.5e-15, \ + weak_zone1:4.5e-15, \ + weak_zone2:4.5e-15 + + set Stress exponents for diffusion creep = 1.0 + + set Activation volumes for diffusion creep = background:8.2e-6, \ + porosity:8.2e-6, \ + bound_fluid:8.2e-6, \ + peridotite:8.2e-6, \ + gabbro:8.2e-6, \ + MORB:8.2e-6, \ + sediment:8.2e-6, \ + weak_zone1:8.2e-6, \ + weak_zone2:8.2e-6 + + set Activation energies for diffusion creep = background:375e3, \ + porosity:375e3, \ + bound_fluid:375e3, \ + peridotite:375e3, \ + gabbro:375e3, \ + MORB:375e3, \ + sediment:375e3, \ + weak_zone1:375e3 + + set Prefactors for dislocation creep = background:7.4e-15, \ + porosity:7.4e-15, \ + bound_fluid:7.4e-15, \ + peridotite:7.4e-15, \ + gabbro:7.4e-15, \ + MORB:7.4e-15, \ + sediment:7.4e-15, \ + weak_zone1:1e-50, \ + weak_zone2:1e-50 + + set Stress exponents for dislocation creep = 3.5 + + set Activation volumes for dislocation creep = background:14e-6, \ + porosity:14e-6, \ + bound_fluid:14e-6, \ + peridotite:14e-6, \ + gabbro:14e-6, \ + MORB:14e-6, \ + sediment:14e-6, \ + weak_zone1:14e-6, \ + weak_zone2:14e-6 + + set Activation energies for dislocation creep = background:530e3, \ + porosity:530e3, \ + bound_fluid:530e3, \ + peridotite:530e3, \ + gabbro:530e3, \ + MORB:530e3, \ + sediment:530e3, \ + weak_zone1:530e3, \ + weak_zone2:530e3 + + set Angles of internal friction = 0 + + set Cohesions = background:1e50, \ + porosity:1e50, \ + bound_fluid:1e50, \ + peridotite:1e50, \ + gabbro:1e50, \ + MORB:1e50, \ + sediment:1e50, \ + weak_zone1:1e6, \ + weak_zone2:1e6 + + set Minimum viscosity = 5e19 + set Maximum viscosity = 5e23 + # set Use plastic damper = true + # set Plastic damper viscosity = 1e21 + + # set Strain weakening mechanism = plastic weakening with plastic strain only + # set Friction strain weakening factors = 0.5 + # set Cohesion strain weakening factors = 0.5 + # set Start plasticity strain weakening intervals = 0 + # set End plasticity strain weakening intervals = 1 + # set Thermal expansivities = 0 + end +end + +subsection Mesh refinement + set Coarsening fraction = 0.05 + set Refinement fraction = 0.8 + set Initial adaptive refinement = 2 + set Initial global refinement = 4 + set Strategy = isosurfaces, composition threshold, minimum refinement function + set Time steps between mesh refinement = 2 + subsection Isosurfaces + set Isosurfaces = max, max, bound_fluid: 0.005|0.04 + end + + # minimum of 4 global refinements + subsection Minimum refinement function + set Function expression = 4 + end + + # refine where the porosity is bigger than 1e-6 + subsection Composition threshold + set Compositional field thresholds = 1e-6, 10.0, 10.0, 10.0, 10.0, 10.0, 10.0, 10.0 + end +end + +subsection Boundary velocity model + set Tangential velocity boundary indicators = top, bottom, right, left +end + +subsection Melt settings + set Include melt transport = false +end + +subsection Postprocess + set List of postprocessors = visualization, composition statistics, velocity statistics + subsection Visualization + set List of output variables = density, viscosity, thermal expansivity, named additional outputs + set Output format = vtu + set Time between graphical output = 0 + set Interpolate output = true + set Number of grouped files = 0 + end +end diff --git a/doc/sphinx/user/cookbooks/geophysical-setups.md b/doc/sphinx/user/cookbooks/geophysical-setups.md index 6146b3d9cfc..f62c596ffae 100644 --- a/doc/sphinx/user/cookbooks/geophysical-setups.md +++ b/doc/sphinx/user/cookbooks/geophysical-setups.md @@ -114,6 +114,7 @@ cookbooks/tian_parameterization_kinematic_slab/doc/tian_parameterization_kinemat cookbooks/mantle_convection_with_continents_in_annulus/doc/mantle_convection_in_annulus.md cookbooks/inclusions/doc/inclusions.md cookbooks/subduction_initiation/doc/subduction_initiation.md +cookbooks/dynamic_subduction_with_two_phases/doc/dynamic_subduction_with_two_phases.md cookbooks/vankeken_subduction/doc/vankeken_subduction.md cookbooks/2d_annulus_visualization/doc/2d_annulus_visualization.md cookbooks/tomography_based_plate_motions/doc/tomography_based_plate_motions.md