Added options for soil resistance (#193)

Provides a choice for setting the "sellers" soil resistance model, or the default "sakagucki_zeng".

Note that tests don't pass because of a different problem in Amanzi's TPLs, not because of this PR.  Therefore I am overriding the CI...
---------

Co-authored-by: Ethan Coon <coonet@ornl.gov>

src/pks/surface_balance/constitutive_relations/land_cover/LandCover.cc

@@ -40,6 +40,7 @@ LandCover::LandCover(Teuchos::ParameterList& plist)
     water_transition_depth(plist.get<double>("water transition depth [m]", NAN)),
     dessicated_zone_thickness(plist.get<double>("dessicated zone thickness [m]", NAN)),
     clapp_horn_b(plist.get<double>("Clapp and Hornberger b [-]", NAN)),
+    rs_method(plist.get<std::string>("Soil resistance method", "sakagucki_zeng")),
     roughness_ground(plist.get<double>("roughness length of bare ground [m]", NAN)),
     roughness_snow(plist.get<double>("roughness length of snow [m]", NAN)),
     mannings_n(plist.get<double>("Manning's n [?]", NAN))

src/pks/surface_balance/constitutive_relations/land_cover/LandCover.hh

@@ -111,6 +111,9 @@ same region-based partitioning.
       and/or refactored, we must provide this value.  When in doubt, just use
       1.
 
+    * `"Soil resistance method`" ``[string]`` **NAN** This denotes what model to use
+      to calculate soil resistance. Currently support {sakagucki_zeng, sellers}.
+
     * `"roughness length of bare ground [m]`" ``[double]`` **NaN** Roughness length of
       the bare soil, used in calculating sensible/latent heat in the
       physics-based SEB model.  A typical value is 0.04.
@@ -177,6 +180,7 @@ struct LandCover {
   double dessicated_zone_thickness; // [m] Thickness over which vapor must diffuse
                                     //  when the soil is dry.
   double clapp_horn_b;              // [-] exponent of the WRM, Clapp & Hornberger eqn 1
+  std::string rs_method;            // soil resistance method {sakagucki_zeng, sellers, TBD}
   double roughness_ground;          // [m] Fetch length for latent/sensible heat fluxes.
   double roughness_snow;            // [m] Fetch length for latent/sensible heat fluxes.
 };

src/pks/surface_balance/constitutive_relations/land_cover/evaporation_downregulation_evaluator.cc

@@ -41,10 +41,11 @@ EvaporationDownregulationEvaluator::InitializeFromPlist_()
   domain_surf_ = Keys::getDomain(my_keys_.front().first);
   domain_sub_ = Keys::readDomainHint(plist_, domain_surf_, "surface", "domain");
 
-  // sat gas and porosity on subsurface
+  // sat gas, sat liquid, and porosity on subsurface
   sat_gas_key_ = Keys::readKey(plist_, domain_sub_, "saturation gas", "saturation_gas");
   dependencies_.insert(KeyTag{ sat_gas_key_, tag });
-
+  sat_liq_key_ = Keys::readKey(plist_, domain_sub_, "saturation liquid", "saturation_liquid");
+  dependencies_.insert(KeyTag{ sat_liq_key_, tag });
   poro_key_ = Keys::readKey(plist_, domain_sub_, "porosity", "porosity");
   dependencies_.insert(KeyTag{ poro_key_, tag });
 
@@ -62,6 +63,8 @@ EvaporationDownregulationEvaluator::Evaluate_(const State& S,
   Tag tag = my_keys_.front().second;
   const Epetra_MultiVector& sat_gas =
     *S.Get<CompositeVector>(sat_gas_key_, tag).ViewComponent("cell", false);
+  const Epetra_MultiVector& sat_liq =
+    *S.Get<CompositeVector>(sat_liq_key_, tag).ViewComponent("cell", false);
   const Epetra_MultiVector& poro =
     *S.Get<CompositeVector>(poro_key_, tag).ViewComponent("cell", false);
   const Epetra_MultiVector& pot_evap =
@@ -78,7 +81,7 @@ EvaporationDownregulationEvaluator::Evaluate_(const State& S,
     for (AmanziMesh::Entity_ID sc : lc_ids) {
       auto c = sub_mesh.cells_of_column(sc)[0];
       surf_evap[0][sc] =
-        region_model.second->Evaporation(sat_gas[0][c], poro[0][c], pot_evap[0][sc]);
+        region_model.second->Evaporation(sat_gas[0][c], poro[0][c], pot_evap[0][sc], sat_liq[0][c]);
     }
   }
 }
@@ -95,6 +98,8 @@ EvaporationDownregulationEvaluator::EvaluatePartialDerivative_(
   if (wrt_key == pot_evap_key_) {
     const Epetra_MultiVector& sat_gas =
       *S.Get<CompositeVector>(sat_gas_key_, tag).ViewComponent("cell", false);
+    const Epetra_MultiVector& sat_liq =
+      *S.Get<CompositeVector>(sat_liq_key_, tag).ViewComponent("cell", false);
     const Epetra_MultiVector& poro =
       *S.Get<CompositeVector>(poro_key_, tag).ViewComponent("cell", false);
     const Epetra_MultiVector& pot_evap =
@@ -112,7 +117,7 @@ EvaporationDownregulationEvaluator::EvaluatePartialDerivative_(
       for (AmanziMesh::Entity_ID sc : lc_ids) {
         auto c = sub_mesh.cells_of_column(sc)[0];
         surf_evap[0][sc] = region_model.second->DEvaporationDPotentialEvaporation(
-          sat_gas[0][c], poro[0][c], pot_evap[0][sc]);
+          sat_gas[0][c], poro[0][c], pot_evap[0][sc], sat_liq[0][c]);
       }
     }
   }
@@ -136,6 +141,9 @@ EvaporationDownregulationEvaluator::EnsureCompatibility_ToDeps_(State& S)
     S.Require<CompositeVector, CompositeVectorSpace>(sat_gas_key_, tag)
       .SetMesh(S.GetMesh(domain_sub_))
       ->AddComponent("cell", AmanziMesh::Entity_kind::CELL, 1);
+    S.Require<CompositeVector, CompositeVectorSpace>(sat_liq_key_, tag)
+      .SetMesh(S.GetMesh(domain_sub_))
+      ->AddComponent("cell", AmanziMesh::Entity_kind::CELL, 1);
     S.Require<CompositeVector, CompositeVectorSpace>(pot_evap_key_, tag)
       .SetMesh(S.GetMesh(domain_surf_))
       ->AddComponent("cell", AmanziMesh::Entity_kind::CELL, 1);

src/pks/surface_balance/constitutive_relations/land_cover/evaporation_downregulation_evaluator.hh

@@ -23,6 +23,7 @@ Hornberger b.
    KEYS:
 
    - `"saturation gas`" **DOMAIN_SUB-saturation_gas**
+   - `"saturation liquid`" **DOMAIN_SUB-saturation_liquid**
    - `"porosity`" **DOMAIN_SUB-porosity**
    - `"potential evaporation`" **DOMAIN_SUB-potential_evaporation**
 
@@ -60,6 +61,7 @@ class EvaporationDownregulationEvaluator : public EvaluatorSecondaryMonotypeCV {
   void InitializeFromPlist_();
 
   Key sat_gas_key_;
+  Key sat_liq_key_;
   Key poro_key_;
   Key pot_evap_key_;
 

src/pks/surface_balance/constitutive_relations/land_cover/evaporation_downregulation_model.cc

@@ -22,43 +22,72 @@ EvaporationDownregulationModel::EvaporationDownregulationModel(Teuchos::Paramete
 {
   dess_dz_ = plist.get<double>("dessicated zone thickness [m]", 0.1);
   Clapp_Horn_b_ = plist.get<double>("Clapp and Hornberger b of surface soil [-]", 1.0);
+  rs_method_ = plist.get<std::string>("Soil resistance method", "sakagucki_zeng");
 }
 
 // Constructor from LandCover
 EvaporationDownregulationModel::EvaporationDownregulationModel(const LandCover& lc)
 {
   dess_dz_ = lc.dessicated_zone_thickness;
   Clapp_Horn_b_ = lc.clapp_horn_b;
+  rs_method_ = lc.rs_method;
 }
 
 // main method
 double
-EvaporationDownregulationModel::Evaporation(double sg, double poro, double pot_evap) const
+EvaporationDownregulationModel::Evaporation(double sg,
+                                            double poro,
+                                            double pot_evap,
+                                            double sl) const
 {
-  double rsoil = Relations::EvaporativeResistanceCoef(sg, poro, dess_dz_, Clapp_Horn_b_);
+  double rsoil;
+  if (rs_method_ == "sellers") {
+    rsoil = Relations::EvaporativeResistanceCoefSellers(sl);
+  } else {
+    rsoil = Relations::EvaporativeResistanceCoefSakaguckiZeng(sg, poro, dess_dz_, Clapp_Horn_b_);
+  }
   return pot_evap / (1. + rsoil);
 }
 
 double
 EvaporationDownregulationModel::DEvaporationDSaturationGas(double sg,
                                                            double poro,
-                                                           double pot_evap) const
+                                                           double pot_evap,
+                                                           double sl) const
 {
   return 0.;
 }
 
 double
-EvaporationDownregulationModel::DEvaporationDPorosity(double sg, double poro, double pot_evap) const
+EvaporationDownregulationModel::DEvaporationDPorosity(double sg,
+                                                      double poro,
+                                                      double pot_evap,
+                                                      double sl) const
+{
+  return 0.;
+}
+
+double
+EvaporationDownregulationModel::DEvaporationDSaturationLiquid(double sg,
+                                                              double poro,
+                                                              double pot_evap,
+                                                              double sl) const
 {
   return 0.;
 }
 
 double
 EvaporationDownregulationModel::DEvaporationDPotentialEvaporation(double sg,
                                                                   double poro,
-                                                                  double pot_evap) const
+                                                                  double pot_evap,
+                                                                  double sl) const
 {
-  return 1. / (1. + Relations::EvaporativeResistanceCoef(sg, poro, dess_dz_, Clapp_Horn_b_));
+  if (rs_method_ == "sellers") {
+    return 1. / (1. + Relations::EvaporativeResistanceCoefSellers(sl));
+  } else {
+    return 1. / (1. + Relations::EvaporativeResistanceCoefSakaguckiZeng(
+                        sg, poro, dess_dz_, Clapp_Horn_b_));
+  }
 }
 
 } // namespace Relations

