Add Brooks-Corey model

src/pks/flow/constitutive_relations/wrm/rel_perm_brooks_corey_freezing_coeff.hh

@@ -0,0 +1,41 @@
+/*
+  Copyright 2010-202x held jointly by participating institutions.
+  ATS is released under the three-clause BSD License.
+  The terms of use and "as is" disclaimer for this license are
+  provided in the top-level COPYRIGHT file.
+
+  Authors: Bo Gao (gaob@ornl.gov)
+*/
+
+#pragma once
+#include <cmath>
+
+namespace Amanzi {
+namespace Flow {
+namespace BrooksCoreyFrzCoef {
+
+
+inline double
+frzcoef(double sl, double sg, double omega)
+{
+  return 1. - std::exp(-omega * (sl + sg)) + std::exp(-omega);
+}
+
+
+inline double
+d_frzcoef_dsl(double sl, double sg, double omega)
+{
+  return omega * std::exp(-omega * (sl + sg));
+}
+
+
+inline double
+d_frzcoef_dsg(double sl, double sg, double omega)
+{
+  return omega * std::exp(-omega * (sl + sg));
+}
+
+
+} // namespace BrooksCoreyFrzCoef
+} // namespace Flow
+} // namespace Amanzi

src/pks/flow/constitutive_relations/wrm/rel_perm_frzBC_evaluator.cc

@@ -9,6 +9,7 @@
 */
 
 //! Evaluates relative permeability using an empirical model for frozen conditions.
+#include "rel_perm_brooks_corey_freezing_coeff.hh"
 #include "rel_perm_frzBC_evaluator.hh"
 
 namespace Amanzi {
@@ -103,7 +104,6 @@ RelPermFrzBCEvaluator::InitializeFromPlist_()
   min_val_ = plist_.get<double>("minimum rel perm cutoff", 0.);
   perm_scale_ = plist_.get<double>("permeability rescaling");
   omega_ = plist_.get<double>("omega [-]", 2.0);
-  b_ = plist_.get<double>("Clapp and Hornberger b [-]", 2.0);
 }
 
 
