Add Tian Parametrization to Reactive Fluid Transport

doc/modules/changes/20240224_danieldouglas92

@@ -0,0 +1,5 @@
+Implemented a new fluid-rock interaction scheme for the reactive fluid
+transport material model based on a parametrization by Tian et al., 2019 G3,
+https://doi.org/10.1029/2019GC008488.
+<br>
+ (Daniel Douglas, 2024/02/24)

include/aspect/material_model/reactive_fluid_transport.h

@@ -62,6 +62,18 @@ namespace aspect
          */
         virtual double reference_darcy_coefficient () const override;
 
+
+        /**
+         * Compute the maximum allowed bound water content at the input
+         * pressure and temperature conditions. This is used to determine
+         * how free water interacts with the solid phase.
+         * @param in Object that contains the current conditions.
+         * @param q unsigned int from 0-3 indexing which rock phase the equilbrium
+         * bound water content is being calculated for
+         */
+        std::vector<double> tian_equilibrium_bound_water_content(const MaterialModel::MaterialModelInputs<dim> &in,
+                                                                 unsigned int q) const;
+
         /**
          * Compute the free fluid fraction that can be present in the material based on the
          * fluid content of the material and the fluid solubility for the given input conditions.
@@ -138,6 +150,46 @@ namespace aspect
         // Time scale for fluid release and absorption.
         double fluid_reaction_time_scale;
 
+        // The maximum water content for each of the 4 rock types in the tian approximation
+        // method. These are important for keeping the polynomial bounded within reasonable
+        // values.
+        double tian_max_peridotite_water;
+        double tian_max_gabbro_water;
+        double tian_max_MORB_water;
+        double tian_max_sediment_water;
+
+        // The following coefficients are taken from a publication from Tian et al., 2019, and can be found
+        // in Table 3 (Gabbro), Table B1 (MORB), Table B2 (Sediments) and Table B3 (peridotite).
+        // LR refers to the effective enthalpy change for devolatilization reactions,
+        // csat is the saturated mass fraction of water in the solid, and Td is the
+        // onset temperature of devolatilization for water.
+        std::vector<double> LR_peridotite_poly_coeffs {-19.0609, 168.983, -630.032, 1281.84, -1543.14, 1111.88, -459.142, 95.4143, 1.97246};
+        std::vector<double> csat_peridotite_poly_coeffs {0.00115628, 2.42179};
+        std::vector<double> Td_peridotite_poly_coeffs {-15.4627, 94.9716, 636.603};
+
+        std::vector<double> LR_gabbro_poly_coeffs {-1.81745, 7.67198, -10.8507, 5.09329, 8.14519};
+        std::vector<double> csat_gabbro_poly_coeffs {-0.0176673, 0.0893044, 1.52732};
+        std::vector<double> Td_gabbro_poly_coeffs {-1.72277, 20.5898, 637.517};
+
+        std::vector<double> LR_MORB_poly_coeffs {-1.78177, 7.50871, -10.4840, 5.19725, 7.96365};
+        std::vector<double> csat_MORB_poly_coeffs {0.0102725, -0.115390, 0.324452, 1.41588};
+        std::vector<double> Td_MORB_poly_coeffs {-3.81280, 22.7809, 638.049};
+
+        std::vector<double> LR_sediment_poly_coeffs {-2.03283, 10.8186, -21.2119, 18.3351, -6.48711, 8.32459};
+        std::vector<double> csat_sediment_poly_coeffs {-0.150662, 0.301807, 1.01867};
+        std::vector<double> Td_sediment_poly_coeffs {2.83277, -24.7593, 85.9090, 524.898};
+
+        std::vector<std::vector<double>> devolatilization_enthalpy_changes {LR_peridotite_poly_coeffs, LR_gabbro_poly_coeffs, \
+                                                                             LR_MORB_poly_coeffs, LR_sediment_poly_coeffs
+                                                                            };
+
+        std::vector<std::vector<double>> water_mass_fractions {csat_peridotite_poly_coeffs, csat_gabbro_poly_coeffs, \
+                                                                csat_MORB_poly_coeffs, csat_sediment_poly_coeffs
+                                                               };
+
+        std::vector<std::vector<double>> devolatilization_onset_temperatures {Td_peridotite_poly_coeffs, Td_gabbro_poly_coeffs, \
+                                                                               Td_MORB_poly_coeffs, Td_sediment_poly_coeffs
+                                                                              };
         /**
          * Enumeration for selecting which type of scheme to use for
          * reactions between fluids and solids. The available
@@ -164,12 +216,26 @@ namespace aspect
          * base model to be incompressible, as otherwise the advection
          * equation would only be valid for mass and not volume
          * fractions.
+         *
+         * The tian approximation model implements parametrized phase
+         * diagrams from Tian et al., 2019 G3, https://doi.org/10.1029/2019GC008488
+         * and calculates the fluid-solid reactions for four different rock types:
+         * sediments, MORB, gabbro and peridotite. This is achieved by calculating the
+         * maximum allowed bound water content for each composition at the current
+         * Pressure-Temperature conditions, and releasing bound water as free water if:
+         * (maximum bound water content < current bound water content)
+         * or incorporating free water (if present) into the solid phase as bound water:
+         * maximum bound water content > current bound water content
+         * This model requires that 4 compositional fields named after the 4 different rock
+         * types exist in the input file.
          */
         enum ReactionScheme
         {
           no_reaction,
-          zero_solubility
-        } fluid_solid_reaction_scheme;
+          zero_solubility,
+          tian_approximation
+        }
+        fluid_solid_reaction_scheme;
     };
   }
 }

