Add Viscosity Limit to Reactive Fluid Transport

doc/modules/changes/20260305_danieldouglas92

@@ -0,0 +1,4 @@
+Added: A parameter to the reactive fluid transport model that limits
+the amount of fluid weakening on the solid viscosity.
+<br>
+(Daniel Douglas, 2026/03/05)

doc/sphinx/parameters/Material_20model.md

@@ -200,7 +200,7 @@ Viscous stress may also be limited by a non-linear stress limiter that has a for
 
  The value for the components of this formula and additional parameters are read from the parameter file in subsection  ’Material model/Visco Plastic’.
 
-‘viscoelastic’: An implementation of a simple linear viscoelastic rheology that only includes the deviatoric components of elasticity. Specifically, the viscoelastic rheology only takes into account the elastic shear strength (e.g., shear modulus), while the tensile and volumetric strength (e.g., Young’s and bulk modulus) are not considered. The model is incompressible and allows specifying an arbitrary number of compositional fields, where each field represents a different rock type or component of the viscoelastic stress tensor. The stress tensor in 2d and 3d, respectively, contains 3 or 6 components. The compositional fields representing these components must be named and listed in a very specific format, which is designed to minimize mislabeling stress tensor components as distinct ’compositional rock types’ (or vice versa). For 2d models, the first six compositional fields must be labeled ’stress\_xx’, ’stress\_yy’ and ’stress\_xy’, ’stress\_xx\_old’, ’stress\_yy\_old’ and ’stress\_xy\_old’, In 3d, the first twelve compositional fields must be labeled ’stress\_xx’, ’stress\_yy’, ’stress\_zz’, ’stress\_xy’, ’stress\_xz’, ’stress\_yz’, ’stress\_xx\_old’, ’stress\_yy\_old’, ’stress\_zz\_old’, ’stress\_xy\_old’, ’stress\_xz\_old’, ’stress\_yz\_old’.
+‘viscoelastic’: An implementation of a simple linear viscoelastic rheology that only includes the deviatoric components of elasticity. Specifically, the viscoelastic rheology only takes into account the elastic shear strength (e.g., shear modulus), while the tensile and volumetric strength (e.g., Young’s and bulk modulus) are not considered. The model is incompressible and allows specifying an arbitrary number of compositional fields, where each field represents a different rock type or component of the viscoelastic stress tensor. The stress tensor in 2d and 3d, respectively, contains 3 or 6 components. The compositional fields representing these components must be named and listed in a very specific format, which is designed to minimize mislabeling stress tensor components as distinct ’compositional rock types’ (or vice versa). For 2d models, the first six compositional fields of type stress must be labeled ’ve\_stress\_xx’, ’ve\_stress\_yy’ and ’ve\_stress\_xy’, ’ve\_stress\_xx\_old’, ’ve\_stress\_yy\_old’ and ’ve\_stress\_xy\_old’, In 3d, the first twelve compositional fields of type stress must be labeled ’ve\_stress\_xx’, ’ve\_stress\_yy’, ’ve\_stress\_zz’, ’ve\_stress\_xy’, ’ve\_stress\_xz’, ’ve\_stress\_yz’, ’ve\_stress\_xx\_old’, ’ve\_stress\_yy\_old’, ’ve\_stress\_zz\_old’,  ’ve\_stress\_xy\_old’, ’ve\_stress\_xz\_old’, ’ve\_stress\_yz\_old’.
 
  Expanding the model to include non-linear viscous flow (e.g., diffusion/dislocation creep) and plasticity would produce a constitutive relationship commonly referred to as partial elastoviscoplastic (e.g., pEVP) in the geodynamics community. While extensively discussed and applied within the geodynamics literature, notable references include: Moresi et al. (2003), J. Comp. Phys., v. 184, p. 476-497. Gerya and Yuen (2007), Phys. Earth. Planet. Inter., v. 163, p. 83-105. Gerya (2010), Introduction to Numerical Geodynamic Modeling. Kaus (2010), Tectonophysics, v. 484, p. 36-47. Choi et al. (2013), J. Geophys. Res., v. 118, p. 2429-2444. Keller et al. (2013), Geophys. J. Int., v. 195, p. 1406-1442.
 