@@ -157,9 +157,8 @@ RelPermFrzBCEvaluator::Evaluate_(const State& S, const std::vector<CompositeVect
   for (unsigned int c = 0; c != ncells; ++c) {
     int index = (*wrms_->first)[c];
     double sat_res = wrms_->second[index]->residualSaturation();
-    double coef = 1 - std::exp(-omega_ * (sat_c[0][c] + sat_gas_c[0][c])) + std::exp(-omega_);
-    res_c[0][c] = std::max(
-      coef * std::pow((1 - sat_gas_c[0][c] - sat_res) / (1 - sat_res), 3 + 2 * b_), min_val_);
+    double coef = BrooksCoreyFrzCoef::frzcoef(sat_c[0][c], sat_gas_c[0][c], omega_);
+    res_c[0][c] = std::max(wrms_->second[index]->k_relative(1. - sat_gas_c[0][c]) * coef,  min_val_);
   }
 
   // -- Potentially evaluate the model on boundary faces as well.
@@ -186,34 +185,29 @@ RelPermFrzBCEvaluator::Evaluate_(const State& S, const std::vector<CompositeVect
       int index = (*wrms_->first)[cells[0]];
       double sat_res = wrms_->second[index]->residualSaturation();
       double krel;
-      double coef_b =
-        1 - std::exp(-omega_ * (sat_bf[0][bf] + sat_gas_bf[0][bf])) + std::exp(-omega_);
-      double coef_c =
-        1 - std::exp(-omega_ * (sat_c[0][cells[0]] + sat_gas_c[0][cells[0]])) + std::exp(-omega_);
+
+      double coef_b = BrooksCoreyFrzCoef::frzcoef(sat_bf[0][bf], sat_gas_bf[0][bf], omega_);
+      double coef_c =BrooksCoreyFrzCoef::frzcoef(sat_c[0][cells[0]], sat_gas_c[0][cells[0]], omega_);
+
       if (boundary_krel_ == BoundaryRelPerm::HARMONIC_MEAN) {
         double krelb =
-          std::max(coef_b * std::pow((1 - sat_gas_bf[0][bf] - sat_res) / (1 - sat_res), 3 + 2 * b_),
-                   min_val_);
+          std::max(wrms_->second[index]->k_relative(1. - sat_gas_bf[0][bf]) * coef_b, min_val_);
         double kreli = std::max(
-          coef_c * std::pow((1 - sat_gas_c[0][cells[0]] - sat_res) / (1 - sat_res), 3 + 2 * b_),
-          min_val_);
+          wrms_->second[index]->k_relative(1. - sat_gas_c[0][cells[0]]) * coef_c, min_val_);
         krel = 1.0 / (1.0 / krelb + 1.0 / kreli);
       } else if (boundary_krel_ == BoundaryRelPerm::ARITHMETIC_MEAN) {
         double krelb =
-          std::max(coef_b * std::pow((1 - sat_gas_bf[0][bf] - sat_res) / (1 - sat_res), 3 + 2 * b_),
-                   min_val_);
+          std::max(wrms_->second[index]->k_relative(1. - sat_gas_bf[0][bf]) * coef_b, min_val_);
         double kreli = std::max(
-          coef_c * std::pow((1 - sat_gas_c[0][cells[0]] - sat_res) / (1 - sat_res), 3 + 2 * b_),
-          min_val_);
+          wrms_->second[index]->k_relative(1. - sat_gas_c[0][cells[0]]) * coef_c, min_val_);
         krel = (krelb + kreli) / 2.0;
       } else if (boundary_krel_ == BoundaryRelPerm::INTERIOR_PRESSURE) {
-        krel = std::max(
-          coef_c * std::pow((1 - sat_gas_c[0][cells[0]] - sat_res) / (1 - sat_res), 3 + 2 * b_),
-          min_val_);
+        krel = std::max(wrms_->second[index]->k_relative(1. - sat_gas_c[0][cells[0]]) * coef_c,
+                        min_val_);
       } else if (boundary_krel_ == BoundaryRelPerm::ONE) {
         krel = 1.;
       } else {
-        krel = coef_b * std::pow((1 - sat_gas_bf[0][bf] - sat_res) / (1 - sat_res), 3 + 2 * b_);
+        krel = wrms_->second[index]->k_relative(1. - sat_gas_bf[0][bf]) * coef_b;
       }
       res_bf[0][bf] = std::max(krel, min_val_);
     }
