Gaob/relperm (#191)

Significant changes to improve hillslope permafrost runs:

(1) fixes a minor bug in depression storage evaluator
(2) add the "magic omega" kr model from Niu and Yang (2006) which allows more water to infiltrate during frozen conditions, resulting in higher infiltration and less runoff during snowmelt
(3) revised air resistance, but set default coeffs to 0. So while this does not change the default model, there are more options there now.
(4) Adds options to use the Sellers model of evaporative resistance vs the previous model of Sakagucki & Zeng model, keeping S&Z as the default. The Sellers model has no free parameter, but may work well in some cases.

---------

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

src/pks/flow/constitutive_relations/elevation/volumetric_snow_ponded_depth_evaluator.cc

@@ -194,7 +194,8 @@ VolumetricSnowPondedDepthEvaluator::EnsureCompatibility_ToDeps_(State& S)
 
     // Loop over my dependencies, ensuring they meet the requirements.
     for (const auto& key_tag : dependencies_) {
-      if (key_tag.first == delta_max_key_ || key_tag.first == delta_ex_key_) {
+      if (key_tag.first == delta_max_key_ || key_tag.first == delta_ex_key_ ||
+          key_tag.first == sd_key_) {
         S.Require<CompositeVector, CompositeVectorSpace>(key_tag.first, key_tag.second)
           .Update(*no_bf_dep_fac);
       } else {

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

@@ -0,0 +1,128 @@
+/*
+  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: Ethan Coon (ecoon@lanl.gov)
+           Bo Gao (gaob@ornl.gov)
+*/
+
+//! Evaluates relative permeability using an empirical model for frozen conditions.
+/*!
+
+This is an empirical relative permeability model according to Niu and Yang (2006). 
+This model is based on Brooks-Corey relative permeability model and an additional
+coefficient term is added to account for the effect of soil ice. This model is 
+used for freezing conditions to make snowmelt water infiltrate deeper. See paper
+Agnihotri et al. (2023) for discussions about the influence of relative permeability 
+model on discharge under freezing conditions.
+
+.. math::
+   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.
+
+.. _rel-perm-evaluator-spec
+.. admonition:: rel-perm-evaluator-spec
+
+   * `"use density on viscosity in rel perm`" ``[bool]`` **true**
+
+   * `"boundary rel perm strategy`" ``[string]`` **boundary pressure** Controls
+     how the rel perm is calculated on boundary faces.  Note, this may be
+     overwritten by upwinding later!  One of:
+
+      - `"boundary pressure`" Evaluates kr of pressure on the boundary face, upwinds normally.
+      - `"interior pressure`" Evaluates kr of the pressure on the interior cell (bad idea).
+      - `"harmonic mean`" Takes the harmonic mean of kr on the boundary face and kr on the interior cell.
+      - `"arithmetic mean`" Takes the arithmetic mean of kr on the boundary face and kr on the interior cell.
+      - `"one`" Sets the boundary kr to 1.
+      - `"surface rel perm`" Looks for a field on the surface mesh and uses that.
+
+   * `"minimum rel perm cutoff`" ``[double]`` **0.** Provides a lower bound on rel perm.
+
+   * `"permeability rescaling`" ``[double]`` Typically rho * kr / mu is very big
+     and K_sat is very small.  To avoid roundoff propagation issues, rescaling
+     this quantity by offsetting and equal values is encourage.  Typically 10^7 or so is good.
+
+   * `"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.
+
+   KEYS:
+
+   - `"rel perm`"
+   - `"saturation_liquid`"
+   - `"saturation_gas`"
+   - `"density`" (if `"use density on viscosity in rel perm`" == true)
+   - `"viscosity`" (if `"use density on viscosity in rel perm`" == true)
+   - `"surface relative permeability`" (if `"boundary rel perm strategy`" == `"surface rel perm`")
+
+*/
+
+#pragma once
+
+#include "wrm.hh"
+#include "wrm_partition.hh"
+#include "rel_perm_evaluator.hh"
+#include "EvaluatorSecondaryMonotype.hh"
+#include "Factory.hh"
+
+namespace Amanzi {
+namespace Flow {
+
+class RelPermFrzBCEvaluator : public EvaluatorSecondaryMonotypeCV {
+ public:
+  // constructor format for all derived classes
+  explicit RelPermFrzBCEvaluator(Teuchos::ParameterList& plist);
+
+  RelPermFrzBCEvaluator(Teuchos::ParameterList& plist, const Teuchos::RCP<WRMPartition>& wrms);
+
+  RelPermFrzBCEvaluator(const RelPermFrzBCEvaluator& other) = default;
+  virtual Teuchos::RCP<Evaluator> Clone() const override;
+
+  Teuchos::RCP<WRMPartition> get_WRMs() { return wrms_; }
+
+ protected:
+  virtual void EnsureCompatibility_ToDeps_(State& S) override;
+
+  // Required methods from EvaluatorSecondaryMonotypeCV
+  virtual void Evaluate_(const State& S, const std::vector<CompositeVector*>& result) override;
+  virtual void EvaluatePartialDerivative_(const State& S,
+                                          const Key& wrt_key,
+                                          const Tag& wrt_tag,
+                                          const std::vector<CompositeVector*>& result) override;
+
+ protected:
+  void InitializeFromPlist_();
+
+  Teuchos::RCP<WRMPartition> wrms_;
+  Key sat_key_;
+  Key sat_gas_key_;
+  Key dens_key_;
+  Key visc_key_;
+  Key surf_rel_perm_key_;
+
+  bool is_dens_visc_;
+  Key surf_domain_;
+  BoundaryRelPerm boundary_krel_;
+
+  double perm_scale_;
+  double min_val_;
+  double omega_;
+  double b_;
+
+ private:
+  static Utils::RegisteredFactory<Evaluator, RelPermFrzBCEvaluator> factory_;
+};
+
+} // namespace Flow
+} // namespace Amanzi
+

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

@@ -0,0 +1,20 @@
+/*
+  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: Ethan Coon (ecoon@lanl.gov), Bo Gao (gaob@ornl.gov)
+*/
+
+#include "rel_perm_frzBC_evaluator.hh"
+
+namespace Amanzi {
+namespace Flow {
+
+// registry of method
+Utils::RegisteredFactory<Evaluator, RelPermFrzBCEvaluator>
+  RelPermFrzBCEvaluator::factory_("Brooks-Corey based high frozen rel perm");
+
+} // namespace Flow
+} // namespace Amanzi

src/pks/flow/permafrost_pk.cc

@@ -21,6 +21,7 @@
 #include "wrm_permafrost_evaluator.hh"
 #include "rel_perm_evaluator.hh"
 #include "rel_perm_sutraice_evaluator.hh"
+#include "rel_perm_frzBC_evaluator.hh"
 #include "pk_helpers.hh"
 
 #include "permafrost.hh"

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

@@ -40,7 +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")),
+    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,7 +111,7 @@ 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
+    * `"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

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

@@ -44,22 +44,32 @@ static const double c_R_ideal_gas = 461.52;                 // [Pa m^3 kg^-1 K^-
 // grid-cell by grid-cell basis.
 struct ModelParams {
   ModelParams()
-    : density_air(1.275),          // [kg/m^3]
-      density_freshsnow(100.),     // [kg/m^3]
-      density_frost(200.),         // [kg/m^3]
-      thermalK_freshsnow(0.029),   // thermal conductivity of fresh snow [W/m-K]
-      thermalK_snow_exp(2),        // exponent in thermal conductivity of snow model [-]
-      H_fusion(333500.0),          // Heat of fusion for melting snow -- [J/kg]
-      H_sublimation(2834000.),     // Latent heat of sublimation ------- [J/kg]
-      H_vaporization(2497848.),    // Latent heat of vaporization ------ [J/kg]
-      Cp_air(1004.0),              // Specific heat of air @ const pres- [J/K kg]
-      Cv_water(4218.636),          // Specific heat of water ----------- [J/K kg]
-      P_atm(101325.),              // atmospheric pressure ------------- [Pa]
-      gravity(9.807),              // gravity [kg m / s^2]
-      evap_transition_width(100.), // transition on evaporation from surface to
-                                   // evaporation from subsurface [m],
-                                   // THIS IS DEPRECATED
-      KB(0.)
+    : density_air(1.275),            // [kg/m^3]
+      dynamic_viscosity_air(1.7e-5), // [Pa s]
+      density_freshsnow(100.),       // [kg/m^3]
+      density_frost(200.),           // [kg/m^3]
+      thermalK_freshsnow(0.029),     // thermal conductivity of fresh snow [W/m-K]
+      thermalK_snow_exp(2),          // exponent in thermal conductivity of snow model [-]
+      H_fusion(333500.0),            // Heat of fusion for melting snow -- [J/kg]
+      H_sublimation(2834000.),       // Latent heat of sublimation ------- [J/kg]
+      H_vaporization(2497848.),      // Latent heat of vaporization ------ [J/kg]
+      Cp_air(1004.0),                // Specific heat of air @ const pres- [J/K kg]
+      Cv_water(4218.636),            // Specific heat of water ----------- [J/K kg]
+      P_atm(101325.),                // atmospheric pressure ------------- [Pa]
+      gravity(9.807),                // gravity [kg m / s^2]
+      evap_transition_width(100.),   // transition on evaporation from surface to
+                                     // evaporation from subsurface [m],
+                                     // THIS IS DEPRECATED
+      KB(0.),                        // a paramter denoting the log ratio of momentum roughness length
+                                     // to vapor roughness length for vapor flux or 
+                                     // to thermal roughness length for sensible heat flux.
+      // Da0_a, Da0_b, Cd0_c, Cd0_d are fitting parameters in the formulization of the log ratio 
+      // of momentum roughness length to the vapor roughness length for vapor flux [-]. Setting
+      // Da0_a = Cd0_c = Cd0_d = 0 means momentum roughness length = vapor roughness length.
+      Da0_a(0.),
+      Da0_b(0.1),
+      Cd0_c(0.),
+      Cd0_d(0.)
   {}
 
   ModelParams(Teuchos::ParameterList& plist) : ModelParams()
@@ -72,19 +82,25 @@ struct ModelParams {
       plist.get<double>("evaporation transition width [Pa]", evap_transition_width);
 
     KB = plist.get<double>("log ratio between z0m and z0h [-]", KB);
+    Da0_a = plist.get<double>("coefficient a for Da0 [-]", Da0_a);
+    Da0_b = plist.get<double>("exponent b for Da0 [-]", Da0_b);
+    Cd0_c = plist.get<double>("coefficient c for Cd0 [-]", Cd0_c);
+    Cd0_d = plist.get<double>("coefficient d for Cd0 [-]", Cd0_d);
   }
 
   // likely constants
   double gravity;
   double P_atm;
   double density_air;
+  double dynamic_viscosity_air;
   double density_freshsnow;
   double density_frost;
   double thermalK_freshsnow;
   double thermalK_snow_exp;
   double H_fusion, H_sublimation, H_vaporization;
   double Cp_air, Cv_water;
   double KB;
+  double Da0_a, Da0_b, Cd0_c, Cd0_d;
 
   // other parameters
   double evap_transition_width;

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

@@ -290,9 +290,26 @@ UpdateEnergyBalanceWithSnow_Inner(const GroundProperties& surf,
 
   // latent heat
   double vapor_pressure_skin = SaturatedVaporPressure(snow.temp);
-
+  // KB here is to formulize the log ration between momentum roughness length and vapor roughness length
+  // to add the effect of surface microtopography. Details see Gao et al. 2021 (WRR) Eq.(17) and similar
+  // studies by Mölder & Lindroth 2001 (Agricultural and Forest Meteorology) but in heat transport. 
+  // Theories see Chapter 4 of book [Brutsaert, W. (1982). Evaporation into the atmosphere: Theory, 
+  // history, and applications]. Da0_a, Da0_b, Cd0_c, Cd0_d here are four fitting coefficients in the 
+  // formulization. There are several different values proposed by different studies (see Brutsaert 1982). 
+  // In Gao et al. (2021), these four coefficients are determined by lab experiments. 
+  // If set Da0_a, Cd0_c, Cd0_d to 0, it means that vapor roughness length is assumed equal to momentum 
+  // roughness length. Currently, this approach still needs testing or other formulization approaches 
+  // may also be added later for testing. So keep Da0_a = Cd0_c = Cd0_d = 0 for users.
+  double u_star =
+    met.Us * c_von_Karman /
+    std::log(met.Z_Us / CalcRoughnessFactor(snow.height, surf.roughness, snow.roughness));
+  double Re0 = params.density_air * u_star *
+               CalcRoughnessFactor(snow.height, surf.roughness, snow.roughness) /
+               params.dynamic_viscosity_air;
+  double KB =
+    (params.Da0_a * std::pow(Re0, params.Da0_b) - (params.Cd0_c * std::log(Re0) + params.Cd0_d)) * c_von_Karman;
   double Dhe_latent = WindFactor(
-    met.Us, met.Z_Us, CalcRoughnessFactor(snow.height, surf.roughness, snow.roughness), 0.);
+    met.Us, met.Z_Us, CalcRoughnessFactor(snow.height, surf.roughness, snow.roughness), KB);
   eb.fQe = LatentHeat(Dhe_latent * Sqig,
                       params.density_air,
                       params.H_sublimation,
@@ -366,7 +383,11 @@ UpdateEnergyBalanceWithoutSnow(const GroundProperties& surf,
 
   // latent heat
   double vapor_pressure_skin = VaporPressureGround(surf, params);
-  double Dhe_latent = WindFactor(met.Us, met.Z_Us, surf.roughness, 0.);
+  double u_star = met.Us * c_von_Karman / std::log(met.Z_Us / surf.roughness);
+  double Re0 = params.density_air * u_star * surf.roughness / params.dynamic_viscosity_air;
+  double KB =
+    (params.Da0_a * std::pow(Re0, params.Da0_b) - (params.Cd0_c * std::log(Re0) + params.Cd0_d)) * c_von_Karman;
+  double Dhe_latent = WindFactor(met.Us, met.Z_Us, surf.roughness, KB);
   double Rsoil = EvaporativeResistanceGround(surf, met, params, vapor_pressure_skin);
   double coef = 1.0 / (Rsoil + 1.0 / (Dhe_latent * Sqig));
 

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

@@ -85,7 +85,7 @@ BeersLawAbsorptivity(double k_extinction, double lai);
 // Wind speed term D_he
 // ------------------------------------------------------------------------------------------
 double
-WindFactor(double Us, double Z_Us, double Z_rough, double c_von_Karman, double KB);
+WindFactor(double Us, double Z_Us, double Z_rough, double KB);
 
 //
 // Stability of convective overturning term Zeta AKA Sqig

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

@@ -436,9 +436,9 @@ SEBThreeComponentEvaluator::Evaluate_(const State& S, const std::vector<Composit
         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.ponded_depth = 0;    // does not matter
+        surf.saturation_gas = 0.; // does not matter
+        surf.saturation_liq = sat_liq[0][cells[0]];
         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;