Advanced 2D-Subduction cookbook with reactive fluid transport

cookbooks/2D_subduction_with_two_phase_flow/doc/2D_subduction_two_phase_flow.md

@@ -0,0 +1,104 @@
+(sec:cookbooks:2D-subduction-two-phase-flow)=
+# 2D Subduction Model with Parameterized Solid-Fluid Reactions
+
+*This section was contributed by Daniel Douglas.*
+
+In this cookbook we expand on what was demonstrated in the {ref}`sec:cookbooks:tian_parameterization_kinematic_slab` cookbook by incrementally building towards a dynamic model of subduction that includes reactive fluid transport at varying degrees of fluid--solid coupling.
+
+This cookbook requires that ASPECT is compiled with the Geodynamic World Builder (GWB), which is enabled by setting `ASPECT_WITH_WORLD_BUILDER=ON`
+when configuring ASPECT with CMake (this is the default setting). The GWB is a powerful tool that allows ASPECT users to create complex initial
+conditions. In this example, we will use it to define the temperature and hydration state of a two-dimensional subduction zone. To use GWB with
+ASPECT, you must specify the path to a World Builder (.wb) file in the ASPECT input file, and indicate that the initial temperature and composition
+are generated using GWB. These settings look like this:
+```{literalinclude} input_world_builder.part.prm
+```
+
+```{figure-md} fig:model-overview
+<img src="model_overview.png" />
+
+ The model domain coloured by the model temperature. The subducting plate descends into the mantle at a constant dip of 45&deg;. The temperature of the mantle is fixed to 1573 K, the overriding plate is defined with a linear geotherm, and the subducting plate with a plate cooling model with a convergence rate of 3 cm/yr.
+```
+
+The World Builder file can be found at [cookbooks/2D_subduction_with_two_phase_flow](https://www.github.com/geodynamics/aspect/blob/main/cookbooks/
+2D_subduction_two_phase_flow/). However, this cookbook will only focus on the ASPECT side of the model. For more details on the World Builder and how
+to use it, specifically in the context of using geodynamic software like ASPECT for modeling subduction zones, please refer to the GWB manual. There
+are two comprehensive guides that are relevant to this cookbook, the first focuses on defining a [complex slab geometry and initial thermal
+distribution](https://gwb.readthedocs.io/en/latest/user_manual/cookbooks/simple_subduction_2d_cartesian/doc/README.html), and the second demonstrates
+how to define an [initial hydration state of a subducting plate](https://gwb.readthedocs.io/en/latest/user_manual/cookbooks/
+2d_cartesian_hydrated_slab/doc/README.html).
+
+The model domain is a rectangular box spanning 8700 km × 2900 km {numref}`fig:model-overview`. The trench is located at x = 4000 km. The subducting
+plate is 3000 km long and 120 km thick, while the overriding plate is 2500 km long and 80 km thick. The slab geometry is relatively simple: beginning
+at the trench, the slab bends to a dip of 45&deg; over a (slab) length of 300 km, then continues into the mantle at a constant 45&deg; dip for an
+additional 800 km (measured along the slab, not by depth). The subducting plate forms at a spreading center located 3000 km from the trench, and the
+temperature of the plate is initialized using a plate cooling model assuming a convergence rate of 3 cm/yr. This results in a 100 Myr old plate just
+before subduction. The overriding plate has a linear temperature gradient from a surface temperature of 273 K to 1573 K at the base of the plate. The
+subducting plate consists of multiple lithological layers: a 10 km thick sediment layer at the top, followed by a 10 km thick mid-ocean ridge basalt
+(MORB) layer, underlain by a 10 km thick gabbro layer. The remainder of the subducting plate, along with the mantle and the overriding plate, is
+comprised of peridotite. The initial hydration states of the layers within the subducting plate are as follows: sediment contains up to 2 wt% bound
+water, MORB up to 1 wt%, gabbro up to 0.5 wt%, and peridotite within the subducting plate up to 1 wt%. These values are then multiplied by a factor
+of 1.05, resulting in a 5% excess of bound water relative to the equilibrium value {numref}`fig:initial-bound-water`. This disequilibrium triggers
+continuous fluid release from the subducting plate, and is meant to represent the continuous input of the slab into the trench.
+
+```{figure-md} fig:initial-bound-water
+<img src="bound_water.png" />
+
+ The initial water content within the subducting plate. The model includes layered lithologies. From top to bottom: 10 km of sediment with 2 wt% H$_2$O, 10 km of MORB with 1 wt% H$_2$O, 10 km of gabbro with 0.5 wt% H$_2$O, and 90 km of peridotite with 1 wt% H$_2$O. The white contours show isotherms at 200 K intervals spanning 300 K to 1300 K, and the black contours show depths at 100 km intervals.
+```
+
+Once water is released from the subducting plate, it is advected according to either the Darcy velocity or the fluid velocity (from the fully coupled
+McKenzie equations (to be implemented)). The solid velocity is still computed each time step, as the presence of bound or free fluids act to reduce
+the solid viscosity and thereby impact the solid velocity. Since the fluid velocity depends on the solid velocity in both advection cases, updating
+the solid velocity every time step is important, even if we do not advect the solid in this setup {numref}`fig:velocity-temperature`.
+The total solid viscosity is determined via:
+
+```{math}
+:label: eq:creep-viscosity
+\eta_i = \frac{1}{2}A_i^{-\frac{1}{n_i}} \dot{\epsilon}_{ii}^{\frac{1 - n_i}{n_i}} d^{m_i} f_{H_2O}^{-r_i} \text{exp}\left(-\alpha \phi \right) \text{exp}\left( \frac{E_i + PV_i}{n_i RT} \right)
+```
+
+```{figure-md} fig:velocity-temperature
+<img src="solid_velocity.png" />
+
+ The temperature of the model within the immediate vicinity of the subduction zone. The vectors show the direction of and are scaled by the velocity of the solid phase. The cool, dense slab generates corner flow style convection patterns in the mantle wedge.
+```
+
+where $i$ denote the deformation mechanism (diffusion or dislocation creep), $n_i$ is the stress exponent, $m_i$ is the grain size exponent, $r_i$ is
+the water fugacity exponent, $A_i$ is a prefactor, $E_i$ is the activation energy, $V_i$ is the activation volume, $d$ is the grain size, $P$ is
+pressure, $T$ is temperature, and $R$ is the gas constant. The presence of water reduces the viscosity via the terms
+$\text{exp}\left(-\alpha \phi \right)$, where $\alpha$ is an exponent constant and $\phi$ is the volume fraction of free water,
+and through $f_{H_2O}$, which is the water fugacity determined based on the amount of water within the solid {numref}`fig:viscosity-comparison`. In
+the input file, enabling this behaviour looks like this:
+```{literalinclude} input_water_viscous_weakening.part.prm
+```
+
+```{figure-md} fig:viscosity-comparison
+<img src="viscosity_profiles.png" />
+
+ The viscosity of the model within the immediate vicinity of the subduction zone just after the beginning of the model run-time (top) and at the end of the model run-time (bottom). The viscosity within the mantle wedge and the overriding plate is notably reduced due to the presence of free and bound-fluid (see {numref}`fig:fluid-pathways`)
+```
+
+In these models, the water density (1000 kg/m$^3$) is significantly lower than the solid density (3300 kg/m$^3$), so the fluid velocity is dominantly
+vertical due to the high buoyancy force experienced by the fluid. However, corner flow in the mantle wedge does impose a trench-ward horizontal
+component to the fluid velocity. As the fluid ascends through the hot peridotite mantle wedge, the PT conditions do not allow for the free fluid to
+be reabsorbed into the solid phase {numref}`fig:fluid-pathways`. However, when the free fluid starts ascending through the cooler overriding plate,
+the temperature is sufficiently low that hydration of the overriding plate begins to occur {numref}`fig:fluid-pathways`. The free fluid in the mantle
+wedge leads to a reduction of wedge viscosities, and the combination of free fluid and bound fluid in the overriding plate leads localized reductions
+of the overriding plate viscosity {numref}`fig:viscosity-comparison`. Within the model where the fluid is advected with the Darcy velocity, two
+distinct bands of free water can be seen seeping out of the subducting plate. The first band produces a larger flux of water, is centered around a
+distance of ~175 km landward of the trench, and is sourced from the subducting peridotite layer. The second band is lower in flux magnitude, is
+centered around a distance of ~400 km landward from the trench, and is sourced from the gabbro and MORB layers.
+
+```{figure-md} fig:fluid-pathways
+<img src="fluid_figures.png" />
+
+ The viscosity distribution near the start of the model run (top) and at end of the model run (bottom), with the free fluid overlain on top. White contour shows where the bound water makes up 0.5 wt% of the solid.
+```
+
+## Extending the Model
+There are several parameters which will dictate the magnitude of fluid flux and the fluid pathways within this cookbook. To investigate the
+importance of these parameters, here are some suggestions for how to tweak the .prm files:
+
+-   Change the `Disequilibrium percentage` parameter. The value used in this cookbook is 5%, and increasing this value will lead to a larger flux of fluid off the subducting plate.
+-   Change the `Fluid reaction time scale for operator splitting` parameter. The value used in this cookbook is 10,000 years, and increasing this will lead to a decrease in the fluid flux off the subducting plate.
+-   Changing the `Reference permeability` parameter. Decreasing the reference permeability will decrease the magnitude of the vertical fluid velocity, leading to more lateral advection of the fluid within the mantle wedge as it ascends.

cookbooks/2D_subduction_with_two_phase_flow/doc/input_water_viscous_weakening.part.prm

@@ -0,0 +1,20 @@
+subsection Material model
+  subsection Reactive Fluid Transport Model
+    # This parameter sets alpha and controls weakening based on free water
+    set Exponential fluid weakening factor            = 30
+    set Base model                                    = visco plastic
+  end
+
+ subsection Visco Plastic
+   # Setting the viscosity prefactor scheme to HK04 olivine hydration enables weakening of the solid
+   # when bound water is present
+   set Viscosity prefactor scheme                     = HK04 olivine hydration
+
+   set Stress exponents for diffusion creep           = 1.0
+   set Stress exponents for dislocation creep         = 3.5
+
+   # These parameters control the amount of weakening, with higher exponents resulting in a larger weakening
+   set Water fugacity exponents for diffusion creep   = 0.7 # 0.7 / 1.0 = 0.7
+   set Water fugacity exponents for dislocation creep = 0.34285714285714286 # 1.2 / 3.5 = 0.34285714285714286
+  end
+end

cookbooks/2D_subduction_with_two_phase_flow/doc/input_world_builder.part.prm

@@ -0,0 +1,13 @@
+set World builder file = world_builder_files/static_2D_cross_section_MOD.wb
+
+subsection Initial temperature model
+  set Model name = world builder
+end
+
+subsection Initial composition model
+  set List of model names = world builder disequilibrium
+  subsection World builder
+    set List of relevant compositions = porosity, bound_fluid, peridotite, gabbro, MORB, sediment
+    set Disequilibrium percentage     = 5
+  end
+end

cookbooks/2D_subduction_with_two_phase_flow/fixed_slab.wb

@@ -0,0 +1,95 @@
+// Composition 1 is bound fluid
+// Composition 2 is peridotite
+// Composition 3 is gabbro
+// Composition 4 is MORB
+// Composition 5 is sediment
+
+{
+  "version": "1.0",
+  "gravity model":{"model":"uniform", "magnitude":9.81},
+  "cross section":[[0,0],[50e3,0]],
+  "surface temperature":273,
+  "thermal expansion coefficient":3.0e-5, 
+  "thermal diffusivity":1.0e-6,
+  "potential mantle temperature":1573,
+  "features":
+  [
+    // Define a mantle layer, which has a uniform temperature of 1573 K. Following the order of compositional fields
+    // listed in the cookbook PRM file, the mantle layer is defined as the third compositional field (index 2), 
+    // which represents peridotite.
+    {"model": "mantle layer", "name": "peridotite mantle", "coordinates": [[-500e3, 100e3], [-500e3, -100e3], [16000e3, -100e3], [16000e3, 100e3]],
+     "min depth": 0, "max depth":10000e3, 
+     "composition models": 
+     [{"model": "uniform", "min depth":0.0, "max depth":10000e3, "compositions": [2]}],
+     "temperature models":
+     [{"model": "uniform", "min depth":0.0, "max depth":10000e3, "temperature": 1573}]},
+
+     // Define an overriding plate, which is 80 km thick, 2500 km long, has a linear temperature profile, and is assigned the peridotite composition
+     {"model": "oceanic plate", "name": "Overriding Plate",
+     "coordinates": [[4000e3, -100e3], [4000e3, 100e3], [7500e3, 100e3], [7500e3, -100e3]],
+     "min depth": 0,
+     "max depth": 80e3,
+     "composition models": [{"model": "uniform", "compositions": [2], "min depth":0, "max depth": 100e3}],
+     "temperature models": [{"model": "linear", "top temperature":273, "bottom temperature":1573, "min depth": 0, "max depth":80e3}]
+   },
+
+    // Define the unsubducted portion of the subducting plate. The plate is 3000 km long, 120 km thick, has a plate cooling geotherm,
+    // and is assigned a layered composition. From top to bottom, the layers are a 10 km thick sediment, a 10 km thick mid-ocean ridge
+    // basalt, a 10 km thick gabbro, and a 90 km thick peridotite, with and initial bound water equal to 2 wt%, 1 wt%, 0.5 wt%, and 1 wt%,
+    // respectively.
+    {"model": "oceanic plate", "name": "Subducting Plate",
+      "coordinates": [[4000e3, -100e3], [4000e3, 100e3], [1000e3, 100e3], [1000e3, -100e3]],
+      "min depth": 0,
+      "max depth": 300e3,
+      "composition models": 
+      [
+        {"model": "tian water content", "compositions": [1], "min depth":0.0, "max depth": 10e3, "lithology": "sediment", "initial water content": 2, "cutoff pressure": 1},
+        {"model": "tian water content", "compositions": [1], "min depth":10e3, "max depth": 20e3, "lithology": "MORB", "initial water content": 1, "cutoff pressure": 16},
+        {"model": "tian water content", "compositions": [1], "min depth":20e3, "max depth": 30e3, "lithology": "gabbro", "initial water content": 0.5, "cutoff pressure": 26},
+        {"model": "tian water content", "compositions": [1], "min depth":30e3, "max depth": 120e3, "lithology": "peridotite", "initial water content": 1, "cutoff pressure": 10},
+        {"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"},
+        {"model": "uniform", "compositions": [2], "min depth":30e3, "max depth": 120e3, "operation":"add"}],
+
+      "temperature models": [{"model": "plate model", "ridge coordinates": [[[1000e3,-100e3],[1000e3,100e3]]],
+      "spreading velocity": 0.03
+      , "min depth": 0, "max depth": 300e3, "bottom temperature": 1573}]
+    },
+
+    // Define the subducted portion of the subducting plate. The slab has the same compositional properties as the 
+    // unsubducted portion defined in the feature above. The slab is 1100 km long in total, first curving to a dip 
+    // of 45 degrees over 300 km, and then continuing down at a constant 45 degree dip for another 800 km. The 
+    // temperature is defined using a plate cooling model for the base of the plate, and an infinite half-space
+    // model for the top of the plate. 
+    {"model":"subducting plate", "name":"Slab",   
+
+      "coordinates":[[4000e3, -100e3],[4000e3, 100e3]], 
+      "dip point":[1e10,0],"max depth":1000e3,
+	    "segments":[{"length":300e3,"thickness":[300e3],"top truncation":[-50e3],"angle":[0,45]},
+                        {"length":500e3,"thickness":[300e3],"top truncation":[-50e3],"angle":[45,45]}, 
+                        {"length":300e3,"thickness":[300e3],"top truncation":[-50e3],"angle":[45,45]}],
+
+	    "composition models":[
+        {"model":"tian water content", "compositions":[1], "min distance slab top":0, "max distance slab top":10e3, "lithology":"sediment", "initial water content": 2, "cutoff pressure": 1},
+        {"model":"tian water content", "compositions":[1], "min distance slab top":10e3, "max distance slab top":20e3, "lithology":"MORB", "initial water content": 1, "cutoff pressure": 16},
+        {"model":"tian water content", "compositions":[1], "min distance slab top":20e3, "max distance slab top":30e3, "lithology":"gabbro", "initial water content": 0.5, "cutoff pressure": 26},
+        {"model":"tian water content", "compositions":[1], "min distance slab top":30e3, "max distance slab top":120e3, "lithology":"peridotite", "initial water content": 1, "cutoff pressure": 10},
+        {"model":"uniform", "compositions":[5], "min distance slab top":0.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"},
+        {"model":"uniform", "compositions":[2], "min distance slab top":30e3, "max distance slab top":120e3, "operation":"add"}],
+
+        "temperature models":[{"model":"mass conserving", 
+	      "reference model name": "plate model",
+		    "density":3300, "thermal conductivity":3.3, "adiabatic heating":false,
+		    "subducting velocity":0.03,
+                    "spreading velocity":0.03,
+		    "ridge coordinates":[[[1000e3,-100e3],[1000e3,100e3]]],
+	 	    "coupling depth":80e3, 
+		    "forearc cooling factor":1.0,
+		    "taper distance":150e3, 
+	 	    "min distance slab top":-200e3, "max distance slab top":300e3}] 
+	 }
+  ]
+}