@@ -301,9 +295,8 @@ RelPermFrzBCEvaluator::EvaluatePartialDerivative_(const State& S,
     for (unsigned int c = 0; c != ncells; ++c) {
       int index = (*wrms_->first)[c];
       double sat_res = wrms_->second[index]->residualSaturation();
-      res_c[0][c] = omega_ *
-                    std::pow((1. - sat_gas_c[0][c] - sat_res) / (1 - sat_res), 2 * b_ + 3) *
-                    std::exp(-omega_ * (sat_c[0][c] + sat_gas_c[0][c]));
+      double dcoef_dsl = BrooksCoreyFrzCoef::d_frzcoef_dsl(sat_c[0][c], sat_gas_c[0][c], omega_);
+      res_c[0][c] = dcoef_dsl * wrms_->second[index]->k_relative(1. - sat_gas_c[0][c]);
       AMANZI_ASSERT(res_c[0][c] >= 0.);
     }
 
@@ -345,13 +338,10 @@ RelPermFrzBCEvaluator::EvaluatePartialDerivative_(const State& S,
     for (unsigned int c = 0; c != ncells; ++c) {
       int index = (*wrms_->first)[c];
       double sat_res = wrms_->second[index]->residualSaturation();
-      double dbc_dsg = -(2 * b_ + 3) / (1 - sat_res) *
-                       std::pow((1 - sat_gas_c[0][c] - sat_res) / (1 - sat_res), 2 * b_ + 2);
-      double coef = 1 - std::exp(-omega_ * (sat_c[0][c] + sat_gas_c[0][c])) + std::exp(-omega_);
-      double dcoef_dsg = omega_ * std::exp(-omega_ * (sat_c[0][c] + sat_gas_c[0][c]));
-      res_c[0][c] =
-        dbc_dsg * coef +
-        dcoef_dsg * std::pow((1. - sat_gas_c[0][c] - sat_res) / (1 - sat_res), 2 * b_ + 3);
+      double coef = BrooksCoreyFrzCoef::frzcoef(sat_c[0][c], sat_gas_c[0][c], omega_);
+      double dcoef_dsg = BrooksCoreyFrzCoef::d_frzcoef_dsg(sat_c[0][c], sat_gas_c[0][c], omega_);
+      res_c[0][c] = -coef * wrms_->second[index]->d_k_relative(1. - sat_gas_c[0][c]) +
+                    dcoef_dsg * wrms_->second[index]->k_relative(1. - sat_gas_c[0][c]);
     }
 
     // -- Potentially evaluate the model on boundary faces as well.

src/pks/flow/constitutive_relations/wrm/rel_perm_frzBC_evaluator.hh

@@ -22,10 +22,13 @@ model on discharge under freezing conditions.
    k_{rel} = ( 1 - F_{frz} ) \times ( \frac{1 - s_{g} - s_r}{1 - s_r} )^{2*b + 3} \\
    F_{frz} = \mathrm{exp}( -\omega \times ( s_{l} + s_{g} ) ) - \mathrm{exp}( -\omega )
 
-Note this implementation is currently a bit inconsistent in that it uses WRMs
-to get the residual saturation, and uses a global value for b, the Clapp and
-Hornberger parameter (equivalent to 1/lambda in Brooks & Corey).  This will be
-fixed shortly (see ticket #196) to make b spatially variable by WRM.
+Under freezing conditions, it is recommended to call Brooks-Corey based relative 
+permeability corrected by ice content. This model needs Brooks-Corey parameters:
+Brooks-Corey lambda, Brooks-Corey saturated matric suction (Pa), and residual 
+saturation. The reciprocal of Brooks-Corey lambda is Clapp-Hornberger b. Use tool
+`"convert_paramters_vg2bc.py`" to convert van Genuchten parameters to Brooks-Corey 
+paramters. The conversion method is referred to Lenhard et al. (1989) or Ma et al. (1999)
+method 2. 
 
 .. _rel-perm-evaluator-spec
 .. admonition:: rel-perm-evaluator-spec
@@ -52,13 +55,10 @@ fixed shortly (see ticket #196) to make b spatially variable by WRM.
    * `"omega [-]`" ``[double]`` **2.0** A scale dependent parameter in the relative permeability model. 
      See paper Niu & Yang (2006) for details about the model. Typical values range from 2-3.
 
-   * `"Clapp and Hornberger b [-]`" ``[double]`` **2.0** Clapp-Hornberger b parameter.
-
-   * `"WRM parameters`" ``[wrm-typedinline-spec-list]``  List (by region) of WRM specs.
+   * `"WRM parameters`" ``[wrm-typedinline-spec-list]`` List (by region) of WRM specs.
 
    KEYS:
 
-   - `"rel perm`"
    - `"saturation_liquid`"
    - `"saturation_gas`"
    - `"density`" (if `"use density on viscosity in rel perm`" == true)
@@ -117,7 +117,6 @@ class RelPermFrzBCEvaluator : public EvaluatorSecondaryMonotypeCV {
   double perm_scale_;
   double min_val_;
   double omega_;
-  double b_;
 
  private:
   static Utils::RegisteredFactory<Evaluator, RelPermFrzBCEvaluator> factory_;

src/pks/flow/constitutive_relations/wrm/wrm_brooks_corey.cc

@@ -0,0 +1,134 @@
+/*
+  Copyright 2010-202x held jointly by participating institutions.
+  ATS is released under the three-clause BSD License.
+  The terms of use and "as is" disclaimer for this license are
+  provided in the top-level COPYRIGHT file.
+
+  Authors: Konstantin Lipnikov (lipnikov@lanl.gov)
+           Bo Gao (gaob@ornl.gov)
+*/
+
+#include <cmath>
+#include "dbc.hh"
+#include "errors.hh"
+#include "Spline.hh"
+
+#include "wrm_brooks_corey.hh"
+
+namespace Amanzi {
+namespace Flow {
+
+/* ******************************************************************
+ * Setup fundamental parameters for this model.
+ ****************************************************************** */
+WRMBrooksCorey::WRMBrooksCorey(Teuchos::ParameterList& plist) : plist_(plist)
+{
+  InitializeFromPlist_();
+};
+
+
+void 
+WRMBrooksCorey::InitializeFromPlist_() 
+{
+  lambda_ = plist_.get<double>("Brooks Corey lambda [-]");
+  b_ = 1. / lambda_; // clapp hornberger b
+  p_sat_ = plist_.get<double>("Brooks Corey saturated matric suction [Pa]");
+  sr_ = plist_.get<double>("residual saturation [-]", 0.0);
+
+  s0_ = 1.0 - plist_.get<double>("smoothing interval width [saturation]", 0.0);
+  if (s0_ < 1.) { fit_kr_.Setup(s0_, k_relative(s0_), d_k_relative(s0_), 1.0, 1.0, 0.0); }
+}
+
+
+/* ******************************************************************
+ * Relative permeability formula: input is saturation.
+ * The original curve is regulized on interval (s0, 1) using the
+ * Hermite interpolant of order 3. Formulas (3.11)-(3.12).
+ ****************************************************************** */
+double 
+WRMBrooksCorey::k_relative(double s)
+{
+  if (s <= s0_) {
+    double se = (s - sr_) / (1 - sr_);
+    return pow(se, 2 * b_ + 3);
+  } else if (s == 1.) {
+    return 1.;
+  } else {
+    return fit_kr_(s);
+  }
+}
+
+
+/* ******************************************************************
+ * D Relative permeability / D saturation
+ ****************************************************************** */
+double 
+WRMBrooksCorey::d_k_relative(double s)
+{
+  if (s <= s0_) {
+    double se = (s - sr_) / (1 - sr_);
+    double dkdse = (2 * b_ + 3) * pow(se, 2 * b_ + 2);
+    return dkdse / (1. - sr_);
+  } else if (s == 1.) {
+    return 0.0;
+  } else{
+    return fit_kr_.Derivative(s);
+  }
+}
+
+
+/* ******************************************************************
+ * Saturation formula (3.5)-(3.8).
+ ****************************************************************** */
+double 
+WRMBrooksCorey::saturation(double pc)
+{
+  if (pc <= p_sat_) {
+    return 1.;
+  } else {
+    return pow(p_sat_ / pc, lambda_) * (1. - sr_) + sr_;
+  }
+}
+
+
+/* ******************************************************************
+ * Derivative of the saturation formula w.r.t. capillary pressure.
+ ****************************************************************** */
+double 
+WRMBrooksCorey::d_saturation(double pc)
+{
+  if ( pc <= p_sat_) {
+    return 0.;
+  } else {
+    return -(1. - sr_) * lambda_ * pow(p_sat_ / pc, lambda_) / pc;
+  }
+}
+
+
+/* ******************************************************************
+ * Pressure as a function of saturation, formula (3.9).
+ ****************************************************************** */
+double 
+WRMBrooksCorey::capillaryPressure(double s)
+{
+  double se = (s - sr_) / (1. - sr_);
+  se = std::min<double>(se, 1.0);
+  se = std::max<double>(se, 1.e-40);
+  return pow(se, -b_) * p_sat_;
+}
+
+
+/* ******************************************************************
+ * Derivative of pressure formulat w.r.t. saturation.
+ ****************************************************************** */
+double 
+WRMBrooksCorey::d_capillaryPressure(double s)
+{
+  double se = (s - sr_) / (1. - sr_);
+  se = std::min<double>(se, 1.0);
+  se = std::max<double>(se, 1.e-40);
+  return -b_ * p_sat_ / (1. - sr_) * pow(se, -b_ - 1.);
+}
+
+} // namespace Flow
+} // namespace Amanzi