@@ -3821,6 +3821,15 @@ Also note that the fluid reaction time scale has to be larger than or equal to t
 **Documentation:** Lower cutoff for the compaction viscosity. Units: $\text{Pa}\text{s}$.
 ::::
 
+::::{dropdown} __Parameter:__ {ref}`Minimum weakened viscosity<parameters:Material_20model/Reactive_20Fluid_20Transport_20Model/Minimum_20weakened_20viscosity>`
+:name: parameters:Material_20model/Reactive_20Fluid_20Transport_20Model/Minimum_20weakened_20viscosity
+**Default value:** 1e17
+
+**Pattern:** [Double 0...MAX_DOUBLE (inclusive)]
+
+**Documentation:** Lower cutoff for the viscosity of the solid phase in the presence of fluids. Units: $\text{Pa}\text{s}$.
+::::
+
 ::::{dropdown} __Parameter:__ {ref}`Reference fluid density<parameters:Material_20model/Reactive_20Fluid_20Transport_20Model/Reference_20fluid_20density>`
 :name: parameters:Material_20model/Reactive_20Fluid_20Transport_20Model/Reference_20fluid_20density
 **Default value:** 2500

include/aspect/material_model/reactive_fluid_transport.h

@@ -164,6 +164,7 @@ namespace aspect
         double shear_to_bulk_viscosity_ratio;
         double min_compaction_viscosity;
         double max_compaction_viscosity;
+        double min_weakened_viscosity;
         double reference_permeability;
         double alpha_phi;
         double reference_T;

source/material_model/reaction_model/katz2003_mantle_melting.cc

@@ -344,15 +344,6 @@ namespace aspect
                 melt_out->compaction_viscosities[i] *= visc_temperature_dependence;
               }
           }
-
-        if (this->include_melt_transport() && in.requests_property(MaterialProperties::viscosity))
-          {
-            for (unsigned int i=0; i<in.n_evaluation_points(); ++i)
-              {
-                const double porosity = std::min(1.0, std::max(in.composition[i][porosity_idx],0.0));
-                out.viscosities[i] *= std::exp(- alpha_phi * porosity);
-              }
-          }
       }
 
 

source/material_model/reactive_fluid_transport.cc