source/material_model/reactive_fluid_transport.cc

@@ -50,6 +50,63 @@ namespace aspect
 
 
 
+    template <int dim>
+    std::vector<double>
+    ReactiveFluidTransport<dim>::
+    tian_equilibrium_bound_water_content (const MaterialModel::MaterialModelInputs<dim> &in,
+                                          unsigned int q) const
+    {
+      // Pressure, which must be in GPa for the parametrization, or GPa^-1
+      const double pressure = in.pressure[q]<=0 ? 1e-12 : in.pressure[q]/1.e9;
+      const double inverse_pressure = std::pow(pressure, -1);
+
+      // Create arrays that will store the values of the polynomials at the current pressure
+      std::vector<double> LR_values(4);
+      std::vector<double> csat_values(4);
+      std::vector<double> Td_values(4);
+
+      // Loop over the four rock types (peridotite, gabbro, MORB, sediment) and the polynomial
+      // coefficients to fill the vectors defined above. The polynomials for LR are defined in
+      // equations 13, B2, B10, and B18. csat polynomials are defined in equations 14, B1, B9, and B17.
+      // Td polynomials are defined in equations 15, B3, B11, and B19.
+      for (unsigned int i = 0; i<devolatilization_enthalpy_changes.size(); ++i)
+        {
+          for (unsigned int j = 0; j<devolatilization_enthalpy_changes[i].size(); ++j)
+            {
+              LR_values[i] += devolatilization_enthalpy_changes[i][j] * std::pow(inverse_pressure, devolatilization_enthalpy_changes[i].size() - 1 - j);
+            }
+
+          for (unsigned int j = 0; j<water_mass_fractions[i].size(); ++j)
+            {
+              csat_values[i] += i==3 ? water_mass_fractions[i][j] * std::pow(std::log10(pressure), water_mass_fractions[i].size() - 1 - j) :\
+                                water_mass_fractions[i][j] * std::pow(pressure, water_mass_fractions[i].size() - 1 - j);
+            }
+
+          for (unsigned int j = 0; j<devolatilization_onset_temperatures[i].size(); ++j)
+            {
+              Td_values[i] += devolatilization_onset_temperatures[i][j] * std::pow(pressure, devolatilization_onset_temperatures[i].size() - 1 - j);
+            }
+        }
+
+      // Create an array for the equilibrium bound water content that is calculated from these polynomials
+      std::vector<double> eq_bound_water_content(4);
+
+      // Define the maximum bound water content allowed for the four different rock compositions
+      std::vector<double> max_bound_water_content = {tian_max_peridotite_water, tian_max_gabbro_water, tian_max_MORB_water, tian_max_sediment_water};
+
+      // Loop over all rock compositions and fill the equilibrium bound water content, divide by 100 to convert
+      // from percentage to fraction (equation 1)
+      for (unsigned int k = 0; k<LR_values.size(); ++k)
+        {
+          eq_bound_water_content[k] = (std::min(std::exp(csat_values[k]) * \
+                                                std::exp(std::exp(LR_values[k]) * (1/in.temperature[q] - 1/Td_values[k])), \
+                                                max_bound_water_content[k]) / 100.0);
+        }
+      return eq_bound_water_content;
+    }
+
+
+
     template <int dim>
     void
     ReactiveFluidTransport<dim>::
