add more doc

Author: danieldouglas92

cookbooks/2D_subduction_with_two_phase_flow/doc/2D_subduction_two_phase_flow.md

@@ -5,32 +5,45 @@
 
 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.
 
-The cookbook requires that ASPECT is compiled with the Geodynamic WorldBuilder (GWB), which can be done by setting `ASPECT_WITH_WORLD_BUILDER=ON` when configuring ASPECT with CMAKE (notably, this is the default setting). The GWB is a powerful tool which allows ASPECT users to generate complicated initial conditions, and we will use it to generate the temperature and hydration state of a two-dimensional subduction zone. Using ASPECT with the GWB requires specifying the path to a worldbuilder (.wb) file in the ASPECT input file, as well as specifying that the initial temperature and composition is being generated with the GWB. These lines look like this:
+
+
+This cookbook requires that ASPECT is compiled with the Geodynamic WorldBuilder (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 WorldBuilder (.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} fixed_slab.part.prm
 ```
 
-The worldbuilder 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 worldbuilder 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 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 worldbuilder 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 worldbuilder 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).
+
+```{figure-md} fig:model-overview
+<img src="model_overview.png" />
+
+ The model geometry coloured by the model temperature.
+```
+
+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 mantle a temperature of 1573 K. 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.1, resulting in a 10% 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" />
 
- Test1
+ The initial water content within the subducting plate. White contours shown isotherms at 200 K intervals from 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). The solid velocity is still computed each time step, as the presence of fluid--either bound within the solid or as free water--reduce the solid viscosity which can 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. 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)
+```
+
+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. In these models, the water density (1000 kg/m$^3$) is significantly lower than the solid density (3000 kg/m$^3$), so the fluid velocity is dominantly vertical due to the high buoyancy force exerted on the fluid. However, corner flow in the mantle wedge does impose a trench-ward horizontal component to the fluid velocity.
+
 ```{figure-md} fig:final-water
 <img src="final_water.png" />
 
- Test1
+ The viscosity distribution at the end of the model run, with the free fluid overlain on top.
 ```
 
 ```{figure-md} fig:initial-viscosity
 <img src="slab_viscosity.png" />
 
- Test1
-```
-
-```{figure-md} fig:model-overview
-<img src="model_overview.png" />
-
- Test1
+ The viscosity 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.
 ```

cookbooks/2D_subduction_with_two_phase_flow/fixed_slab.wb

@@ -14,13 +14,15 @@
   "potential mantle temperature":1573,
   "features":
   [
+    // Define a mantle layer, which has a uniform temperature of 1573 K. The mantle layer has a composition 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,
@@ -29,6 +31,10 @@
      "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,
@@ -49,6 +55,11 @@
       , "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]], 

cookbooks/2D_subduction_with_two_phase_flow/fixed_slab_darcy_with_fugacity.prm

@@ -12,7 +12,7 @@
 # In this parameter file, when fluid is contained within the solid phase the viscosity is reduced
 # according to a calculated water fugacity term. This water fugacity term assumes a uniform olivine
 # composition of 90% forsterite and 10% fayalite throughout the entire model domain. While this is
-# an oversimplification, it is sufficient for demonstrating the functionaly within the scope of
+# an oversimplification, it is sufficient for demonstrating the functionaity within the scope of
 # this cookbook.
 
 # Load the plugin library and the world builder file, which are both required for this cookbook

cookbooks/2D_subduction_with_two_phase_flow/fixed_slab_mckenzie_with_fugacity.prm

@@ -33,6 +33,6 @@ end
 # Remove the darcy velocity from the list of output variables.
 subsection Postprocess
   subsection Visualization
-    set List of output variables      = material properties, strain rate, named additional outputs
+    set List of output variables      = material properties, strain rate, named additional outputs, melt material properties
   end
 end