src/pks/surface_balance/constitutive_relations/land_cover/evaporation_downregulation_model.hh

@@ -25,6 +25,9 @@ Requires two parameters,
   Genuchten curves that we typically use, but it doesn't appear to be the most
   important parameter.
 
+* `"Soil resistance method`" ``[string]`` Soil resistance model, currently support 
+  {sakagucki_zeng, sellers}, default is always sakagucki_zeng.
+
 These may be provided via parameter list or LandCover type.
 
 */
@@ -43,18 +46,21 @@ class EvaporationDownregulationModel {
   EvaporationDownregulationModel(Teuchos::ParameterList& plist);
   EvaporationDownregulationModel(const LandCover& lc);
 
-  double Evaporation(double sg, double poro, double pot_evap) const;
+  double Evaporation(double sg, double poro, double pot_evap, double sl) const;
 
-  double DEvaporationDSaturationGas(double sg, double poro, double pot_evap) const;
-  double DEvaporationDPorosity(double sg, double poro, double pot_evap) const;
-  double DEvaporationDPotentialEvaporation(double sg, double poro, double pot_evap) const;
+  double DEvaporationDSaturationGas(double sg, double poro, double pot_evap, double sl) const;
+  double DEvaporationDPorosity(double sg, double poro, double pot_evap, double sl) const;
+  double DEvaporationDSaturationLiquid(double sg, double poro, double pot_evap, double sl) const;
+  double
+  DEvaporationDPotentialEvaporation(double sg, double poro, double pot_evap, double sl) const;
 
  protected:
   void InitializeFromPlist_(Teuchos::ParameterList& plist);
 
  protected:
   double dess_dz_;
   double Clapp_Horn_b_;
+  std::string rs_method_;
 };
 
 } // namespace Relations

