Add cutoff for water fugacity calculation

include/aspect/material_model/rheology/compositional_viscosity_prefactors.h

@@ -103,6 +103,7 @@ namespace aspect
           // Initialize variables for the water fugacity calculation, from HK04
           std::vector<double> diffusion_water_fugacity_exponents;
           std::vector<double> dislocation_water_fugacity_exponents;
+          std::vector<double> weight_percent_water_for_dry_creep_cutoff;
 
           // From Hirth & Kohlstaedt 2004, equation 6
           const double A_H2O = 2.6e-5; // 1/Pa

source/material_model/rheology/compositional_viscosity_prefactors.cc

@@ -61,7 +61,7 @@ namespace aspect
               // olivine assuming a composition of 90 mol% Forsterite and 10 mol% Fayalite from Hirth
               // and Kohlstaedt 2004 10.1029/138GM06.
               const unsigned int bound_fluid_idx = this->introspection().compositional_index_for_name("bound_fluid");
-              const double weight_fraction_H2O = in.composition[q][bound_fluid_idx]; // mass fraction of bound water
+              const double weight_fraction_H2O = std::max(weight_percent_water_for_dry_creep_cutoff[composition_index], in.composition[q][bound_fluid_idx]); // mass fraction of bound water
               const double weight_fraction_olivine = 1 - weight_fraction_H2O; // mass fraction of olivine
               const double COH = (weight_fraction_H2O/molar_mass_H2O) / (weight_fraction_olivine/molar_mass_olivine) * 1e6; // COH in H / Si ppm
               const double point_water_fugacity = COH / A_H2O *
@@ -84,6 +84,14 @@ namespace aspect
       void
       CompositionalViscosityPrefactors<dim>::declare_parameters (ParameterHandler &prm)
       {
+        prm.declare_entry ("Minimum weight percent bound water content for fugacity", "6.15e-6",
+                           Patterns::List(Patterns::Double(0)),
+                           "The minimum water content for the HK04 olivine hydration viscosity "
+                           "prefactor scheme. This acts as the cutoff between 'dry' creep and 'wet' creep "
+                           "for olivine, and the default value is chosen based on the value reported by "
+                           "Hirth & Kohlstaedt 2004. For weight percent bound water beneath this value, "
+                           "this value is used instead to compute the water fugacity. Units: weight %.");
+
         prm.declare_entry ("Water fugacity exponents for diffusion creep", "0.0",
                            Patterns::List(Patterns::Double(0)),
                            "List of water fugacity exponents for diffusion creep for "
@@ -144,6 +152,8 @@ namespace aspect
                                                  options);
             dislocation_water_fugacity_exponents = Utilities::MapParsing::parse_map_to_double_array (prm.get("Water fugacity exponents for dislocation creep"),
                                                    options);
+            weight_percent_water_for_dry_creep_cutoff = Utilities::MapParsing::parse_map_to_double_array (prm.get("Minimum weight percent bound water content for fugacity"),
+                                                        options);
           }
       }
     }

tests/hk04_olivine_composite_hydration_prefactor_cutoff.prm

@@ -0,0 +1,19 @@
+# A test case that checks the water fugacity viscous prefactor
+# multiplication scheme in combination with the visco plastic material model.
+# This test is almost identical to hk04_olivine_composite_hydration_prefactor.prm,
+# with the only difference being that the bound_fluid is 0. This test checks that
+# the HK04 olivine hydration Viscosity prefactor scheme will not divide by 0, but
+# will instead use the user input value set using 'Minimum weight percent bound
+# water content for fugacity'. The default value is 6.15e-6, and so following the
+# same calculation outlined in hk04_olivine_composite_hydration_prefactor.prm, the
+# material model should return 3.1724782150033374e22 Pa s for the viscosity.
+
+include $ASPECT_SOURCE_DIR/tests/hk04_olivine_composite_hydration_prefactor.prm
+
+# Set the bound water everywhere in the model to be 0.
+subsection Initial composition model
+  set Model name = function
+  subsection Function
+    set Function expression = 0.0
+  end
+end

tests/hk04_olivine_composite_hydration_prefactor_cutoff/screen-output

@@ -0,0 +1,30 @@
+-----------------------------------------------------------------------------
+-----------------------------------------------------------------------------
+-----------------------------------------------------------------------------
+
+-----------------------------------------------------------------------------
+-----------------------------------------------------------------------------
+Number of active cells: 100 (on 1 levels)
+Number of degrees of freedom: 1,885 (882+121+441+441)
+
+*** Timestep 0:  t=0 years, dt=0 years
+   Solving temperature system... 0 iterations.
+   Skipping bound_fluid composition solve because RHS is zero.
+   Rebuilding Stokes preconditioner...
+   Solving Stokes system... 14+0 iterations.
+      Relative nonlinear residual (Stokes system) after nonlinear iteration 1: 1
+
+
+   Postprocessing:
+     Average density / Average viscosity / Total mass:  3198 kg/m^3, 3.172e+22 Pa s, 3.198e+13 kg
+     Writing graphical output:                          output-hk04_olivine_composite_hydration_prefactor_cutoff/solution/solution-00000
+
+Termination requested by criterion: end time
+
+
++----------------------------------------------+------------+------------+
++----------------------------------+-----------+------------+------------+
++----------------------------------+-----------+------------+------------+
+
+-----------------------------------------------------------------------------
+-----------------------------------------------------------------------------