@@ -59,7 +116,7 @@ namespace aspect
       for (unsigned int q=0; q<in.temperature.size(); ++q)
         {
           const unsigned int porosity_idx = this->introspection().compositional_index_for_name("porosity");
-
+          const unsigned int bound_fluid_idx = this->introspection().compositional_index_for_name("bound_fluid");
           switch (fluid_solid_reaction_scheme)
             {
               case no_reaction:
@@ -75,18 +132,45 @@ namespace aspect
                 // The fluid volume fraction in equilibrium with the solid
                 // at any point (stored in the melt_fractions vector) is
                 // equal to the sum of the bound fluid content and porosity.
-                const unsigned int bound_fluid_idx = this->introspection().compositional_index_for_name("bound_fluid");
                 melt_fractions[q] = in.composition[q][bound_fluid_idx] + in.composition[q][porosity_idx];
                 break;
               }
+              case tian_approximation:
+              {
+                // The bound fluid content is calculated using parametrized phase
+                // diagrams for four different rock types: sediment, MORB, gabbro, and
+                // peridotite.
+                const unsigned int sediment_idx = this->introspection().compositional_index_for_name("sediment");
+                const unsigned int MORB_idx = this->introspection().compositional_index_for_name("MORB");
+                const unsigned int gabbro_idx = this->introspection().compositional_index_for_name("gabbro");
+                const unsigned int peridotite_idx = this->introspection().compositional_index_for_name("peridotite");
+
+                // Initialize a vector that stores the compositions (mass fractions) for
+                // the four different rock compositions,
+                std::vector<double> tracked_rock_mass_fractions(4);
+                tracked_rock_mass_fractions[0] = (in.composition[q][peridotite_idx]);
+                tracked_rock_mass_fractions[1] = (in.composition[q][gabbro_idx]);
+                tracked_rock_mass_fractions[2] = (in.composition[q][MORB_idx]);
+                tracked_rock_mass_fractions[3] = (in.composition[q][sediment_idx]);
+
+                // The bound water content (water within the solid phase) for the four different rock types
+                std::vector<double> tian_eq_bound_water_content = tian_equilibrium_bound_water_content(in, q);
+
+                // average the water content between the four different rock types
+                double average_eq_bound_water_content = MaterialUtilities::average_value (tracked_rock_mass_fractions, tian_eq_bound_water_content, MaterialUtilities::arithmetic);
+
+                // The fluid volume fraction in equilibrium with the solid (stored in the melt_fractions vector)
+                // is equal to the sum of the porosity and the change in bound fluid content
+                // (current bound fluid - updated average bound fluid).
+                melt_fractions[q] = std::max(in.composition[q][bound_fluid_idx] + in.composition[q][porosity_idx] - average_eq_bound_water_content, 0.0);
+                break;
+              }
               default:
               {
                 AssertThrow(false, ExcNotImplemented());
                 break;
               }
             }
-
-
         }
     }
 
@@ -246,12 +330,26 @@ namespace aspect
                              "computed. If the model does not use operator splitting, this parameter is not used. "
                              "Units: yr or s, depending on the ``Use years "
                              "in output instead of seconds'' parameter.");
+          prm.declare_entry ("Maximum weight percent water in sediment", "3",
+                             Patterns::Double (0),
+                             "The maximum allowed weight percent that the sediment composition can hold.");
+          prm.declare_entry ("Maximum weight percent water in MORB", "2",
+                             Patterns::Double (0),
+                             "The maximum allowed weight percent that the sediment composition can hold.");
+          prm.declare_entry ("Maximum weight percent water in gabbro", "1",
+                             Patterns::Double (0),
+                             "The maximum allowed weight percent that the sediment composition can hold.");
+          prm.declare_entry ("Maximum weight percent water in peridotite", "8",
+                             Patterns::Double (0),
+                             "The maximum allowed weight percent that the sediment composition can hold.");
           prm.declare_entry ("Fluid-solid reaction scheme", "no reaction",
-                             Patterns::Selection("no reaction|zero solubility"),
+                             Patterns::Selection("no reaction|zero solubility|tian approximation"),
                              "Select what type of scheme to use for reactions between fluid and solid phases. "
                              "The current available options are models where no reactions occur between "
                              "the two phases, or the solid phase is insoluble (zero solubility) and all "
-                             "of the bound fluid is released into the fluid phase.");
+                             "of the bound fluid is released into the fluid phase, tian approximation "
+                             "use polynomials to describe hydration and dehydration reactions for four different "
+                             "rock compositions as defined in Tian et al., 2019.");
         }
         prm.leave_subsection();
       }
