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.

Sign up for GitHub

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] Add 2D Dynamic Subduction Cookbook with Two-Phase Reactive Fluid Transport #5785

Open
wants to merge 4 commits into
base: main
Choose a base branch
from
+867 −0
Open

[WIP] Add 2D Dynamic Subduction Cookbook with Two-Phase Reactive Fluid Transport #5785

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
d6fb760
Add cookbook for water cyclying in a fixed 2D slab
danieldouglas92 Jun 4, 2024
deeea8f
fix stuff
danieldouglas92 Jun 4, 2024
1f22ea9
Add dynamic 2D two-phase fluid subduction cookbook
danieldouglas92 Jun 5, 2024
187066f
Merge branch 'geodynamics:main' into add_fixed_slab_subduction_cookbook
danieldouglas92 Sep 17, 2024
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
80 changes: 80 additions & 0 deletions cookbooks/dynamic_subduction_with_two_phases/2D_tian_slab.wb
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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}]
}
]
}
12 changes: 12 additions & 0 deletions ...ks/dynamic_subduction_with_two_phases/doc/dynamic_subduction_with_two_phases.md
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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.

Check warning on line 6 in cookbooks/dynamic_subduction_with_two_phases/doc/dynamic_subduction_with_two_phases.md

View workflow job for this annotation

GitHub Actions / Check for new typos 

perhaps "boundaires" should be "boundaries".

```{figure-md} fig:initial-temperature-hydration
<img src="fixed_slab_temperature_fluid.png" style="width:90.0%" />

Schematic diagram showing the model region of a subduction zone, and the actual schematic diagram of the model.
```
Binary file added cookbooks/dynamic_subduction_with_two_phases/doc/fixed_slab_temperature_fluid.png
View file
Open in desktop
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
252 changes: 252 additions & 0 deletions cookbooks/dynamic_subduction_with_two_phases/uncoupled-fixed-slab.prm
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
Loading
Loading