Add changelog, fix tests, address comments

cookbooks/tian_parameterization_kinematic_slab/coupled-two-phase-tian-parameterization-kinematic-slab.prm

@@ -4,8 +4,8 @@
 # are prescribed an initial hydration state after a period of equilibration, and
 # subsequently dehydrate as they subduct. This cookbook utlizes the existing melt framework
 # within ASPECT which governs the mechanical coupling of a solid and fluid phase (see
-# https://doi.org/10.1093/gji/ggw329) but introduces solid-fluid reactions
-# (i.e partitioning of fluid into/out of the solid) that are governed by an approximation
+# https://doi.org/10.1093/gji/ggw329) but introduces solid-volatile reactions
+# that are governed by an approximation
 # published by Tian et al., 2019 (https://doi.org/10.1029/2019GC008488).
 # In detail, the model simulates the upper 20 km of a subducting slab with a dip
 # of 45 degrees and a convergence rate of 5 cm/yr. The Cartesian box is rotated
@@ -15,17 +15,19 @@
 # that is already subducted. For stability, the slab first subducts for 100 kyr before
 # the hydration is added, and the slab surface temperature diffuses to the mantle over
 # 10 km perpendicular to the slab surface. From the slab surface to the base of the model,
-# the composition is comprised of a 5km thick layer of sediment with 3 wt% water,
-# below is a 7km thick layer of MORB with 2 wt% water, and below that is an
+# the composition is comprised of a 5 km thick layer of sediment with 3 wt% water,
+# below is a 7 km thick layer of MORB with 2 wt% water, and below that is an
 # 8 km thick layer of gabbro with 1 wt% water.
 
 set Adiabatic surface temperature              = 1600
 set Nonlinear solver scheme                    = iterated Advection and Stokes
 # We choose a low number of iterations and less restrictive nonlinear solver tolerance
 # for model efficiency. During the initial pulse of dehydration (for 3 time steps),
 # solvers do not quite converge to 1e-4, but achieve the convergence criteria for the
-# rest of the model run time. For a production model, one would choose a tolerance
-# of 1e-5, or even stricter if possible.
+# rest of the model run time. For a production model, a tolerance of at least 1e-5
+# should be selected. Reaching these tolerances may require on the order of
+# 50-100 nonlinear iterations, and the maximum number of nonlinear iterations
+# will need to be increased accordingly.
 set Max nonlinear iterations                   = 25
 set Nonlinear solver tolerance                 = 1e-4
 set Dimension                                  = 2
@@ -38,7 +40,7 @@ set Output directory                           = output_coupled_two_phase_tian_p
 
 # Checkpoint the model so that a restart can be done if needed.
 subsection Checkpointing
-  set Steps between checkpoint = 200
+  set Steps between checkpoint = 20
 end
 
 subsection Discretization
@@ -229,6 +231,9 @@ subsection Material model
   subsection Visco Plastic
     # Use the composite flow law to simulate dislocation and diffusion creep
     # of dry olivine with the parameters from Hirth & Kohlstaedt 2004 (10.1029/138GM06).
+    # This will be inconsistent with the physical state of the model (the solid is hydrated),
+    # but ongoing work will introduce a more self-consistent coupling between the hydration
+    # state and the rheology.
     set Viscosity averaging scheme                      = harmonic
     set Viscous flow law                                = composite
 

cookbooks/tian_parameterization_kinematic_slab/doc/tian_parameterization_kinematic_slab.md

@@ -22,8 +22,8 @@ schematic diagram showing where the model is in the context of a generic subduct
 :::{note}
 The low fluid viscosity and high contrast in density between the solid and fluid phases leads to large fluid
 velocities on the order of m/yr. Combined with the need for high resolution around the fluid to resolve compaction