@@ -180,21 +180,22 @@ namespace aspect
     {
       base_model->evaluate(in,out);
 
-      if (fluid_solid_reaction_scheme != katz2003)
-        {
-          const unsigned int porosity_idx = this->introspection().compositional_index_for_name("porosity");
+      const unsigned int porosity_idx = this->introspection().compositional_index_for_name("porosity");
 
-          // Modify the viscosity from the base model based on the presence of fluid.
-          if (in.requests_property(MaterialProperties::viscosity))
+      // Modify the viscosity from the base model based on the presence of fluid.
+      if (in.requests_property(MaterialProperties::viscosity))
+        {
+          // Scale the base model viscosity value based on the porosity.
+          for (unsigned int q=0; q<out.n_evaluation_points(); ++q)
             {
-              // Scale the base model viscosity value based on the porosity.
-              for (unsigned int q=0; q<out.n_evaluation_points(); ++q)
-                {
-                  const double porosity = std::max(in.composition[q][porosity_idx],0.0);
-                  out.viscosities[q] *= (1.0 - porosity) * std::exp(- alpha_phi * porosity);
-                }
+              const double porosity = std::max(in.composition[q][porosity_idx],0.0);
+              const double weakened_viscosity = out.viscosities[q] * (1.0 - porosity) * std::exp(- alpha_phi * porosity);
+              out.viscosities[q] = std::max(weakened_viscosity, min_weakened_viscosity);
             }
+        }
 
+      if (fluid_solid_reaction_scheme != katz2003)
+        {
           // Fill the melt outputs if they exist. Note that the MeltOutputs class was originally
           // designed for two-phase flow material models in ASPECT that model the flow of melt,
           // but can be reused for a geofluid of arbitrary composition.
@@ -377,6 +378,10 @@ namespace aspect
                              boost::lexical_cast<std::string>(std::numeric_limits<double>::max()),
                              Patterns::Double (0),
                              "Upper cutoff for the compaction viscosity. Units: $\\text{Pa}\\text{s}$.");
+          prm.declare_entry ("Minimum weakened viscosity", "1e17",
+                             Patterns::Double (0),
+                             "Lower cutoff for the viscosity of the solid phase in the presence of fluids. "
+                             "Units: $\\text{Pa}\\text{s}$.");
           prm.declare_entry ("Reference fluid viscosity", "10",
                              Patterns::Double (0),
                              "The value of the constant melt/fluid viscosity $\\eta_f$. Units: $\\text{Pa}\\text{s}$.");
@@ -448,6 +453,7 @@ namespace aspect
           shear_to_bulk_viscosity_ratio     = prm.get_double ("Shear to bulk viscosity ratio");
           max_compaction_viscosity          = prm.get_double ("Maximum compaction viscosity");
           min_compaction_viscosity          = prm.get_double ("Minimum compaction viscosity");
+          min_weakened_viscosity            = prm.get_double ("Minimum weakened viscosity");
           eta_f                             = prm.get_double ("Reference fluid viscosity");
           reference_permeability            = prm.get_double ("Reference permeability");
           alpha_phi                         = prm.get_double ("Exponential fluid weakening factor");

tests/reactive_fluid_transport_min_viscosity.prm

@@ -0,0 +1,32 @@
+# This model tests that the minimum viscosity within the reactive fluid transport model
+# correctly stops the viscosity from decreasing below a threshold in the presence
+# of fluids. The model prescribes a porosity of 10 % and an exponential fluid 
+# weakening factor of 30. Without a minimum value for the viscosity, this would 
+# result in the material viscosity of 4.481e17 Pa s, which is below the minimum
+# value that is specified in the Visco Plastic base material model. By setting
+# the parameter Minimum weakened viscosity to 1e19 Pa s, we ensure that the viscosity 
+# of the material does not drop below this value, even in the presence of fluids.
+include $ASPECT_SOURCE_DIR/tests/reactive_fluid_transport_zero_solubility.prm
+
+set Output directory = reactive_fluid_transport_min_viscosity
+set End time         = 0
+
+# 10 % porosity initially.
+subsection Initial composition model
+  set Model name = function
+
+  subsection Function
+    set Function expression = 0.1; 0.0
+  end
+end
+
+subsection Material model
+  subsection Reactive Fluid Transport Model
+    set Exponential fluid weakening factor = 30
+    set Minimum weakened viscosity         = 1e19
+  end
+end
+
+subsection Postprocess
+  set List of postprocessors = material statistics
+end

tests/reactive_fluid_transport_min_viscosity/screen-output

@@ -0,0 +1,22 @@
+
+Number of active cells: 1 (on 1 levels)
+Number of degrees of freedom: 74 (18+7+18+4+9+9+9)
+
+*** Timestep 0:  t=0 years, dt=0 years
+   Solving temperature system... 0 iterations.
+   Solving porosity system ... 0 iterations.
+   Skipping bound_fluid composition solve because RHS is zero.
+   Rebuilding Stokes preconditioner...
+   Solving Stokes system (AMG)... 0+0 iterations.
+   Solving fluid velocity system... 1 iterations.
+      Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.3813e-16, 4.89778e-17, 0, 1.55195e-09
+      Relative nonlinear residual (total system) after nonlinear iteration 1: 1.55195e-09
+
+
+   Postprocessing:
+     Average density / Average viscosity / Total mass:  3000 kg/m^3, 1e+19 Pa s, 1.2e+10 kg
+
+Termination requested by criterion: end time
+
+
+

tests/reactive_fluid_transport_min_viscosity/statistics

@@ -0,0 +1,18 @@
+# 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 Stokes solver
+# 13: Velocity iterations in Stokes preconditioner
+# 14: Schur complement iterations in Stokes preconditioner
+# 15: Average density (kg/m^3)
+# 16: Average viscosity (Pa s)
+# 17: Total mass (kg)
+0 0.000000000000e+00 0.000000000000e+00 1 22 9 18 1 0 0 0 0 0 0 3.00000000e+03 1.00000000e+19 1.20000000e+10