src/pks/surface_balance/constitutive_relations/land_cover/radiation_balance_evaluator.cc

@@ -167,10 +167,9 @@ RadiationBalanceEvaluator::Evaluate_(const State& S, const std::vector<Composite
       rad_bal_surf[0][c] = (1 - albedo[0][c]) * sw_atm_surf + lw_down - lw_up_surf;
       rad_bal_snow[0][c] = (1 - albedo[1][c]) * sw_atm_surf + lw_down - lw_up_snow;
 
-      rad_bal_can[0][c] =  (1 - lc.second.albedo_canopy) * sw_atm_can
-        + lw_atm_can - 2 * lw_can
-        + area_frac[0][c] * e_can_lw * lw_up_surf
-        + area_frac[1][c] * e_can_lw * lw_up_snow;
+      rad_bal_can[0][c] = (1 - lc.second.albedo_canopy) * sw_atm_can + lw_atm_can - 2 * lw_can +
+                          area_frac[0][c] * e_can_lw * lw_up_surf +
+                          area_frac[1][c] * e_can_lw * lw_up_snow;
     }
   }
 }

src/pks/surface_balance/constitutive_relations/land_cover/seb_physics_defs.hh

@@ -100,9 +100,11 @@ struct GroundProperties {
   double density_w;         // density [kg/m^3]
   double dz;                // [m]
   double clapp_horn_b;      // [-]
+  std::string rs_method;    // options: {sakagucki_zeng, sellers}
   double albedo;            // [-]
   double emissivity;        // [-]
   double saturation_gas;    // [-]
+  double saturation_liq;    // [-]
   double roughness;         // [m] surface roughness of a bare domain
   double snow_death_rate;   // [kg/m^2/s] snow that must die this timestep, make it melt!
   double unfrozen_fraction; // [-] fraction of ground water that is unfrozen
@@ -118,8 +120,10 @@ struct GroundProperties {
       albedo(NaN),
       emissivity(NaN),
       saturation_gas(NaN),
+      saturation_liq(NaN),
       roughness(NaN),
       clapp_horn_b(3),
+      rs_method("sakagucki_zeng"),
       snow_death_rate(0.),
       unfrozen_fraction(0.),
       water_transition_depth(0.01)

src/pks/surface_balance/constitutive_relations/land_cover/seb_physics_funcs.cc

@@ -31,7 +31,6 @@ double
 CalcAlbedoSnow(double density_snow)
 {
   double AlSnow;
-  //  432.23309912785146
   if (density_snow <= 432.23309912785146) {
     AlSnow = 1.0 - 0.247 * std::pow(0.16 + 110 * std::pow(density_snow / 1000, 4), 0.5);
   } else {
@@ -183,16 +182,22 @@ EvaporativeResistanceGround(const GroundProperties& surf,
   if (met.vp_air > vapor_pressure_ground) { // condensation
     return 0.;
   } else {
-    return EvaporativeResistanceCoef(
-      surf.saturation_gas, surf.porosity, surf.dz, surf.clapp_horn_b);
+    if (surf.rs_method == "sakagucki_zeng") {
+      return EvaporativeResistanceCoefSakaguckiZeng(
+        surf.saturation_gas, surf.porosity, surf.dz, surf.clapp_horn_b);
+    } else if (surf.rs_method == "sellers") {
+      return EvaporativeResistanceCoefSellers(surf.saturation_liq);
+    } else {
+      throw("Currently support {sakagucki_zeng, sellers}");
+    }
   }
 }
 
 double
-EvaporativeResistanceCoef(double saturation_gas,
-                          double porosity,
-                          double dessicated_zone_thickness,
-                          double Clapp_Horn_b)
+EvaporativeResistanceCoefSakaguckiZeng(double saturation_gas,
+                                       double porosity,
+                                       double dessicated_zone_thickness,
+                                       double Clapp_Horn_b)
 {
   double Rsoil;
   if (saturation_gas == 0.) {
@@ -220,6 +225,11 @@ EvaporativeResistanceCoef(double saturation_gas,
   return Rsoil;
 }
 
+double
+EvaporativeResistanceCoefSellers(double saturation_liq)
+{
+  return std::exp(8.206 - 4.255 * saturation_liq);
+}
 
 double
 SensibleHeat(double resistance_coef,

src/pks/surface_balance/constitutive_relations/land_cover/seb_physics_funcs.hh

@@ -126,11 +126,13 @@ EvaporativeResistanceGround(const GroundProperties& surf,
                             double vapor_pressure_ground);
 
 double
-EvaporativeResistanceCoef(double saturation_gas,
-                          double porosity,
-                          double dessicated_zone_thickness,
-                          double Clapp_Horn_b);
+EvaporativeResistanceCoefSakaguckiZeng(double saturation_gas,
+                                       double porosity,
+                                       double dessicated_zone_thickness,
+                                       double Clapp_Horn_b);
 
+double
+EvaporativeResistanceCoefSellers(double saturation_liq);
 
 //
 // Basic sensible heat.

src/pks/surface_balance/constitutive_relations/land_cover/seb_threecomponent_evaluator.cc

@@ -156,6 +156,8 @@ SEBThreeComponentEvaluator::SEBThreeComponentEvaluator(Teuchos::ParameterList& p
   // -- subsurface properties for evaporating bare soil
   sat_gas_key_ = Keys::readKey(plist, domain_ss_, "gas saturation", "saturation_gas");
   dependencies_.insert(KeyTag{ sat_gas_key_, tag });
+  sat_liq_key_ = Keys::readKey(plist, domain_ss_, "liquid saturation", "saturation_liquid");
+  dependencies_.insert(KeyTag{ sat_liq_key_, tag });
   poro_key_ = Keys::readKey(plist, domain_ss_, "porosity", "porosity");
   dependencies_.insert(KeyTag{ poro_key_, tag });
   ss_pres_key_ = Keys::readKey(plist, domain_ss_, "subsurface pressure", "pressure");
@@ -208,6 +210,7 @@ SEBThreeComponentEvaluator::Evaluate_(const State& S, const std::vector<Composit
 
   // collect subsurface properties
   const auto& sat_gas = *S.Get<CompositeVector>(sat_gas_key_, tag).ViewComponent("cell", false);
+  const auto& sat_liq = *S.Get<CompositeVector>(sat_liq_key_, tag).ViewComponent("cell", false);
   const auto& poro = *S.Get<CompositeVector>(poro_key_, tag).ViewComponent("cell", false);
   const auto& ss_pres = *S.Get<CompositeVector>(ss_pres_key_, tag).ViewComponent("cell", false);
 
@@ -284,11 +287,13 @@ SEBThreeComponentEvaluator::Evaluate_(const State& S, const std::vector<Composit
         surf.density_w = mass_dens[0][c];
         surf.dz = lc.second.dessicated_zone_thickness;
         surf.clapp_horn_b = lc.second.clapp_horn_b;
+        surf.rs_method = lc.second.rs_method;
         surf.albedo = sg_albedo[0][c];
         surf.emissivity = emissivity[0][c];
         surf.ponded_depth = 0.; // by definition
         surf.porosity = poro[0][cells[0]];
         surf.saturation_gas = sat_gas[0][cells[0]];
+        surf.saturation_liq = sat_liq[0][cells[0]];
         surf.unfrozen_fraction = unfrozen_fraction[0][c];
         surf.water_transition_depth = lc.second.water_transition_depth;
 
@@ -355,11 +360,13 @@ SEBThreeComponentEvaluator::Evaluate_(const State& S, const std::vector<Composit
         surf.density_w = mass_dens[0][c];
         surf.dz = lc.second.dessicated_zone_thickness;
         surf.clapp_horn_b = lc.second.clapp_horn_b;
+        surf.rs_method = lc.second.rs_method;
         surf.emissivity = emissivity[1][c];
         surf.albedo = sg_albedo[1][c];
         surf.ponded_depth = std::max(lc.second.water_transition_depth, ponded_depth[0][c]);
         surf.porosity = 1.;
         surf.saturation_gas = 0.;
+        surf.saturation_liq = sat_liq[0][cells[0]];
         surf.unfrozen_fraction = unfrozen_fraction[0][c];
         surf.water_transition_depth = lc.second.water_transition_depth;
 
@@ -426,10 +433,12 @@ SEBThreeComponentEvaluator::Evaluate_(const State& S, const std::vector<Composit
         surf.density_w = mass_dens[0][c];
         surf.dz = lc.second.dessicated_zone_thickness;
         surf.clapp_horn_b = lc.second.clapp_horn_b;
+        surf.rs_method = lc.second.rs_method; // does not matter
         surf.emissivity = emissivity[2][c];
         surf.albedo = sg_albedo[2][c];
         surf.ponded_depth = 0;                            // does not matter
         surf.saturation_gas = 0.;                         // does not matter
+        surf.saturation_liq = sat_liq[0][cells[0]];       // does not matter
         surf.porosity = 1.;                               // does not matter
         surf.unfrozen_fraction = unfrozen_fraction[0][c]; // does not matter
         surf.water_transition_depth = lc.second.water_transition_depth;
@@ -634,7 +643,8 @@ SEBThreeComponentEvaluator::EnsureCompatibility_ToDeps_(State& S)
                                    "roughness_ground",
                                    "water_transition_depth",
                                    "dessicated_zone_thickness",
-                                   "clapp_horn_b" });
+                                   "clapp_horn_b",
+                                   "rs_method" });
 
     // use domain name to set the mesh type
     CompositeVectorSpace domain_fac;