-viscosity gradients means that this cookbook requires more computational power than some other cookbooks within
-ASPECT.
+viscosity gradients means that this cookbook requires more computational resources than some other cookbooks within
+ASPECT. In detail, this simulation was run on 128 processors for approximately 2 hours.
 :::
 
 ```{figure-md} fig:schematic-diagram-overview
@@ -40,19 +40,20 @@ Schematic diagram showing the model design.
 
 ## The input file
 This cookbook includes two input files, one which advects the porosity (free fluid) compositional field through the
-darcy velocity, and one which advects the porosity compositional field according to the coupled two-phase reactions
-solved for in the Reactive Fluid Transport model. Both models feature dehydration and rehydration reactions using
-the Tian parameterization and both input files can be found at [cookbooks/tian_parameterization_kinematic_slab](https://www.github.com/geodynamics/aspect/blob/main/cookbooks/tian_parameterization_kinematic_slab/). One
+darcy velocity, and one which advects the porosity compositional field according to the fluid velocities obtained by solving the fully coupled two-phase flow system of equations. Both models feature dehydration and rehydration reactions using the Tian parameterization and both input files can be found at [cookbooks/tian_parameterization_kinematic_slab](https://www.github.com/geodynamics/aspect/blob/main/cookbooks/tian_parameterization_kinematic_slab/). One
 important problem in models that track the partitioning of fluid into/out of a solid phase is that these
 interactions can be much faster than the time step of the model. To model these type of processes, ASPECT uses
 operator splitting (see also {ref}`sec:benchmark:operator-splitting`): Reactions are solved on a different time
 scale than advection. For this model, this means that at the beginning of each time step, all hydration/dehydration
 reactions are solved using several shorter sub-time steps. In the input file, we have to choose both the size of
-these sub-time steps and the rate (or characteristic time scale) of the solid-fluid reactions, and they have to be
+these sub-time steps and the rate (or characteristic time scale) of the solid-fluid reactions, which must be
 consistent in the sense that the operator splitting time step can not be larger than the reaction time scale. We
 choose a conservative value for the fluid-solid reaction time scale of 50 kyr for model stability and computational
 efficiency. For a production model one would likely want to reduce this value to be as low as possible without
-affecting model convergence.
+affecting model convergence, since in reality the timescale for the fluid-solid reactions is likely much lower than
+the advection time step. For example, at 300$^o$C it has been estimated that to fully serpentinize a 1km line of
+mantle would take ~ 10 kyr, and this timescale seems to decrease with increasing temperature {cite:t}`macdonald_fyfe_1985`. However, smaller time scales would lead to very large changes in porosity that make the
+nonlinear solver converge less reliably unless much smaller advection time steps/higher resolution are enforced.
 
 ## Model Evolution
 As the slab 'subducts' across the base of the model, the hot ambient mantle heats the slab leading to progressive
@@ -109,3 +110,8 @@ can partition into the peridotite mantle.
 this dip will impact the magnitude of the slab-orthogonal component of gravity which has a large control on the
 fluid pathways. Varying the convergence rate will influence how much heat is able to conduct from the mantle into
 the slab and affect the timing of dehydration.
+
+It is worth reiterating that to extend this model towards a production simulation, the nonlinear solver tolerance
+should be stricter (at least 1e-5), and reaction rates should be reduced. This cookbook serves as a base model
+that showcases the use of tian approximation in the reactive two-phase fluid material model, and these
+simplifications were made in this cookbook to make the model more user-friendly.

cookbooks/tian_parameterization_kinematic_slab/uncoupled-two-phase-tian-parameterization-kinematic-slab.prm

@@ -10,8 +10,9 @@ include $ASPECT_SOURCE_DIR/cookbooks/tian_parameterization_kinematic_slab/couple
 set Output directory = output_uncoupled_two_phase_tian_parameterization_kinematic_slab
 
   subsection Stokes solver parameters
-    # Lower solver tolerance to improve efficiency. Since the problem is simpler (i.e
-    # no compaction), we can get away with this.
+    # As the system of equations being solved is simpler, the linear solver tolerance can be reduced
+    # to improve efficiency. Note that in production models this value should be tested and decreased
+    # until it no longer has an effect on the model results when using the "iterated Advection and Stokes"
     set Linear solver tolerance             = 1e-8
     set Number of cheap Stokes solver steps = 1000
   end

doc/modules/changes/20240605_douglas

@@ -0,0 +1,6 @@
+Added: Cookbook which showcases the use of the tian approximation
+in the Reactive Fluid Transport material model. The cookbook features
+a kinematically subducting slab with an initial hydration state
+dehydrating as it advects through a hot mantle wedge.
+<br>
+(Daniel Douglas, 2024/06/05)

doc/sphinx/references.bib

@@ -12214,3 +12214,14 @@ @article{tian_et_al_2019
 doi = {https://doi.org/10.1029/2019GC008488},
 year = {2019}
 }
+
+@article{macdonald_fyfe_1985,
+author = {MacDonald, AH and Fyfe, WS},
+title = {Rate of Serpentinization in Seafloor Experiments},
+journal = {Tectonophysics},
+volume = {116},
+number = {1-2},
+pages = {123--135},
+doi = {10.1016/0040-1951(85)90225-2},
+year = {1985}
+}

tests/coupled_two_phase_tian_parameterization_kinematic_slab/screen-output

@@ -70,4 +70,3 @@ Termination requested by criterion: end time
 
 
 
-

tests/uncoupled_two_phase_tian_parameterization_kinematic_slab/screen-output

@@ -66,4 +66,3 @@ Termination requested by criterion: end time
 
 
 
-