@@ -281,6 +379,11 @@ namespace aspect
           fluid_compressibility             = prm.get_double ("Fluid compressibility");
           fluid_reaction_time_scale         = prm.get_double ("Fluid reaction time scale for operator splitting");
 
+          tian_max_peridotite_water         = prm.get_double ("Maximum weight percent water in peridotite");
+          tian_max_gabbro_water             = prm.get_double ("Maximum weight percent water in gabbro");
+          tian_max_MORB_water               = prm.get_double ("Maximum weight percent water in MORB");
+          tian_max_sediment_water           = prm.get_double ("Maximum weight percent water in sediment");
+
           // Create the base model and initialize its SimulatorAccess base
           // class; it will get a chance to read its parameters below after we
           // leave the current section.
@@ -293,9 +396,29 @@ namespace aspect
 
           // Reaction scheme parameter
           if (prm.get ("Fluid-solid reaction scheme") == "zero solubility")
-            fluid_solid_reaction_scheme = zero_solubility;
+            {
+              fluid_solid_reaction_scheme = zero_solubility;
+            }
           else if (prm.get ("Fluid-solid reaction scheme") == "no reaction")
-            fluid_solid_reaction_scheme = no_reaction;
+            {
+              fluid_solid_reaction_scheme = no_reaction;
+            }
+          else if (prm.get ("Fluid-solid reaction scheme") == "tian approximation")
+            {
+              AssertThrow(this->introspection().compositional_name_exists("sediment"),
+                          ExcMessage("The Tian approximation only works "
+                                     "if there is a compositional field called sediment."));
+              AssertThrow(this->introspection().compositional_name_exists("MORB"),
+                          ExcMessage("The Tian approximation only works "
+                                     "if there is a compositional field called MORB."));
+              AssertThrow(this->introspection().compositional_name_exists("gabbro"),
+                          ExcMessage("The Tian approximation only works "
+                                     "if there is a compositional field called gabbro."));
+              AssertThrow(this->introspection().compositional_name_exists("peridotite"),
+                          ExcMessage("The Tian approximation only works "
+                                     "if there is a compositional field called peridotite."));
+              fluid_solid_reaction_scheme = tian_approximation;
+            }
           else
             AssertThrow(false, ExcMessage("Not a valid fluid-solid reaction scheme"));
 
@@ -311,6 +434,12 @@ namespace aspect
                           ExcMessage("The Fluid-reaction scheme zero solubility must be used with operator splitting."));
             }
 
+          if (fluid_solid_reaction_scheme == tian_approximation)
+            {
+              AssertThrow(this->get_parameters().use_operator_splitting,
+                          ExcMessage("The Fluid-reaction scheme tian approximation must be used with operator splitting."));
+            }
+
           if (this->get_parameters().use_operator_splitting)
             {
               AssertThrow(fluid_reaction_time_scale >= this->get_parameters().reaction_time_step,
@@ -338,7 +467,6 @@ namespace aspect
       // After parsing the parameters for this model, parse parameters related to the base model.
       base_model->parse_parameters(prm);
       this->model_dependence = base_model->get_model_dependence();
-
       if (fluid_solid_reaction_scheme == zero_solubility)
         {
           AssertThrow(this->get_material_model().is_compressible() == false,

tests/tian_MORB.prm

@@ -0,0 +1,25 @@
+# This model generates a representation of a phase diagram for the equilibrium
+# wt% water that can be stored in a MORB as shown in Figure 5 of Tian et. al., 2019.
+# This model reproduces that phase diagram (very low res by default), with an important
+# caveat being that the pressure axis is inverted relative to the ASPECT output.
+# A surface pressure of 600 MPa is imposed to reproduce the pressure axis,
+# and a lateral temperature gradient is imposed to reproduce the temperature axis.
+include $ASPECT_SOURCE_DIR/tests/tian_peridotite.prm
+
+subsection Initial composition model
+  set Model name = function
+  subsection Function
+    set Function expression = initial_porosity; \
+                              initial_bound_water; \
+                              0; \
+                              0; \
+                              1; \
+                              0
+  end
+end
+
+subsection Material model
+  subsection Reactive Fluid Transport Model
+    set Maximum weight percent water in MORB                    = 5.3
+  end
+end