From 7610da77e26a132007fa5b865092a4c85220c611 Mon Sep 17 00:00:00 2001 From: Bo Gao <83304414+gaobhub@users.noreply.github.com> Date: Fri, 30 Jun 2023 16:00:35 -0400 Subject: [PATCH] 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 --- .../land_cover/LandCover.cc | 1 + .../land_cover/LandCover.hh | 4 ++ .../evaporation_downregulation_evaluator.cc | 16 ++++++-- .../evaporation_downregulation_evaluator.hh | 2 + .../evaporation_downregulation_model.cc | 41 ++++++++++++++++--- .../evaporation_downregulation_model.hh | 14 +++++-- .../land_cover/radiation_balance_evaluator.cc | 7 ++-- .../land_cover/seb_physics_defs.hh | 4 ++ .../land_cover/seb_physics_funcs.cc | 24 +++++++---- .../land_cover/seb_physics_funcs.hh | 10 +++-- .../seb_threecomponent_evaluator.cc | 12 +++++- .../seb_threecomponent_evaluator.hh | 3 +- .../land_cover/seb_twocomponent_evaluator.cc | 11 ++++- .../land_cover/seb_twocomponent_evaluator.hh | 3 +- 14 files changed, 119 insertions(+), 33 deletions(-) diff --git a/src/pks/surface_balance/constitutive_relations/land_cover/LandCover.cc b/src/pks/surface_balance/constitutive_relations/land_cover/LandCover.cc index 2b60a49f8..2cfb44348 100644 --- a/src/pks/surface_balance/constitutive_relations/land_cover/LandCover.cc +++ b/src/pks/surface_balance/constitutive_relations/land_cover/LandCover.cc @@ -40,6 +40,7 @@ LandCover::LandCover(Teuchos::ParameterList& plist) water_transition_depth(plist.get("water transition depth [m]", NAN)), dessicated_zone_thickness(plist.get("dessicated zone thickness [m]", NAN)), clapp_horn_b(plist.get("Clapp and Hornberger b [-]", NAN)), + rs_method(plist.get("Soil resistance method", "sakagucki_zeng")), roughness_ground(plist.get("roughness length of bare ground [m]", NAN)), roughness_snow(plist.get("roughness length of snow [m]", NAN)), mannings_n(plist.get("Manning's n [?]", NAN)) diff --git a/src/pks/surface_balance/constitutive_relations/land_cover/LandCover.hh b/src/pks/surface_balance/constitutive_relations/land_cover/LandCover.hh index d05d97cea..42d4c1502 100644 --- a/src/pks/surface_balance/constitutive_relations/land_cover/LandCover.hh +++ b/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. }; diff --git a/src/pks/surface_balance/constitutive_relations/land_cover/evaporation_downregulation_evaluator.cc b/src/pks/surface_balance/constitutive_relations/land_cover/evaporation_downregulation_evaluator.cc index 38dfce8fa..7b7399f12 100644 --- a/src/pks/surface_balance/constitutive_relations/land_cover/evaporation_downregulation_evaluator.cc +++ b/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(sat_gas_key_, tag).ViewComponent("cell", false); + const Epetra_MultiVector& sat_liq = + *S.Get(sat_liq_key_, tag).ViewComponent("cell", false); const Epetra_MultiVector& poro = *S.Get(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(sat_gas_key_, tag).ViewComponent("cell", false); + const Epetra_MultiVector& sat_liq = + *S.Get(sat_liq_key_, tag).ViewComponent("cell", false); const Epetra_MultiVector& poro = *S.Get(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(sat_gas_key_, tag) .SetMesh(S.GetMesh(domain_sub_)) ->AddComponent("cell", AmanziMesh::Entity_kind::CELL, 1); + S.Require(sat_liq_key_, tag) + .SetMesh(S.GetMesh(domain_sub_)) + ->AddComponent("cell", AmanziMesh::Entity_kind::CELL, 1); S.Require(pot_evap_key_, tag) .SetMesh(S.GetMesh(domain_surf_)) ->AddComponent("cell", AmanziMesh::Entity_kind::CELL, 1); diff --git a/src/pks/surface_balance/constitutive_relations/land_cover/evaporation_downregulation_evaluator.hh b/src/pks/surface_balance/constitutive_relations/land_cover/evaporation_downregulation_evaluator.hh index d2fde3c68..457b690f0 100644 --- a/src/pks/surface_balance/constitutive_relations/land_cover/evaporation_downregulation_evaluator.hh +++ b/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_; diff --git a/src/pks/surface_balance/constitutive_relations/land_cover/evaporation_downregulation_model.cc b/src/pks/surface_balance/constitutive_relations/land_cover/evaporation_downregulation_model.cc index f0ac58321..8a3995ca9 100644 --- a/src/pks/surface_balance/constitutive_relations/land_cover/evaporation_downregulation_model.cc +++ b/src/pks/surface_balance/constitutive_relations/land_cover/evaporation_downregulation_model.cc @@ -22,6 +22,7 @@ EvaporationDownregulationModel::EvaporationDownregulationModel(Teuchos::Paramete { dess_dz_ = plist.get("dessicated zone thickness [m]", 0.1); Clapp_Horn_b_ = plist.get("Clapp and Hornberger b of surface soil [-]", 1.0); + rs_method_ = plist.get("Soil resistance method", "sakagucki_zeng"); } // Constructor from LandCover @@ -29,26 +30,48 @@ EvaporationDownregulationModel::EvaporationDownregulationModel(const LandCover& { 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.; } @@ -56,9 +79,15 @@ EvaporationDownregulationModel::DEvaporationDPorosity(double sg, double poro, do 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 diff --git a/src/pks/surface_balance/constitutive_relations/land_cover/evaporation_downregulation_model.hh b/src/pks/surface_balance/constitutive_relations/land_cover/evaporation_downregulation_model.hh index 462e6d7d2..f241a4095 100644 --- a/src/pks/surface_balance/constitutive_relations/land_cover/evaporation_downregulation_model.hh +++ b/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,11 +46,13 @@ 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); @@ -55,6 +60,7 @@ class EvaporationDownregulationModel { protected: double dess_dz_; double Clapp_Horn_b_; + std::string rs_method_; }; } // namespace Relations diff --git a/src/pks/surface_balance/constitutive_relations/land_cover/radiation_balance_evaluator.cc b/src/pks/surface_balance/constitutive_relations/land_cover/radiation_balance_evaluator.cc index 1c672821c..f5bb21ad0 100644 --- a/src/pks/surface_balance/constitutive_relations/land_cover/radiation_balance_evaluator.cc +++ b/src/pks/surface_balance/constitutive_relations/land_cover/radiation_balance_evaluator.cc @@ -167,10 +167,9 @@ RadiationBalanceEvaluator::Evaluate_(const State& S, const std::vector 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, diff --git a/src/pks/surface_balance/constitutive_relations/land_cover/seb_physics_funcs.hh b/src/pks/surface_balance/constitutive_relations/land_cover/seb_physics_funcs.hh index fef2bc3d4..df8dc2647 100644 --- a/src/pks/surface_balance/constitutive_relations/land_cover/seb_physics_funcs.hh +++ b/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. diff --git a/src/pks/surface_balance/constitutive_relations/land_cover/seb_threecomponent_evaluator.cc b/src/pks/surface_balance/constitutive_relations/land_cover/seb_threecomponent_evaluator.cc index 10c35e22d..d81e55a1c 100644 --- a/src/pks/surface_balance/constitutive_relations/land_cover/seb_threecomponent_evaluator.cc +++ b/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(sat_gas_key_, tag).ViewComponent("cell", false); + const auto& sat_liq = *S.Get(sat_liq_key_, tag).ViewComponent("cell", false); const auto& poro = *S.Get(poro_key_, tag).ViewComponent("cell", false); const auto& ss_pres = *S.Get(ss_pres_key_, tag).ViewComponent("cell", false); @@ -284,11 +287,13 @@ SEBThreeComponentEvaluator::Evaluate_(const State& S, const std::vector(sat_gas_key_, tag).ViewComponent("cell", false); + const auto& sat_liq = *S.Get(sat_liq_key_, tag).ViewComponent("cell", false); const auto& poro = *S.Get(poro_key_, tag).ViewComponent("cell", false); const auto& ss_pres = *S.Get(ss_pres_key_, tag).ViewComponent("cell", false); @@ -283,11 +286,13 @@ SEBTwoComponentEvaluator::Evaluate_(const State& S, const std::vectorSetGhosted()->AddComponent("cell", AmanziMesh::CELL, 1); diff --git a/src/pks/surface_balance/constitutive_relations/land_cover/seb_twocomponent_evaluator.hh b/src/pks/surface_balance/constitutive_relations/land_cover/seb_twocomponent_evaluator.hh index deae17fab..bf7c4475d 100644 --- a/src/pks/surface_balance/constitutive_relations/land_cover/seb_twocomponent_evaluator.hh +++ b/src/pks/surface_balance/constitutive_relations/land_cover/seb_twocomponent_evaluator.hh @@ -71,6 +71,7 @@ the ground from the atmosphere. - `"temperature`" **DOMAIN-temperature** [K] surface skin temperature. - `"pressure`" **DOMAIN-pressure** [Pa] surface skin pressure. - `"gas saturation`" **DOMAIN_SS-saturation_gas** [-] subsurface gas saturation + - `"liquid saturation`" **DOMAIN_SS-saturation_liquid** [-] subsurface liquid saturation - `"porosity`" [-] subsurface porosity - `"subsurface pressure`" **DOMAIN_SS-pressure** [Pa] - `"molar density liquid`" **DOMAIN-molar_density_liquid** [mol m^-3] @@ -147,7 +148,7 @@ class SEBTwoComponentEvaluator : public EvaluatorSecondaryMonotypeCV { Key ponded_depth_key_, unfrozen_fraction_key_; Key sg_albedo_key_, sg_emissivity_key_, area_frac_key_; Key surf_temp_key_, surf_pres_key_; - Key sat_gas_key_, poro_key_, ss_pres_key_; + Key sat_gas_key_, sat_liq_key_, poro_key_, ss_pres_key_; Key mol_dens_key_, mass_dens_key_; Key melt_key_, evap_key_;