diff --git a/cookbooks/2D_subduction_with_two_phase_flow/doc/2D_subduction_two_phase_flow.md b/cookbooks/2D_subduction_with_two_phase_flow/doc/2D_subduction_two_phase_flow.md
new file mode 100644
index 00000000000..998a4eab27d
--- /dev/null
+++ b/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
+
+
+ The model domain coloured by the model temperature. The subducting plate descends into the mantle at a constant dip of 45°. 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° over a (slab) length of 300 km, then continues into the mantle at a constant 45° 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
+
+
+ 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
+
+
+ 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
+
+
+ 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
+
+
+ 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.
diff --git a/cookbooks/2D_subduction_with_two_phase_flow/doc/bound_water.png b/cookbooks/2D_subduction_with_two_phase_flow/doc/bound_water.png
new file mode 100644
index 00000000000..b864d2a5fa5
Binary files /dev/null and b/cookbooks/2D_subduction_with_two_phase_flow/doc/bound_water.png differ
diff --git a/cookbooks/2D_subduction_with_two_phase_flow/doc/fluid_figures.png b/cookbooks/2D_subduction_with_two_phase_flow/doc/fluid_figures.png
new file mode 100644
index 00000000000..5f9529eb1b2
Binary files /dev/null and b/cookbooks/2D_subduction_with_two_phase_flow/doc/fluid_figures.png differ
diff --git a/cookbooks/2D_subduction_with_two_phase_flow/doc/input_water_viscous_weakening.part.prm b/cookbooks/2D_subduction_with_two_phase_flow/doc/input_water_viscous_weakening.part.prm
new file mode 100644
index 00000000000..1342a38ee44
--- /dev/null
+++ b/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
diff --git a/cookbooks/2D_subduction_with_two_phase_flow/doc/input_world_builder.part.prm b/cookbooks/2D_subduction_with_two_phase_flow/doc/input_world_builder.part.prm
new file mode 100644
index 00000000000..8bff7d118cc
--- /dev/null
+++ b/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
diff --git a/cookbooks/2D_subduction_with_two_phase_flow/doc/model_overview.png b/cookbooks/2D_subduction_with_two_phase_flow/doc/model_overview.png
new file mode 100644
index 00000000000..465b814d9ee
Binary files /dev/null and b/cookbooks/2D_subduction_with_two_phase_flow/doc/model_overview.png differ
diff --git a/cookbooks/2D_subduction_with_two_phase_flow/doc/solid_velocity.png b/cookbooks/2D_subduction_with_two_phase_flow/doc/solid_velocity.png
new file mode 100644
index 00000000000..c589fbc66b9
Binary files /dev/null and b/cookbooks/2D_subduction_with_two_phase_flow/doc/solid_velocity.png differ
diff --git a/cookbooks/2D_subduction_with_two_phase_flow/doc/viscosity_profiles.png b/cookbooks/2D_subduction_with_two_phase_flow/doc/viscosity_profiles.png
new file mode 100644
index 00000000000..af8b6fabfec
Binary files /dev/null and b/cookbooks/2D_subduction_with_two_phase_flow/doc/viscosity_profiles.png differ
diff --git a/cookbooks/2D_subduction_with_two_phase_flow/fixed_slab.wb b/cookbooks/2D_subduction_with_two_phase_flow/fixed_slab.wb
new file mode 100644
index 00000000000..18e3b3127bd
--- /dev/null
+++ b/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}]
+ }
+ ]
+}
diff --git a/cookbooks/2D_subduction_with_two_phase_flow/fixed_slab_darcy_with_fugacity.prm b/cookbooks/2D_subduction_with_two_phase_flow/fixed_slab_darcy_with_fugacity.prm
new file mode 100644
index 00000000000..bd71b228229
--- /dev/null
+++ b/cookbooks/2D_subduction_with_two_phase_flow/fixed_slab_darcy_with_fugacity.prm
@@ -0,0 +1,342 @@
+# This parameter file is for a cookbook which simulates fluid migration within a 2D subduction
+# system. The initial conditions (for both temperature and composition) are defined using the
+# Geodynamic WorldBuilder. The initial composition for the world builder is modified to manually
+# increase the bound water content within the subducting plate by 10%, leaving the bound fluid in
+# disequilibrium. However, this cookbook prevents the bound fluid from decreasing within the
+# subducting plate, creating a fixed fluid source that is meant to mimic the steady input of water
+# into the subduction factory by continuous subduction. The free fluid released by the fixed source
+# is advected through the system according to the Darcy velocity, and is allowed to be reabsorbed by
+# the solid phase if the reactions determined by the Tian 2019 reaction model predict that this is
+# stable.
+#
+# 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 functionaity within the scope of
+# this cookbook.
+
+# Load the plugin library and the world builder file, which are both required for this cookbook
+set Additional shared libraries = plugin/lib2D_subduction_with_two_phase_flow.so
+set World builder file = $ASPECT_SOURCE_DIR/cookbooks/2D_subduction_with_two_phase_flow/fixed_slab.wb
+set Nonlinear solver scheme = iterated Advection and Stokes
+set Nonlinear solver failure strategy = continue with next timestep
+set Output directory = output-fixed_slab_darcy_fluid_advection
+# Use a loose nonlinear tolerance of 1e-3 so that this cookbook can run with a relatively low
+# computational expense.
+set Max nonlinear iterations = 100
+set Nonlinear solver tolerance = 1e-3
+set Dimension = 2
+set End time = 1e6
+set CFL number = 0.5 # 0.5 for stability
+set Pressure normalization = surface
+set Adiabatic surface temperature = 1573
+set Surface pressure = 0
+
+# Limit the first time step to prevent numerical instabilities
+set Maximum first time step = 100
+
+# The Tian 2019 reaction model requires the use of Operator splitting.
+set Use operator splitting = true
+# Automatically resume computation, if a restart file is present.
+set Resume computation = auto
+
+# We use the block AMG solver instead of block GMG to more directly compare the
+# behaviour of this model with the model that advects the fluid according to the
+# fully coupled McKenzie equations, which does not support the use of block GMG.
+# Additionally, the Tian 2019 reaction model requires operator splitting with a
+# fixed reaction solver type.
+subsection Solver parameters
+ subsection Stokes solver parameters
+ set Stokes solver type = block AMG
+ set GMRES solver restart length = 1000
+ set Number of cheap Stokes solver steps = 20000
+ set Use full A block as preconditioner = true
+ set Linear solver tolerance = 1e-5
+ set Maximum number of expensive Stokes solver steps = 0
+ end
+ subsection Operator splitting parameters
+ set Reaction solver type = fixed step
+ set Reaction time step = 50
+ end
+end
+
+# Increase the smoothing parameters beta and cR for numerical stability and
+# computational efficiency.
+subsection Discretization
+ subsection Stabilization parameters
+ set beta = 0.52
+ set cR = 1.1
+ end
+end
+
+# Checkpoint the model every 50 time steps
+subsection Checkpointing
+ set Steps between checkpoint = 50
+end
+
+# The Tian 2019 reaction model requires 6 compositional fields:
+# porosity, bound_fluid, peridotite, gabbro, MORB, and sediment.
+# We want to keep all of the solid compositions static, and only
+# advect the porosity composition, therefore we set the solid field
+# advection methods to static and the porosity to darcy field.
+subsection Compositional fields
+ set Number of fields = 6
+ set Names of fields = porosity, bound_fluid, peridotite, gabbro, MORB, sediment
+ set Compositional field methods = darcy field, static, static, static, static, static
+ set Types of fields = porosity, generic, chemical composition, chemical composition, chemical composition, chemical composition
+end
+
+# We also want to keep the temperature distribution fixed.
+subsection Temperature field
+ set Temperature method = static
+end
+
+# The initial composition is defined using a modified version of the world builder
+# initial composition model compiled for this plugin. The `world builder disequilibrium`
+# initial composition model works exactly like the regular world builder initial composition
+# model, except that it increases the value of the bound_fluid composition by 10%.
+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
+
+# The initial temperature is also defined using the world builder.
+subsection Initial temperature model
+ set Model name = world builder
+end
+
+# The geometry is a very large 2D box, and the cells have an aspect ratio of 1:1.
+# The geometry is large to ensure that the model boundaries are not significantly
+# impacting the velocities within the vicinty of the slab and in the mantle wedge.
+# Prior to additional refinements, the base grid size is 1450 x 1450 km.
+subsection Geometry model
+ set Model name = box
+ subsection Box
+ set X extent = 8700e3
+ set Y extent = 2900e3
+
+ set X repetitions = 6
+ set Y repetitions = 2
+
+ set Box origin X coordinate = 0
+ set Box origin Y coordinate = 0
+ end
+end
+
+# Downward gravity
+subsection Gravity model
+ set Model name = function
+ subsection Function
+ set Variable names = x,y
+ set Function expression = 0; -9.81
+ end
+end
+
+# The model is incompressible.
+subsection Formulation
+ set Formulation = Boussinesq approximation
+end
+
+# Fix the temperature on the top and the bottom of the model.
+subsection Boundary temperature model
+ set Fixed temperature boundary indicators = top, bottom
+ set List of model names = box
+ subsection Box
+ set Bottom temperature = 1573
+ set Top temperature = 273
+ end
+end
+
+# Reactions are computed using the Tian 2019 reaction model, and are composited
+# with a visco plastic solid rheology. The Tian 2019 reaction model is modified
+# within this cookbook to prevent the bound_fluid composition from ever decreasing.
+subsection Material model
+
+ set Model name = reactive fluid transport bound fluid source
+ # We use a geometric averaging scheme for the viscosity instead of harmonic to
+ # avoid oversmoothing of the viscosity.
+ set Material averaging = geometric average only viscosity
+
+ subsection Reactive Fluid Transport Model
+ set Base model = visco plastic
+ set Reference fluid density = 1000 # kg/m^3
+ set Shear to bulk viscosity ratio = 0.1
+ set Reference fluid viscosity = 1 # Pa s
+ set Reference permeability = 5e-8
+ # This parameter controls how much the viscosity is reduced in the presence
+ # of free water.
+ set Exponential fluid weakening factor = 30
+ set Fluid compressibility = 0 # incompressible fluid
+ # This parameter controls how quickly the water is released or absorbed by the
+ # solid phase. A time scale of 10,000 years ensures that these reactions are
+ # happening rapdily relative to the timescale of solid advection.
+ set Fluid reaction time scale for operator splitting = 1e4
+ # We use the Tian 2019 reaction model, which parameterizes the solubility using
+ # polynomials for the four solid compositions used in this cookbook.
+ set Fluid-solid reaction scheme = tian approximation
+ subsection Tian 2019 model
+ # The polynomials for each lithology blow up to extremely large values at low
+ # temperatures, we cap the bound water at these values to prevent this.
+ set Maximum weight percent water in peridotite = 4
+ set Maximum weight percent water in gabbro = 2
+ set Maximum weight percent water in MORB = 3
+ set Maximum weight percent water in sediment = 4
+ end
+ end
+
+ # This subsection defines how the solid deforms using a composite visco-plastic
+ # rheology. The rheology uses dislocation creep, diffusion creep, drucker-prager
+ # plasticity, and introduces a viscous weakening term that is dependent on the
+ # modeled bound water content.
+ subsection Visco Plastic
+ # Use a geometric averaging scheme for the viscosity
+ set Viscosity averaging scheme = geometric
+ set Viscous flow law = composite
+ # The HK04 olivine hydration prefactor scales the solid viscosity based on the
+ # modeled water-fugacity, assuming an olivine composition that is 90% forsterite
+ # and 10% fayalite. This implementation is based off of Hirth & Kohlstaedt 2004.
+ set Viscosity prefactor scheme = HK04 olivine hydration
+ set Minimum strain rate = 1e-20
+ set Reference strain rate = 1e-20
+
+ # We introduce phase transitions as a method for changing the rheology of the model
+ # within the lower mantle, and to limit the viscosity of the sediment rheology. Within
+ # the lower mantle (beneath 660 km depth), the solid only deforms via diffusion creep.
+ # Shallower than 660 km, all deformation mechanisms are active. For the sediment
+ # composition only, the viscosity is prescribed to be 1e19 Pa s, at depths shallower
+ # than 125 km.
+ set Phase transition depths = background: 660e3, porosity: 660e3, bound_fluid: 660e3, peridotite: 660e3, gabbro: 660e3, MORB: 660e3, sediment: 125e3|660e3
+ set Phase transition widths = background: 1e3, porosity: 1e3, bound_fluid: 1e3, peridotite: 1e3, gabbro: 1e3, MORB: 1e3, sediment: 10e3|1e3
+ set Phase transition temperatures = background: 1573, porosity: 1573, bound_fluid: 1573, peridotite: 1573, gabbro: 1573, MORB: 1573, sediment: 1573|1573
+ set Phase transition Clapeyron slopes = background: 0, porosity: 0, bound_fluid: 0, peridotite: 0, gabbro: 0, MORB: 0, sediment: 0|0
+ set Densities = background: 3300, porosity: 3300, bound_fluid: 3300, peridotite: 3300, gabbro: 3300, MORB: 3300, sediment: 3300
+ set Minimum viscosity = background: 2.5e18|2.5e18, \
+ porosity: 2.5e18|2.5e18, \
+ bound_fluid: 2.5e18|2.5e18, \
+ peridotite: 2.5e18|2.5e18, \
+ gabbro: 2.5e18|2.5e18, \
+ MORB: 2.5e18|2.5e18, \
+ sediment: 9.999e18|2.5e18|2.5e18
+
+ set Maximum viscosity = background: 2.5e23|2.5e23, \
+ porosity: 2.5e23|2.5e23, \
+ bound_fluid: 2.5e23|2.5e23, \
+ peridotite: 2.5e23|2.5e23, \
+ gabbro: 2.5e23|2.5e23, \
+ MORB: 2.5e23|2.5e23, \
+ sediment: 1.001e19|2.5e23|2.5e23
+
+ # Diffusion creep parameters from Hirth & Kohlstaedt 2004
+ set Prefactors for diffusion creep = 1.577e-21
+
+ set Water fugacity exponents for diffusion creep = 0.7
+
+ set Stress exponents for diffusion creep = 1.0
+
+ set Activation volumes for diffusion creep = 10e-6
+
+ set Grain size = 1e-3
+
+ set Activation energies for diffusion creep = 375e3
+
+ # Dislocation creep parameters from Hirth & Kohlstaedt 2004, disable dislocation creep
+ # within the lower mantle (beneath 660 km depth).
+ set Prefactors for dislocation creep = background: 1.01e-25|1e-50, \
+ porosity: 1.01e-25|1e-50, \
+ bound_fluid: 1.01e-25|1e-50, \
+ peridotite: 1.01e-25|1e-50, \
+ gabbro: 1.01e-25|1e-50, \
+ MORB: 1.01e-25|1e-50, \
+ sediment: 1.01e-25|1.01e-25|1e-50
+
+ set Stress exponents for dislocation creep = 3.5
+
+ set Water fugacity exponents for dislocation creep = 0.34285714285714286
+
+ set Activation volumes for dislocation creep = 22e-6
+
+ set Activation energies for dislocation creep = 520e3
+
+ # Drucker-Prager plasticity parameters
+ set Angles of internal friction = 30
+
+ set Cohesions = 60e6
+
+ set Maximum yield stress = 500e6
+
+ # For stability and efficiency, do not use the full pressure.
+ set Use adiabatic pressure in plasticity = true
+
+ # Thermal parameters.
+ set Thermal expansivities = 3.0e-5
+ set Thermal diffusivities = 1e-6
+ set Heat capacities = 750
+ end
+end
+
+# Adapt the mesh using a combination of strategies. Every 5 time steps,
+# the mesh is refined based on the temperature, viscosity, and the location
+# of free water. We ensure that mesh cannot go below 4 levels of refinement
+# using a minimum refinement function. The highest refinement level is 8
+# which results in a grid size of ~2.75 x 2.75 km.
+subsection Mesh refinement
+ set Coarsening fraction = 0.3
+ set Refinement fraction = 0.5
+ set Initial adaptive refinement = 4
+ set Initial global refinement = 5
+ set Strategy = temperature, viscosity, composition threshold, minimum refinement function
+ set Time steps between mesh refinement = 5
+ set Refinement criteria scaling factors = 2, 1.5, 1.5, 1
+ set Refinement criteria merge operation = max
+ set Run postprocessors on initial refinement = false
+ set Skip solvers on initial refinement = true
+
+ # minimum of 6 levels of refinement within the lithosphere, 5 levels
+ # of refinement at the phase transition, and 4 levels of refinement
+ # everywhere else.
+ subsection Minimum refinement function
+ set Function constants = litho_thickness=120e3, ymax=2900e3
+ set Coordinate system = cartesian
+ set Variable names = x, y
+ set Function expression = if(y>(ymax-litho_thickness),6, \
+ if(y>=(ymax - 662e3) && y<=(ymax - 658e3), 5, 4))
+ end
+
+ # refine where the porosity is bigger than 1e-6, the bound fluid is larger
+ # than 0.001, and where the sediment composition is greater than 0.5. Set
+ # the peridotite, gabbro, and MORB compositions to 10, which means that we
+ # will never refine based on these compositions since they cannot be larger
+ # than 1.0.
+ subsection Composition threshold
+ set Compositional field thresholds = 1e-6, 0.001, 10.0, 10.0, 10.0, 0.5
+ end
+end
+
+# Free slip boundaries on all sides.
+subsection Boundary velocity model
+ set Tangential velocity boundary indicators = left, right, bottom, top
+end
+
+# For this model, we advect the fluid according to the Darcy velocity, so we
+# disable melt transport.
+subsection Melt settings
+ set Include melt transport = false
+end
+
+# Postprocessing, output every 10,000 years.
+subsection Postprocess
+ set List of postprocessors = visualization, composition statistics, velocity statistics, fluid velocity statistics
+ subsection Visualization
+ set List of output variables = material properties, strain rate, named additional outputs, darcy velocity
+ set Output format = vtu
+ set Time between graphical output = 10e3
+ set Interpolate output = true
+ set Number of grouped files = 0
+ subsection Material properties
+ set List of material properties = density, viscosity
+ end
+ end
+end
diff --git a/cookbooks/2D_subduction_with_two_phase_flow/fixed_slab_mckenzie_with_fugacity.prm b/cookbooks/2D_subduction_with_two_phase_flow/fixed_slab_mckenzie_with_fugacity.prm
new file mode 100644
index 00000000000..4f44653302a
--- /dev/null
+++ b/cookbooks/2D_subduction_with_two_phase_flow/fixed_slab_mckenzie_with_fugacity.prm
@@ -0,0 +1,41 @@
+# This model uses the same foundation as the other parameter file used in this cookbook,
+# so we include this parameter file and only change what is needed. The major difference
+# between this parameter file and fixed_slab_darcy_with_fugacity.prm is that here the
+# fluid phase is advected according to the fully coupled McKenzie equations.
+include $ASPECT_SOURCE_DIR/cookbooks/2D_subduction_with_two_phase_flow/fixed_slab_darcy_with_fugacity.prm
+
+# Change the output directory.
+set Output directory = output-fixed_slab_mckenzie_fluid_advection
+
+# We need to make sure that the porosity field is no longer advected with the darcy velocity.
+# So change the advection method back to field.
+subsection Compositional fields
+ set Names of fields = porosity, bound_fluid, peridotite, gabbro, MORB, sediment
+ set Compositional field methods = field, static, static, static, static, static
+end
+
+# Setting Include melt transport = true allows the free fluid to be advected with the
+# McKenzie equations.
+subsection Melt settings
+ set Include melt transport = true
+end
+
+# When Include melt transport = true, the effects of compaction are taken into account
+# when determining the solid and fluid velocities. We limit the compaction viscosities
+# to improve numerical efficiency and stability.
+subsection Material model
+ subsection Reactive Fluid Transport Model
+ set Minimum compaction viscosity = 1e19
+ set Maximum compaction viscosity = 1e23
+ end
+end
+
+# Remove the darcy velocity from the list of output variables, output the melt material properties.
+subsection Postprocess
+ subsection Visualization
+ set List of output variables = material properties, strain rate, named additional outputs, melt material properties
+ subsection Melt material properties
+ set List of properties = compaction length, permeability, compaction viscosity
+ end
+ end
+end
diff --git a/cookbooks/2D_subduction_with_two_phase_flow/plugin/CMakeLists.txt b/cookbooks/2D_subduction_with_two_phase_flow/plugin/CMakeLists.txt
new file mode 100644
index 00000000000..46f90f20022
--- /dev/null
+++ b/cookbooks/2D_subduction_with_two_phase_flow/plugin/CMakeLists.txt
@@ -0,0 +1,11 @@
+CMAKE_MINIMUM_REQUIRED(VERSION 3.13.4)
+
+FIND_PACKAGE(Aspect REQUIRED HINTS ${Aspect_DIR} ../ ../../ $ENV{ASPECT_DIR})
+DEAL_II_INITIALIZE_CACHED_VARIABLES()
+
+SET(TARGET "2D_subduction_with_two_phase_flow")
+PROJECT(${TARGET})
+
+# ADD_LIBRARY(${TARGET} SHARED world_builder_too_high.cc fixed_reactions.cc )
+ADD_LIBRARY(${TARGET} SHARED bound_fluid_source.cc world_builder_disequilibrium_fluid.cc)
+ASPECT_SETUP_PLUGIN(${TARGET})
diff --git a/cookbooks/2D_subduction_with_two_phase_flow/plugin/bound_fluid_source.cc b/cookbooks/2D_subduction_with_two_phase_flow/plugin/bound_fluid_source.cc
new file mode 100644
index 00000000000..01eaba94ea7
--- /dev/null
+++ b/cookbooks/2D_subduction_with_two_phase_flow/plugin/bound_fluid_source.cc
@@ -0,0 +1,73 @@
+/*
+ Copyright (C) 2025 - by the authors of the ASPECT code.
+
+ This file is part of ASPECT.
+
+ ASPECT is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2, or (at your option)
+ any later version.
+
+ ASPECT is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with ASPECT; see the file LICENSE. If not see
+ .
+*/
+
+#include
+#include
+#include
+#include
+
+namespace aspect
+{
+ namespace MaterialModel
+ {
+ template
+ class BoundFluidSource: public MaterialModel::ReactiveFluidTransport
+ {
+ public:
+
+ virtual void evaluate(const MaterialModel::MaterialModelInputs &in,
+ MaterialModel::MaterialModelOutputs &out) const override;
+ };
+
+
+ template
+ void
+ BoundFluidSource::
+ evaluate(const MaterialModel::MaterialModelInputs &in,
+ MaterialModel::MaterialModelOutputs &out) const
+ {
+ ReactiveFluidTransport::evaluate(in, out);
+ const unsigned int bound_fluid_idx = this->introspection().compositional_index_for_name("bound_fluid");
+ const std::shared_ptr> reaction_rate_out
+ = out.template get_additional_output_object>();
+
+ if (this->get_parameters().use_operator_splitting && reaction_rate_out != nullptr
+ && in.requests_property(MaterialProperties::reaction_rates))
+ {
+ for (unsigned int q=0; q < in.n_evaluation_points(); ++q)
+ if (reaction_rate_out->reaction_rates[q][bound_fluid_idx] <= 0.0)
+ reaction_rate_out->reaction_rates[q][bound_fluid_idx] = 0.0;
+ }
+ }
+ }
+}
+
+// explicit instantiations
+namespace aspect
+{
+ namespace MaterialModel
+ {
+ ASPECT_REGISTER_MATERIAL_MODEL(BoundFluidSource,
+ "reactive fluid transport bound fluid source",
+ "A material model that is like the 'reactive "
+ "fluid transport' model, but modifies the reactions "
+ "to prevent the bound fluid content from decreasing.")
+ }
+}
diff --git a/cookbooks/2D_subduction_with_two_phase_flow/plugin/world_builder_disequilibrium_fluid.cc b/cookbooks/2D_subduction_with_two_phase_flow/plugin/world_builder_disequilibrium_fluid.cc
new file mode 100644
index 00000000000..37fd17ee4b8
--- /dev/null
+++ b/cookbooks/2D_subduction_with_two_phase_flow/plugin/world_builder_disequilibrium_fluid.cc
@@ -0,0 +1,225 @@
+/*
+ Copyright (C) 2025 - by the authors of the ASPECT code.
+
+ This file is part of ASPECT.
+
+ ASPECT is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2, or (at your option)
+ any later version.
+
+ ASPECT is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with ASPECT; see the file LICENSE. If not see
+ .
+*/
+
+#include
+
+#ifdef ASPECT_WITH_WORLD_BUILDER
+#include
+#include
+#include
+#include
+#include
+#include
+
+namespace WorldBuilder
+{
+ class World;
+}
+
+
+namespace aspect
+{
+ namespace InitialComposition
+ {
+ /**
+ * A class that implements initial conditions for the compositional fields
+ * based on a functional description provided in the input file through the
+ * World builder.
+ *
+ * @ingroup InitialCompositionModels
+ */
+ template
+ class WorldBuilderDisequilibrium : public Interface,
+ public SimulatorAccess
+ {
+ public:
+ /**
+ * Initialization function. This function is called once at the
+ * beginning of the program after parse_parameters is run and after
+ * the SimulatorAccess (if applicable) is initialized.
+ */
+ void
+ initialize () override;
+
+ /**
+ * Return the initial composition as a function of position and index
+ * of the compositional field. If c is equal to the index of the bound_fluid
+ * compositional field, then this function will increase the value of the
+ * bound_fluid field by a user defined percentage.
+ */
+ double initial_composition (const Point &position, const unsigned int c) const override;
+
+ /**
+ * Declare the parameters this class takes through input files.
+ */
+ static
+ void
+ declare_parameters (ParameterHandler &prm);
+
+ /**
+ * Read the parameters this class declares from the parameter file.
+ */
+ void
+ parse_parameters (ParameterHandler &prm) override;
+
+ private:
+ /**
+ * A vector that specifies for each compositional field index if the world builder
+ * should be evaluated for this compositional field.
+ */
+ std::vector compositional_field_defined_by_world_builder;
+
+ /**
+ * A pointer to the WorldBuilder object. Keeping this pointer ensures
+ * that the object doesn't go away while we still need it.
+ */
+ std::shared_ptr world_builder;
+
+ /**
+ * A value between 0 and 100, that defines the percent increase
+ * in the bound_fluid from it's equilibrium value.
+ */
+ double percentage_increase;
+ };
+ }
+}
+
+
+namespace aspect
+{
+ namespace InitialComposition
+ {
+ template
+ void
+ WorldBuilderDisequilibrium::
+ initialize()
+ {
+ CitationInfo::add("GWB");
+ world_builder = this->get_world_builder_pointer();
+ }
+
+
+
+ template
+ double
+ WorldBuilderDisequilibrium::
+ initial_composition (const Point &position, const unsigned int c) const
+ {
+ const unsigned int bound_fluid_idx = this->introspection().compositional_index_for_name("bound_fluid");
+ if (compositional_field_defined_by_world_builder[c] == true)
+ {
+ double comp_value = world_builder->composition(Utilities::convert_point_to_array(position),
+ -this->get_geometry_model().height_above_reference_surface(position),
+ c);
+ if (c == bound_fluid_idx &&
+ this->get_parameters().compositional_field_methods[bound_fluid_idx] == Parameters::AdvectionFieldMethod::static_field)
+ comp_value *= (1 + percentage_increase/100);
+ return comp_value;
+ }
+
+ return 0.0;
+ }
+
+
+
+ template
+ void
+ WorldBuilderDisequilibrium::declare_parameters (ParameterHandler &prm)
+ {
+ prm.enter_subsection("Initial composition model");
+ {
+ prm.enter_subsection("World builder");
+ {
+ prm.declare_entry ("List of relevant compositions", "",
+ Patterns::Anything(),
+ "A list of names of compositional fields for "
+ "which to determine the initial composition using "
+ "the World Builder. As World Builder evaluations can "
+ "be expensive, this parameter allows to only evaluate "
+ "the fields that are relevant. This plugin returns 0.0 "
+ "for all compositions that are not selected in the list. "
+ "By default the list is empty and the world builder is "
+ "evaluated for all compositional fields.");
+ prm.declare_entry ("Disequilibrium percentage", "5",
+ Patterns::Double(0),
+ "A percentage value that defines the increase in the "
+ "bound_fluid compositional field from its equilibrium value. "
+ "This is only applied if the bound_fluid compositional field "
+ "is defined by the World Builder. The default value is 10%, "
+ "which means that the initial value of the bound_fluid field "
+ "is 10 percent larger than its equilibrium value.");
+ }
+ prm.leave_subsection();
+ }
+ prm.leave_subsection();
+ }
+
+
+ template
+ void
+ WorldBuilderDisequilibrium::parse_parameters (ParameterHandler &prm)
+ {
+ prm.enter_subsection("Initial composition model");
+ {
+ prm.enter_subsection("World builder");
+ {
+ const std::vector composition_names = Utilities::split_string_list(prm.get("List of relevant compositions"));
+ percentage_increase = prm.get_double("Disequilibrium percentage");
+
+ compositional_field_defined_by_world_builder.resize(this->n_compositional_fields(),false);
+
+ if (composition_names.size() == 0)
+ for (unsigned int i=0; in_compositional_fields(); ++i)
+ compositional_field_defined_by_world_builder[i] = true;
+
+ for (const auto &composition_name: composition_names)
+ {
+ AssertThrow(this->introspection().compositional_name_exists (composition_name),
+ ExcMessage("All fields in \"List of relevant compositions\" must match names of compositional "
+ "fields as assigned in the \"Compositional fields/Names of fields\" parameter."));
+
+ compositional_field_defined_by_world_builder[this->introspection().compositional_index_for_name(composition_name)] = true;
+ }
+ }
+
+ prm.leave_subsection();
+ }
+ prm.leave_subsection();
+ }
+
+ }
+}
+
+
+// explicit instantiations
+namespace aspect
+{
+ namespace InitialComposition
+ {
+ ASPECT_REGISTER_INITIAL_COMPOSITION_MODEL(WorldBuilderDisequilibrium,
+ "world builder disequilibrium",
+ "Specify the initial composition through the World Builder. "
+ "If the composition is not named bound_fluid, the composition "
+ "is set to the value returned by the World Builder. If the "
+ "composition is named bound_fluid, the value returned by the "
+ "World Builder is increased by a user defined percentage. ")
+ }
+}
+#endif
diff --git a/doc/modules/changes/20250623_danieldouglas92 b/doc/modules/changes/20250623_danieldouglas92
new file mode 100644
index 00000000000..aef737def78
--- /dev/null
+++ b/doc/modules/changes/20250623_danieldouglas92
@@ -0,0 +1,4 @@
+Added: A cookbook which demonstrates two-phase reactive
+fluid transport in a fixed subduction system.
+
+(Daniel Douglas, 2025/06/23)
diff --git a/doc/sphinx/user/cookbooks/geophysical-setups.md b/doc/sphinx/user/cookbooks/geophysical-setups.md
index 20c994f2385..863afaff0a9 100644
--- a/doc/sphinx/user/cookbooks/geophysical-setups.md
+++ b/doc/sphinx/user/cookbooks/geophysical-setups.md
@@ -112,6 +112,7 @@ cookbooks/transform_fault_behn_2007/doc/transform_fault_behn_2007.md
cookbooks/kinematically_driven_subduction_2d/doc/kinematically_driven_subduction_2d.md
cookbooks/allken_et_al_2012_rift_interaction/doc/allken.md
cookbooks/tian_parameterization_kinematic_slab/doc/tian_parameterization_kinematic_slab.md
+cookbooks/2D_subduction_with_two_phase_flow/doc/2D_subduction_two_phase_flow.md
cookbooks/phase_transition_kinetics/doc/phase-transition-kinetics.md
cookbooks/mantle_convection_with_continents_in_annulus/doc/mantle_convection_in_annulus.md
cookbooks/inclusions/doc/inclusions.md
diff --git a/tests/2D_subduction_two_phase_flow_cookbook.cc b/tests/2D_subduction_two_phase_flow_cookbook.cc
new file mode 100644
index 00000000000..ff6017bde60
--- /dev/null
+++ b/tests/2D_subduction_two_phase_flow_cookbook.cc
@@ -0,0 +1,2 @@
+#include "../cookbooks/2D_subduction_with_two_phase_flow/plugin/bound_fluid_source.cc"
+#include "../cookbooks/2D_subduction_with_two_phase_flow/plugin/world_builder_disequilibrium_fluid.cc"
diff --git a/tests/2D_subduction_two_phase_flow_cookbook.prm b/tests/2D_subduction_two_phase_flow_cookbook.prm
new file mode 100644
index 00000000000..7bb48c12297
--- /dev/null
+++ b/tests/2D_subduction_two_phase_flow_cookbook.prm
@@ -0,0 +1,18 @@
+# This is a test that checks that the cookbook "2D_subduction_two_phase_flow" runs correctly. It also checks
+# that the two .cc files required for this cookbook work correctly by outputting the composition statistics.
+# The bound_fluid should be initialized with a max value of 0.021 (due to world_builder_disequilibrium_fluid.cc)
+# and this value should not decrease (due to bound_fluid_source.cc) even though the porosity becomes non-zero
+# after 100 years.
+include $ASPECT_SOURCE_DIR/cookbooks/2D_subduction_with_two_phase_flow/fixed_slab_darcy_with_fugacity.prm
+
+set End time = 100
+
+subsection Mesh refinement
+ set Initial adaptive refinement = 0
+ set Initial global refinement = 1
+ set Time steps between mesh refinement = 0
+end
+
+subsection Postprocess
+ set List of postprocessors = composition statistics
+end
diff --git a/tests/2D_subduction_two_phase_flow_cookbook/screen-output b/tests/2D_subduction_two_phase_flow_cookbook/screen-output
new file mode 100644
index 00000000000..d6dcc6392f8
--- /dev/null
+++ b/tests/2D_subduction_two_phase_flow_cookbook/screen-output
@@ -0,0 +1,63 @@
+
+Loading shared library <./lib2D_subduction_two_phase_flow_cookbook.debug.so>
+
+Number of active cells: 48 (on 2 levels)
+Number of degrees of freedom: 2,090 (450+65+225+225+225+225+225+225+225)
+
+*** Timestep 0: t=0 years, dt=0 years
+ Skipping porosity composition solve because RHS is zero.
+ Rebuilding Stokes preconditioner...
+ Solving Stokes system (AMG)... 19+0 iterations.
+ Relative nonlinear residuals (temperature, compositional fields, Stokes system): 0, 0, 0, 0, 0, 0, 0, 1
+ Relative nonlinear residual (total system) after nonlinear iteration 1: 1
+
+ Skipping porosity composition solve because RHS is zero.
+ Rebuilding Stokes preconditioner...
+ Solving Stokes system (AMG)... 16+0 iterations.
+ Relative nonlinear residuals (temperature, compositional fields, Stokes system): 0, 0, 0, 0, 0, 0, 0, 0.0325226
+ Relative nonlinear residual (total system) after nonlinear iteration 2: 0.0325226
+
+ Skipping porosity composition solve because RHS is zero.
+ Rebuilding Stokes preconditioner...
+ Solving Stokes system (AMG)... 15+0 iterations.
+ Relative nonlinear residuals (temperature, compositional fields, Stokes system): 0, 0, 0, 0, 0, 0, 0, 0.0122776
+ Relative nonlinear residual (total system) after nonlinear iteration 3: 0.0122776
+
+ Skipping porosity composition solve because RHS is zero.
+ Rebuilding Stokes preconditioner...
+ Solving Stokes system (AMG)... 14+0 iterations.
+ Relative nonlinear residuals (temperature, compositional fields, Stokes system): 0, 0, 0, 0, 0, 0, 0, 0.00470844
+ Relative nonlinear residual (total system) after nonlinear iteration 4: 0.00470844
+
+ Skipping porosity composition solve because RHS is zero.
+ Rebuilding Stokes preconditioner...
+ Solving Stokes system (AMG)... 13+0 iterations.
+ Relative nonlinear residuals (temperature, compositional fields, Stokes system): 0, 0, 0, 0, 0, 0, 0, 0.00179746
+ Relative nonlinear residual (total system) after nonlinear iteration 5: 0.00179746
+
+ Skipping porosity composition solve because RHS is zero.
+ Rebuilding Stokes preconditioner...
+ Solving Stokes system (AMG)... 12+0 iterations.
+ Relative nonlinear residuals (temperature, compositional fields, Stokes system): 0, 0, 0, 0, 0, 0, 0, 0.000685564
+ Relative nonlinear residual (total system) after nonlinear iteration 6: 0.000685564
+
+
+ Postprocessing:
+ Compositions min/max/mass: 0/0/0 // 0/0.021/8.585e+09 // 0/1/2.482e+13 // 0/0/0 // 0/0/0 // 0/1/4.088e+11
+
+*** Timestep 1: t=100 years, dt=100 years
+ Solving composition reactions... in 2 substep(s).
+ Solving porosity system ... 9 iterations.
+ Rebuilding Stokes preconditioner...
+ Solving Stokes system (AMG)... 10+0 iterations.
+ Relative nonlinear residuals (temperature, compositional fields, Stokes system): 0, 0.000159525, 0, 0, 0, 0, 0, 0.000385553
+ Relative nonlinear residual (total system) after nonlinear iteration 1: 0.000385553
+
+
+ Postprocessing:
+ Compositions min/max/mass: -6.174e-08/0.000693/2.448e+08 // 0/0.021/8.585e+09 // 0/1/2.482e+13 // 0/0/0 // 0/0/0 // 0/1/4.088e+11
+
+Termination requested by criterion: end time
+
+
+
diff --git a/tests/2D_subduction_two_phase_flow_cookbook/statistics b/tests/2D_subduction_two_phase_flow_cookbook/statistics
new file mode 100644
index 00000000000..43aa90574f3
--- /dev/null
+++ b/tests/2D_subduction_two_phase_flow_cookbook/statistics
@@ -0,0 +1,38 @@
+# 1: Time step number
+# 2: Time (years)
+# 3: Time step size (years)
+# 4: Number of mesh cells
+# 5: Number of Stokes degrees of freedom
+# 6: Number of temperature degrees of freedom
+# 7: Number of degrees of freedom for all compositions
+# 8: Number of nonlinear iterations
+# 9: Iterations for temperature solver
+# 10: Iterations for composition solver 1
+# 11: Iterations for composition solver 2
+# 12: Iterations for composition solver 3
+# 13: Iterations for composition solver 4
+# 14: Iterations for composition solver 5
+# 15: Iterations for composition solver 6
+# 16: Iterations for Stokes solver
+# 17: Velocity iterations in Stokes preconditioner
+# 18: Schur complement iterations in Stokes preconditioner
+# 19: Minimal value for composition porosity
+# 20: Maximal value for composition porosity
+# 21: Global mass for composition porosity
+# 22: Minimal value for composition bound_fluid
+# 23: Maximal value for composition bound_fluid
+# 24: Global mass for composition bound_fluid
+# 25: Minimal value for composition peridotite
+# 26: Maximal value for composition peridotite
+# 27: Global mass for composition peridotite
+# 28: Minimal value for composition gabbro
+# 29: Maximal value for composition gabbro
+# 30: Global mass for composition gabbro
+# 31: Minimal value for composition MORB
+# 32: Maximal value for composition MORB
+# 33: Global mass for composition MORB
+# 34: Minimal value for composition sediment
+# 35: Maximal value for composition sediment
+# 36: Global mass for composition sediment
+0 0.000000000000e+00 0.000000000000e+00 48 515 225 1350 6 0 0 0 0 0 0 0 89 95 188 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 2.10000000e-02 8.58520833e+09 0.00000000e+00 1.00000000e+00 2.48211806e+13 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 4.08819444e+11
+1 1.000000000000e+02 1.000000000000e+02 48 515 225 1350 1 0 9 0 0 0 0 0 10 11 22 -6.17357269e-08 6.93007567e-04 2.44781320e+08 0.00000000e+00 2.10000000e-02 8.58520833e+09 0.00000000e+00 1.00000000e+00 2.48211806e+13 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 4.08819444e+11