SNOWPACK 3.7.0
SnLaws Class Reference

#include <Laws_sn.h>

Public Types

enum  EventType { event_none , event_wind }
 Types of events for computing new snow density. More...
 
enum  TempDependence { t_term_arrhenius_critical =11 , t_term_arrhenius =911 , t_term_837 =977 , t_term_stk =999 }
 Defines temperature dependence of snow viscosity. More...
 
enum  ViscosityVersion {
  visc_dflt =111 , visc_cal =333 , visc_ant =335 , visc_897 =555 ,
  visc_837 =977 , visc_stk =999
}
 Defines which snow viscosity version to use. More...
 

Static Public Member Functions

static double conductivity_air (void)
 
THERMAL CONDUCTIVITY OF ICE

Based on master thesis of Tobias Hipp, who used relationships by Ling & Zhang (2004).

Version
11.03
Parameters
TemperatureTemperature (K)
Returns
Thermal conductivity of ice
static double conductivity_ice (const double &Temperature)
 
THERMAL CONDUCTIVITY OF WATER

Based on master thesis of Tobias Hipp, who used relationships by Ling & Zhang (2004).

Version
11.03
Parameters
TemperatureTemperature (K)
Returns
Thermal conductivity of water
static double conductivity_water (const double &Temperature)
 
Snow albedo

Computes the density of new snow. The options for SNOW_ALBEDO are:

  • PARAMETERIZED (default is LEHNING_2):
    • LEHNING_[012] : Statistical models of surface snow albedo based on measurements from the Weissfluhjoch study plot (SWin and SWout, K&Z CM21).
    • SCHMUCKI_* : Edgar Schmucki's statistical models (Jan 2014) based on SWin and SWout measurements at 4 stations: Weissfluhjoch study plot (WFJ, 2540 m asl; K&Z CM21), Davos (DAV, 1594 m asl; K&Z CM21), Napf (NAP, 1404 m asl; K&Z CM21), and Payerne (PAY, 490 m asl; K&Z CM21). The variant (*_GSZ) considers grain size as a parameter, while the *_OGS variant replaces grain size with optical equivalent grain size. The user can choose from two average values, ALL_DATA and CUSTOM, the former being the average of all albedo values obtained from all four stations, the latter being the mean of the albedo averages for each single station WFJ, DAV, and PAY; it seems better suited for stations lying below 1500 m asl in Switzerland.
    • NIED : The Japanese version of LEHNING_2
  • MEASURED: Use measured incoming and reflected shortwave radiation fluxes; limited quality checks will be performed. The chosen parameterization will be computed for comparison.
  • FIXED: Use a fixed albedo by assigning ALBEDO-FIXEDVALUE a value between 0.05 and 0.95.
    Parameters
    i_snow_albedotype of albedo computation ()
    i_albedo_parameterization(see above)
    i_albAverageSchmucki
    i_albedo_fixed_valueto use
    Edatacompute albedo for this element
    TssSnow surface temperature (K)
    Mdata
static double compWindPumpingDisplacement (const SnowStation &Xdata)
 Computes the displacement depth in case of ventilation. More...
 
static double compWindPumpingVelocity (const CurrentMeteo &Mdata, const double &d_pump)
 Computes the wind pumping velocity at the surface. More...
 
static double compWindGradientSnow (const ElementData &Edata, double &v_pump)
 Computes the wind gradient near the snow-cover surface. More...
 
static double compSensibleHeatCoefficient (const CurrentMeteo &Mdata, const SnowStation &Xdata, const double &height_of_meteo_values)
 SENSIBLE HEAT EXCHANGE COEFFICIENT (Surface Energy Exchange) More...
 
static double compLatentHeat_Rh (const std::string soil_evaporation, const CurrentMeteo &Mdata, SnowStation &Xdata, const double &height_of_meteo_values)
 Latent heat flux including consideration of soil (one active element) This method uses the Relative Humidity approach: ql = beta*(eA - eS) Latent heat transfer. eA and eS are the vapor pressures of air and snow, respectively. More...
 
static double compLatentHeat (const std::string soil_evaporation, const CurrentMeteo &Mdata, SnowStation &Xdata, const double &height_of_meteo_values)
 LATENT HEAT EXCHANGE (Surface Energy Exchange) David Gustafsson (david.nosp@m.g@kt.nosp@m.h.se) has introduced a resistance approach for latent heat exchange from bare soil as alternative to the relative hyumidity (RH) approach implemented in Snowpack.c (line 473):
An additional resistance, dependent on the relative saturation of the top soil layer, is used to reduce the heat exchange coefficient in the case of evaporation: c = 1/(Ra + Rsoil), where Ra = 1/c as computed above, and Rsoil = 50 [s/m] * field_capacity_soil / theta_soil.
A key SNOWPACK_ADVANCED::soil_evaporation is defined to select method. The resistance formulation originates from van den Hurk et al.(2000) "Offline validation of the ERA40 surface scheme": ECMWF Tech.Memo 295.
A difference from the RH method is that the surface vapour pressure is always assumed to be at saturation at the surface. Thus, some unrealistic effects of the RH method in present form are avoided -> the RH approach tend to estimate negative vapour gradient between the surface and the atmosphere, causing large condensation, even if the top soil layer is saturated, and even if the soil surface is warmer than the atmosphere! If a RH method should work in a discretized model, it is important to consider the difference between vapour pressure at the surface and the average of the top soil layer.
The soil resistance is only used for bare soil layers, when TSS >= 0C and eSurf >= eAtm. More...
 
static double compSoilThermalConductivity (const ElementData &Edata, const double &dvdz, const std::string &soil_thermal_conductivity)
 Heat conduction in soil. The formulation is based on curve fitting, the frozen soil data from Kersten in "Geotechnical Engeneering for Cold Regions" article by Harlan and Nixon, the water influence deduced from deVries and Afgan in "Heat and Mass Transfer in the Biosphere". More...
 
static double soilVaporDiffusivity (const ElementData &Edata)
 Water vapor diffusion coefficient in soil. The formulation is based on Saito et al., 2006 "Numerical analysis of coupled water vapor and heat transport in the vadose zone", see eq. [14] and [15]. It is defined as the product of the tortuosity factor (-) as defined by Millington and Quirck (1961), the air-filled porosity (m3 m-3) and the diffusivity of water vapor in air (m2 s-1). More...
 
static double compEnhanceWaterVaporTransportSoil (const ElementData &Edata, const double &clay_fraction)
 Computes the enhancement factor for water vapor transport in soil. Derived from Cass et al., 1984 "Enhancement of thermal water vapor diffusion in soil", see eq. [19]. Describe the increase in thermal vapor flux as a result of liquid islands and increased temperature gradients in the air phase. More...
 
static double compSoilThermalVaporConductivity (const ElementData &Edata_bot, const ElementData &Edata_top, const double &Te_bot, const double &Te_top, const double &clay_fraction)
 Computes the soil THERMAL vapor hydraulic conductivity. Requires the use of RE to determine pressure head (Edata.h) and saturated water content (theta_s). The THERMAL vapor hydraulic conductivy formulation is based on Saito et al., 2006 "Numerical analysis of coupled water, vapor, and heat transport in the vadose zone", see eq. [13]. It is used to determine the flux density of water vapor in soil due to THERMAL gradient: q_vT = -Kvapor_T*gradT. The enhancement factor is used to describe the increase in the thermal vapor flux as a result of liquid islands and increased temperature gradients in the air phase (Philip and de Vries, 1957) The relative humidity is calculated from the pressure head (h), using a thermodynamic relationship between liquid water and water vapour in soil pores (Philip and de Vries, 1957) More...
 
static double compSoilIsothermalVaporConductivity (const ElementData &Edata_bot, const ElementData &Edata_top, const double &Te_bot, const double &Te_top, const double &T_node)
 Computes the soil ISOTHERMAL vapor hydraulic conductivity. The ISOTHERMAL vapor hydraulic conductivy formulation is based on Saito et al., 2006 "Numerical analysis of coupled water, vapor, and heat transport in the vadose zone", see eq. [12]. It is used to determine the flux density of water vapor in soil due to MOISTURE gradient: q_vh = -Kvapor_h*gradH. More...
 
static double compSnowThermalConductivity (const ElementData &Edata, const double &dvdz, const bool &show_warnings=true)
 Heat conduction in snow Actual version: k = C1*[C2 + C3 + C4 + C5]
Adams/Sato model. The following piece of code was programmed by Perry on a warm June afternoon in Birmensdorf. H.U. Gubler was viciously attacking Ammann, saying nasty things about our beloved chef, (even if his hair is a bit funny), in the DAVOSER ZEITUNG. Perry was dreaming about San Diego and .... and Michael was dreamy about rock climbing in Italy with funky, spunky Italian girls, Borolo Bob was busy with his wife on his way to Mallorca – to watch sexy TOPLESS bikini-slip clad GERMAN secretaries on the BEACH. As usual, if this code does not work it is his FAULT .... (Hey, Bob, when are we going to meet ED ADAMS? And remember PLEDGE 7 of Metamorphism.c. The theory behind this piece of code can be found in Bob's ROUGH DRAFT on MICROSTRUCTURE. More...
 
static double compEnhanceWaterVaporTransportSnow (const SnowStation &Xdata, const size_t &i_e)
 Computes the enhancement factor for water vapor transport in wet snow. More...
 
static double compLWRadCoefficient (const double &t_snow, const double &t_atm, const double &e_atm)
 LONGWAVE RADIATION COEFFICIENT This routine might look a bit unusual: Radiation is treated as a CONVECTIVE boundary condition, similar to the sensible and latent heat exchanges. The exchange coefficient, however, is not a constant, dependent on say the wind velocity, rather it is dependent on the temperature. This routine computes the "pseudo" convective heat exchange coefficient for radiation. More...
 
static double parameterizedSnowAlbedo (const std::string &i_albedo, const std::string &i_albedo_parameterization, const std::string &i_albAverageSchmucki, const double &i_albNIED_av, const double &i_hn_albedo_fixedValue, const ElementData &Edata, const double &Tss, const CurrentMeteo &Mdata, const bool &ageAlbedo=true)
 
static void compShortWaveAbsorption (const std::string &i_sw_absorption_scheme, SnowStation &Xdata, const double &I0)
 Helen LeVesconte's Solution to short wave absorption by the snowpack NOTE on fudge_bohren (fb): Larger values increase extinction --> Energy stays on top. More...
 
static void compAdvectiveHeat (SnowStation &Xdata, const double &advective_heat, const double &depth_begin, const double &depth_end)
 Advective Heat Flux injection (mimicking heat advection by infiltrating water?) More...
 
static double newSnowDensityPara (const std::string &i_hn_model, double TA, double TSS, double RH, double VW, double HH)
 Parameterized new-snow density. More...
 
static double newSnowDensityEvent (const std::string &variant, const SnLaws::EventType &i_event, const CurrentMeteo &Mdata)
 Event driven new-snow density. More...
 
static double newSnowDensityHendrikx (const double ta, const double tss, const double rh, const double vw)
 Jordy Hendrikx' new snow density parameterization for strong winds (> 2.9 m s-1) More...
 
New snow density

Computes the density of new snow. The options for HN_DENSITY are:

  • PARAMETERIZED (default is LEHNING_NEW):
    • ZWART: Costijn Zwart's model (elaborated 2006; in use since 4 Dec 2007
    • LEHNING_NEW: Improved model by M. Lehning, incl. ad-hoc wind & temperature effects (used until 06/07)
    • LEHNING_OLD: First model by M. Lehning
      Note
      {models by M. Lehning can be augmented with a parameterization for winds > 2.9 m s-1 worked out by J. Hendrikx => set jordy_new_snow in Laws_sn.cc}
    • BELLAIRE: Sascha Bellaire's model (elaborated 2007; used summer/fall 2007)
    • PAHAUT: Edmond Pahaut's model, introduced Sep 1995 in CROCUS by G. Giraud
  • EVENT: Driven by event type, that is,
    • event_wind: Implemented 2009 by Christine Groot Zwaaftink for Antarctic variant
  • MEASURED: Use measured new snow density read from meteo input -Note: Set HN_DENSITY_FIXEDVALUE to 1. to use surface snow density as a "measured" value in case of missing values
  • FIXED: Use a fixed new snow density by assigning HN_DENSITY_FIXEDVALUE a value (default: 100 kg m-3, at least min_hn_density)
    Parameters
    i_hn_densitytype of density computation
    i_hn_density_parameterizationto use
    i_hn_density_fixedValueto use
    MdataMeteorological input
    XdataSnow cover data
    tssSnow surface temperature (K)
    variantwhich is currently used
static double compNewSnowDensity (const std::string &i_hn_density, const std::string &i_hn_density_parameterization, const double &i_hn_density_fixedValue, const CurrentMeteo &Mdata, const SnowStation &Xdata, const double &tss, const std::string &variant)
 
static double NewSnowViscosityLehning (const ElementData &Edata)
 NEW SNOW VISCOSITY (dendritic snow, i.e., dd > 0.)
Actual version : ml_lw_VS_Lehning from r7.7
This is Michael's viscosity routine, which is not a function of micro-structure, but which is nonetheless pretty important because it is numerically STABLE and does predict decent settling rates, sometimes a bit too high. This routine was (is) used for NEW or WET snow. More...
 
static double snowViscosityTemperatureTerm (const double &Te)
 Computes the temperature term of viscosity The modifications for POLAR variant are described in: Steger CR, Reijmer CH, van den Broeke MR, Wever N, Forster RR, Koenig LS, Kuipers Munneke P, Lehning M, Lhermitte S, Ligtenberg SRM, Miège C and Noël BPY (2017) Firn Meltwater Retention on the Greenland Ice Sheet: A Model Comparison. Front. Earth Sci. 5:3. doi: 10.3389/feart.2017.00003: "To improve the agreement with observations, the tunable factors in the snow viscosity scheme (Groot Zwaaftink et al., 2013) for the activation energy of snow Qs and the critical exponent β are set to 16,080 J mol−1 and 0.3, respectively.". More...
 
static double compLoadingRateStress (const std::string &variant, ElementData &Edata, const mio::Date &date)
 Computes the additional stress due to loading rate. More...
 
static double loadingRateStressDEFAULT (ElementData &Edata, const mio::Date &date)
 
static double loadingRateStressCALIBRATION (ElementData &Edata, const mio::Date &date)
 
static double snowViscosityFudgeDEFAULT (const ElementData &Edata)
 Determines the fudge factor for viscosity
This fudge factor takes into account bond-ice imperfections and the effect of liquid water. More...
 
static double snowViscosityFudgeCALIBRATION (const ElementData &Edata, const mio::Date &date)
 To calibrate the fudge factor for viscosity. More...
 

Static Public Attributes

static const double field_capacity_soil = 0.15
 Field capacity of the soil (1). Above this levels, water begins to drain. More...
 
static const bool wind_pump = true
 Switch on or off wind pumping in snow. More...
 
static const bool wind_pump_soil = true
 Switch on or off wind pumping in soil. More...
 
static double sn_dt
 
static std::string current_variant = "DEFAULT"
 
static size_t swa_nBands = 5
 Snow extinction coefficients for absorbtion of SW radiation (swa) More...
 
static std::vector< double > swa_k
 mean extinction coefficient for pure ice per band More...
 
static std::vector< double > swa_pc
 fraction of sun power spectrum per band More...
 
static std::vector< double > swa_fb
 fudge_bohren More...
 
static const double rsoilmin = 50.0
 Minimum soil surface resistance, 50 sm-1 (van den Hurk et al, 2000) More...
 
static const double relsatmin = 0.05
 Minimum relative saturation of top soil layer. More...
 
static const double alpha_por_tor_soil = 0.05
 Ratio of porosity to tortuosity for Soil. More...
 
static const double pore_length_soil = 0.01
 Pore length for surface soil for ventilation (m) More...
 
WIND PUMPING
Note
soil parameters
static const double smallest_viscosity = 1.0e6
 Defines the smallest allowable viscosity (Pa s) that a viscosity law will return
Value is DAMM SMALL – smaller values than this are pretty unrealistic. More...
 
static const bool jordy_new_snow = false
 To use J. Hendrikx's parameterization for wind speeds > 2.9 m s-1. More...
 
static const double wind_ext_coef = 0.1
 Used to decribe advective flux attenuation. More...
 
static const double displacement_coef = 0.7
 Integral of snow density corresponding to the maximal displacement depth d_pump (m) More...
 
static const double alpha_por_tor = 0.07
 Ratio of porosity to tortuosity. More...
 
THERMAL CONDUCTIVITY

Defines the constants and parameters for computing snow and soil thermal conductivity

static const double montana_c_fudge = 0.13
 Factor controlling ice to ice conduction. More...
 
static const double montana_vapor_fudge = 2.5
 Factor controlling increase in water vapor transport and thus energy transport in wet snow. More...
 
static const double montana_v_water_fudge = 20.
 Defines the formerly used MONTANA fudge factor for old wet snow settling (up to snowpack_r8.0). More...
 

Event driven density of new snow

Note
{
These static variables are only defined below, if you want to change them for your purposes you need to do so in the function SnLaws::setStaticData where these parameters are set according to the VARIANT used }
  • event : type of event, currently only event_wind implemented
  • event_wind_lowlim : lower limit for event_wind
  • event_wind_highlim : upper limit for event_wind
static double min_hn_density = 30.
 
static double max_hn_density = 250.0
 
static double event_wind_lowlim = 0.0
 
static double event_wind_highlim = 0.0
 
static EventType event = SnLaws::event_none
 
static const bool __init = SnLaws::setStaticData("DEFAULT", "BUCKET")
 
static bool setStaticData (const std::string &variant, const std::string &watertransportmodel)
 This function is used to give default values to a bunch of static members of SnLaws it is called when the helper variable __init is initialized (compile-time) and everytime the user changes the variant and invokes either compAlbedo(..) or compSnowViscosity(...) More...
 

SNOW VISCOSITY

Computes snow viscosity according to the following options:

  • DEFAULT : default model (corresponding to last calibration)
  • KOJIMA : according to parameterization by Kojima
  • CALIBRATION : used for calibration purposes
    Parameters
    variantsee VARIANT
    i_viscosity_modelto use
    Edata
    datecurrent
static TempDependence t_term = SnLaws::t_term_arrhenius_critical
 
static ViscosityVersion visc = SnLaws::visc_cal
 
static double visc_ice_fudge = 1.
 
static double visc_sp_fudge = 1.
 
static double visc_water_fudge = 0.
 
static bool setfix = false
 
static double SnowViscosityMSU (const ElementData &Edata)
 MONTANA snow viscosity (non-dendritic snow, i.e. if dd=0.); From bb_lw_VS_Montana() (Laws.c r7.7). Bob Brown's MICRO-STRUCTURE law. Clearly the BEST law we have right now but also the most UNSTABLE: note that the viscosity is not only a function of the grain dimensions, but also a function of the overburden stress. A series of equations that collectively give the viscosity Vis = S/eDot. The microstructure parameters rb and rg are obtained through Edata pointer and are in mm. The secondary microstructure parameters L and rc are also calculated in mm. This means that the dimensions of rg, rb, L, & rc are in mm since they show up in the following equations as ratios to give dimensionless numbers. NOTE: The theory is NOT valid for NEW SNOW or WET SNOW. However, we try it for wet snow, since it makes sense for wet snow and is related to the wet snow bond growth (Pressure Sintering) If the snow is new, however, then use LEHNING's law, or any other law ... No, do not use any other law, they don't work. More...
 
static double compSnowViscosity (const std::string &variant, const std::string &i_viscosity_model, const std::string &i_watertransport_model, ElementData &Edata, const mio::Date &date)
 
static double snowViscosityDEFAULT (ElementData &Edata)
 SNOW VISCOSITY (all types of snow) More...
 
static double snowViscosityKOJIMA (const ElementData &Edata)
 SNOW VISCOSITY according to formulation by Kojima. More...
 
static double snowViscosityCALIBRATION (ElementData &Edata, const mio::Date &date)
 Calibrate snow viscosity NOTE This is the test or playground version for calibrating settling. More...
 
static double AirEmissivity (mio::MeteoData &md, const std::string &variant)
 Compute the air emissivity It relies on mio::ILWR_parametrized to get a parametrized ILWR if no measured ILWR is available. This in turn relies either on Brutsaert (clear sky) or Omstedt (cloudiness) or Crawford(cloudiness from ISWR). The air temperature and relative humidity should be available for meanigful results. In any case, the returned air emissivity will be between min_air_emissivity and 1, min_air_emissivity depending on the variant. More...
 
static double AirEmissivity (const double &ilwr, const double &ta, const std::string &variant)
 Compute the air emissivity In any case, the returned air emissivity will be between min_air_emissivity and 1, min_air_emissivity depending on the variant. It returns min_air_emissivity if ta==nodata or ilwr==nodata. More...
 
static double ArrheniusLaw (const double ActEnergy, const double T, const double T_ref)
 Computes an Arrhenius-type temperature dependency. More...
 

Member Enumeration Documentation

◆ EventType

Types of events for computing new snow density.

Enumerator
event_none 
event_wind 

Wind driven deposition of snow (Antarctica)

◆ TempDependence

Defines temperature dependence of snow viscosity.

Enumerator
t_term_arrhenius_critical 

Arrhenius type multiplied by critical function near melting point.

t_term_arrhenius 

pure Arrhenius type w/ excitation energy of ice

t_term_837 

as of revision 837 (from r243)

t_term_stk 

calibration 2009 by Walter Steinkogler (MSc thesis)

◆ ViscosityVersion

Defines which snow viscosity version to use.

Enumerator
visc_dflt 

actual version: revision 961, June 2011 by Fierz

visc_cal 

version under test

visc_ant 

if Antarctica needs adaptation of _new

visc_897 

calibration fall 2010 by Fierz

visc_837 

as of revision 837 (deprecated)

visc_stk 

calibration 2009 by Walter Steinkogler (MSc thesis)

Member Function Documentation

◆ AirEmissivity() [1/2]

double SnLaws::AirEmissivity ( const double &  ilwr,
const double &  ta,
const std::string &  variant 
)
static

Compute the air emissivity In any case, the returned air emissivity will be between min_air_emissivity and 1, min_air_emissivity depending on the variant. It returns min_air_emissivity if ta==nodata or ilwr==nodata.

Parameters
ilwrincoming long wave radiation (W/m^2)
taair temperature (K)
variantvariant to use for the minimum allowed air emissivity
Note
observed minimum air emissivities:
  • default: 0.55 (from 1993 data at Weissfluhjoch)
  • Antarctica: 0.31 (from 2006/2007 data of Dome C)
Returns
air emissivity in range [MIN_AIR_EMISSIVITY,1.] (1)

◆ AirEmissivity() [2/2]

double SnLaws::AirEmissivity ( mio::MeteoData &  md,
const std::string &  variant 
)
static

Compute the air emissivity It relies on mio::ILWR_parametrized to get a parametrized ILWR if no measured ILWR is available. This in turn relies either on Brutsaert (clear sky) or Omstedt (cloudiness) or Crawford(cloudiness from ISWR). The air temperature and relative humidity should be available for meanigful results. In any case, the returned air emissivity will be between min_air_emissivity and 1, min_air_emissivity depending on the variant.

Parameters
mdmeteorological conditions
variantvariant to use for the minimum allowed air emissivity
Note
observed minimum air emissivities:
  • default: 0.55 (from 1993 data at Weissfluhjoch)
  • Antarctica: 0.31 (from 2006/2007 data of Dome C)
Returns
air emissivity in range [MIN_AIR_EMISSIVITY,1.] (1)

◆ ArrheniusLaw()

double SnLaws::ArrheniusLaw ( const double  ActEnergy,
const double  T,
const double  T_ref 
)
static

Computes an Arrhenius-type temperature dependency.

Version
9.11
Parameters
ActEnergy(J mol-1)
Tsnow temperature (K)
T_refa reference temperature (K)

◆ compAdvectiveHeat()

void SnLaws::compAdvectiveHeat ( SnowStation Xdata,
const double &  advective_heat,
const double &  depth_begin,
const double &  depth_end 
)
static

Advective Heat Flux injection (mimicking heat advection by infiltrating water?)

Parameters
XdataThe advective heat flux will be added to Xdata.sw_abs
advective_heatheat flux (positive or negative) to add (W m-3) //NO_HACK: input as W/m2 seems more logical but wrong! It is a source term!
depth_begindepth where to begin injecting the heat flux (in m from the soil surface)
depth_enddepth where to stop injecting the heat flux (in m from the soil surface)

◆ compEnhanceWaterVaporTransportSnow()

double SnLaws::compEnhanceWaterVaporTransportSnow ( const SnowStation Xdata,
const size_t &  i_e 
)
static

Computes the enhancement factor for water vapor transport in wet snow.

Version
9Y.mm
Parameters
Xdata
i_eElement index to start from
Returns
Enhancement factor

◆ compEnhanceWaterVaporTransportSoil()

double SnLaws::compEnhanceWaterVaporTransportSoil ( const ElementData Edata,
const double &  clay_fraction 
)
static

Computes the enhancement factor for water vapor transport in soil. Derived from Cass et al., 1984 "Enhancement of thermal water vapor diffusion in soil", see eq. [19]. Describe the increase in thermal vapor flux as a result of liquid islands and increased temperature gradients in the air phase.

Author
Margaux Couttet
Parameters
Edataelement data
clay_fractionfraction of clay in the soil
Returns
Enhancement factor (-)

◆ compLatentHeat()

double SnLaws::compLatentHeat ( const std::string  soil_evaporation,
const CurrentMeteo Mdata,
SnowStation Xdata,
const double &  height_of_meteo_values 
)
static

LATENT HEAT EXCHANGE (Surface Energy Exchange) David Gustafsson (david.nosp@m.g@kt.nosp@m.h.se) has introduced a resistance approach for latent heat exchange from bare soil as alternative to the relative hyumidity (RH) approach implemented in Snowpack.c (line 473):
An additional resistance, dependent on the relative saturation of the top soil layer, is used to reduce the heat exchange coefficient in the case of evaporation: c = 1/(Ra + Rsoil), where Ra = 1/c as computed above, and Rsoil = 50 [s/m] * field_capacity_soil / theta_soil.
A key SNOWPACK_ADVANCED::soil_evaporation is defined to select method. The resistance formulation originates from van den Hurk et al.(2000) "Offline validation of the ERA40 surface scheme": ECMWF Tech.Memo 295.
A difference from the RH method is that the surface vapour pressure is always assumed to be at saturation at the surface. Thus, some unrealistic effects of the RH method in present form are avoided -> the RH approach tend to estimate negative vapour gradient between the surface and the atmosphere, causing large condensation, even if the top soil layer is saturated, and even if the soil surface is warmer than the atmosphere! If a RH method should work in a discretized model, it is important to consider the difference between vapour pressure at the surface and the average of the top soil layer.
The soil resistance is only used for bare soil layers, when TSS >= 0C and eSurf >= eAtm.

Parameters
[in]soil_evaporationThe evaporation method to be used
[in]Mdata
[in]Xdata
[in]height_of_meteo_valuesHeight at which meteo parameters are measured
Returns
Latent heat flux (W m-2)

◆ compLatentHeat_Rh()

double SnLaws::compLatentHeat_Rh ( const std::string  soil_evaporation,
const CurrentMeteo Mdata,
SnowStation Xdata,
const double &  height_of_meteo_values 
)
static

Latent heat flux including consideration of soil (one active element) This method uses the Relative Humidity approach: ql = beta*(eA - eS) Latent heat transfer. eA and eS are the vapor pressures of air and snow, respectively.

Version
9Y.mm
Parameters
soil_evaporationThe evaporation method to be used
Mdata
Xdata
height_of_meteo_valuesHeight at which meteo parameters are measured
Returns
Latent heat flux (W m-2)

◆ compLoadingRateStress()

double SnLaws::compLoadingRateStress ( const std::string &  i_viscosity_model,
ElementData Edata,
const mio::Date &  date 
)
static

Computes the additional stress due to loading rate.

Version
11.06
Parameters
i_viscosity_modelsee compSnowViscosity()
Edata
datecurrent
Returns
Initial stress (Pa); note it is a negative value!

◆ compLWRadCoefficient()

double SnLaws::compLWRadCoefficient ( const double &  t_snow,
const double &  t_atm,
const double &  e_atm 
)
static

LONGWAVE RADIATION COEFFICIENT This routine might look a bit unusual: Radiation is treated as a CONVECTIVE boundary condition, similar to the sensible and latent heat exchanges. The exchange coefficient, however, is not a constant, dependent on say the wind velocity, rather it is dependent on the temperature. This routine computes the "pseudo" convective heat exchange coefficient for radiation.

Version
9Y.mm
Parameters
t_snowSnow surface temperature (K)
t_atmTemperature of the atmosphere, i.e., air (K)
e_atmEmissivity of the atmosphere (1)
Returns
LW radiation coefficient (W m-2 K-1)

◆ compNewSnowDensity()

double SnLaws::compNewSnowDensity ( const std::string &  i_hn_density,
const std::string &  i_hn_density_parameterization,
const double &  i_hn_density_fixedValue,
const CurrentMeteo Mdata,
const SnowStation Xdata,
const double &  tss,
const std::string &  variant 
)
static

◆ compSensibleHeatCoefficient()

double SnLaws::compSensibleHeatCoefficient ( const CurrentMeteo Mdata,
const SnowStation Xdata,
const double &  height_of_meteo_values 
)
static

SENSIBLE HEAT EXCHANGE COEFFICIENT (Surface Energy Exchange)

Version
9Y.mm
Parameters
Mdata
Xdata
height_of_meteo_valuesHeight at which meteo parameters are measured
Returns
Exchange coefficient for sensible heat (1)

◆ compShortWaveAbsorption()

void SnLaws::compShortWaveAbsorption ( const std::string &  i_sw_absorption_scheme,
SnowStation Xdata,
const double &  I0 
)
static

Helen LeVesconte's Solution to short wave absorption by the snowpack NOTE on fudge_bohren (fb): Larger values increase extinction --> Energy stays on top.

Parameters
i_sw_absorption_schemeuse multi band or single band approach
Xdata
I0net shortwave radiation (W m-2)

◆ compSnowThermalConductivity()

double SnLaws::compSnowThermalConductivity ( const ElementData Edata,
const double &  dvdz,
const bool &  show_warnings = true 
)
static

Heat conduction in snow Actual version: k = C1*[C2 + C3 + C4 + C5]
Adams/Sato model. The following piece of code was programmed by Perry on a warm June afternoon in Birmensdorf. H.U. Gubler was viciously attacking Ammann, saying nasty things about our beloved chef, (even if his hair is a bit funny), in the DAVOSER ZEITUNG. Perry was dreaming about San Diego and .... and Michael was dreamy about rock climbing in Italy with funky, spunky Italian girls, Borolo Bob was busy with his wife on his way to Mallorca – to watch sexy TOPLESS bikini-slip clad GERMAN secretaries on the BEACH. As usual, if this code does not work it is his FAULT .... (Hey, Bob, when are we going to meet ED ADAMS? And remember PLEDGE 7 of Metamorphism.c. The theory behind this piece of code can be found in Bob's ROUGH DRAFT on MICROSTRUCTURE.

Much much later, Michael re-programmed this code and introduced the effect of liquid water. He extended the Adams/Sato model and hoped that he would not loose contact with Keegan, even he was not going to fulfill the explicit wishes of his Grandma. The model was now being used from Finland to the US. We will also conquer the southern hemis- sphere.

Version
8.10 (??)
Parameters
Edata
dvdzWind velocity gradient (s-1)
show_warningsIf warnings are produced, show them (default=true)
Returns
Thermal conductivity of snow (W K-1 m-1)

◆ compSnowViscosity()

double SnLaws::compSnowViscosity ( const std::string &  variant,
const std::string &  i_viscosity_model,
const std::string &  i_watertransport_model,
ElementData Edata,
const mio::Date &  date 
)
static

◆ compSoilIsothermalVaporConductivity()

double SnLaws::compSoilIsothermalVaporConductivity ( const ElementData Edata_bot,
const ElementData Edata_top,
const double &  Te_bot,
const double &  Te_top,
const double &  T_node 
)
static

Computes the soil ISOTHERMAL vapor hydraulic conductivity. The ISOTHERMAL vapor hydraulic conductivy formulation is based on Saito et al., 2006 "Numerical analysis of coupled water, vapor, and heat transport in the vadose zone", see eq. [12]. It is used to determine the flux density of water vapor in soil due to MOISTURE gradient: q_vh = -Kvapor_h*gradH.

Author
Margaux Couttet
Parameters
Edata_botelement data
Edata_topelement data
Te_botlower element temperature (K)
Te_topupper element temperature (K)
T_nodenodal temperature (K)
Returns
isothermal vapor hydraulic conductivity (m s-1)

◆ compSoilThermalConductivity()

double SnLaws::compSoilThermalConductivity ( const ElementData Edata,
const double &  dvdz,
const std::string &  soil_thermal_conductivity 
)
static

Heat conduction in soil. The formulation is based on curve fitting, the frozen soil data from Kersten in "Geotechnical Engeneering for Cold Regions" article by Harlan and Nixon, the water influence deduced from deVries and Afgan in "Heat and Mass Transfer in the Biosphere".

Version
11.03: thermal conductivity made temperature dependent.
12.0: thermal conductivity model is now defined by a key SOIL_THERMAL_CONDUCTIVITY in SNOWPACK_ADVANCED
Parameters
[in]Edata
[in]dvdzWind velocity gradient (s-1)
[in]soil_thermal_conductivityThermal conductivity model to use (either "FITTED" or any other string)
Returns
Soil thermal conductivity (W K-1 m-1)

◆ compSoilThermalVaporConductivity()

double SnLaws::compSoilThermalVaporConductivity ( const ElementData Edata_bot,
const ElementData Edata_top,
const double &  Te_bot,
const double &  Te_top,
const double &  clay_fraction 
)
static

Computes the soil THERMAL vapor hydraulic conductivity. Requires the use of RE to determine pressure head (Edata.h) and saturated water content (theta_s). The THERMAL vapor hydraulic conductivy formulation is based on Saito et al., 2006 "Numerical analysis of coupled water, vapor, and heat transport in the vadose zone", see eq. [13]. It is used to determine the flux density of water vapor in soil due to THERMAL gradient: q_vT = -Kvapor_T*gradT. The enhancement factor is used to describe the increase in the thermal vapor flux as a result of liquid islands and increased temperature gradients in the air phase (Philip and de Vries, 1957) The relative humidity is calculated from the pressure head (h), using a thermodynamic relationship between liquid water and water vapour in soil pores (Philip and de Vries, 1957)

Author
Margaux Couttet
Parameters
Edata_botelement data
Edata_topelement data
Te_botlower element temperature (K)
Te_topupper element temperature (K)
clay_fractionfraction of clay in the soil
Returns
thermal vapor hydraulic conductivity (m2 K-1 s-1)

◆ compWindGradientSnow()

double SnLaws::compWindGradientSnow ( const ElementData Edata,
double &  v_pump 
)
static

Computes the wind gradient near the snow-cover surface.

Parameters
Edata
v_pumpWind velocity at element depth (m s-1)
Returns
Wind pumping velocity gradient (s-1)

◆ compWindPumpingDisplacement()

double SnLaws::compWindPumpingDisplacement ( const SnowStation Xdata)
static

Computes the displacement depth in case of ventilation.

Version
10.01
Parameters
Xdata
Returns
Displacement depth (m)

◆ compWindPumpingVelocity()

double SnLaws::compWindPumpingVelocity ( const CurrentMeteo Mdata,
const double &  d_pump 
)
static

Computes the wind pumping velocity at the surface.

Version
10.01
Parameters
Mdata
d_pumpDisplacement depth (m)
Returns
Wind pumping velocity (m s-1)

◆ conductivity_air()

static double SnLaws::conductivity_air ( void  )
static

◆ conductivity_ice()

double SnLaws::conductivity_ice ( const double &  Temperature)
static

◆ conductivity_water()

double SnLaws::conductivity_water ( const double &  Temperature)
static

◆ loadingRateStressCALIBRATION()

double SnLaws::loadingRateStressCALIBRATION ( ElementData Edata,
const mio::Date &  date 
)
static

◆ loadingRateStressDEFAULT()

double SnLaws::loadingRateStressDEFAULT ( ElementData Edata,
const mio::Date &  date 
)
static

◆ newSnowDensityEvent()

double SnLaws::newSnowDensityEvent ( const std::string &  variant,
const SnLaws::EventType i_event,
const CurrentMeteo Mdata 
)
static

Event driven new-snow density.

Parameters
variantSnowpack variant (such as DEFAULT, POLAR...)
i_event- event_wind: rho = 250.3 kg m-3 @ 4 m s-1; rho = 338 kg m-3 @ 7 m s-1 Antarctica
MdataMeteorological input

◆ newSnowDensityHendrikx()

double SnLaws::newSnowDensityHendrikx ( const double  ta,
const double  tss,
const double  rh,
const double  vw 
)
static

Jordy Hendrikx' new snow density parameterization for strong winds (> 2.9 m s-1)

Note
To be used with Lehning's models only!
Parameters
taAir temperature (degC)
tssSnow surface temperature (degC)
rhRelative air humidity (%)
vwMean wind velocity (m s-1)
Returns
New snow density

◆ newSnowDensityPara()

double SnLaws::newSnowDensityPara ( const std::string &  i_hn_model,
double  TA,
double  TSS,
double  RH,
double  VW,
double  HH 
)
static

Parameterized new-snow density.

Parameters
TAAir temperature (K)
TSSSnow surface temperature (K)
RHRelative air humidity (1)
VWMean wind velocity (m s-1)
HHAltitude a.s.l. (m)
i_hn_modelParameterization to be used
Returns
New snow density (kg m-3)

◆ NewSnowViscosityLehning()

double SnLaws::NewSnowViscosityLehning ( const ElementData Edata)
static

NEW SNOW VISCOSITY (dendritic snow, i.e., dd > 0.)
Actual version : ml_lw_VS_Lehning from r7.7
This is Michael's viscosity routine, which is not a function of micro-structure, but which is nonetheless pretty important because it is numerically STABLE and does predict decent settling rates, sometimes a bit too high. This routine was (is) used for NEW or WET snow.

Version
9Y.mm
Returns
Viscosity of new snow (Pa s)

◆ parameterizedSnowAlbedo()

double SnLaws::parameterizedSnowAlbedo ( const std::string &  i_albedo,
const std::string &  i_albedo_parameterization,
const std::string &  i_albAverageSchmucki,
const double &  i_albNIED_av,
const double &  i_hn_albedo_fixedValue,
const ElementData Edata,
const double &  Tss,
const CurrentMeteo Mdata,
const bool &  ageAlbedo = true 
)
static

◆ setStaticData()

bool SnLaws::setStaticData ( const std::string &  variant,
const std::string &  watertransport_model 
)
static

This function is used to give default values to a bunch of static members of SnLaws it is called when the helper variable __init is initialized (compile-time) and everytime the user changes the variant and invokes either compAlbedo(..) or compSnowViscosity(...)

Returns
always true

◆ snowViscosityCALIBRATION()

double SnLaws::snowViscosityCALIBRATION ( ElementData Edata,
const mio::Date &  date 
)
static

Calibrate snow viscosity NOTE This is the test or playground version for calibrating settling.

Version
10.11
Parameters
Edata
datecurrent
Returns
Snow viscosity (Pa s)

◆ snowViscosityDEFAULT()

double SnLaws::snowViscosityDEFAULT ( ElementData Edata)
static

SNOW VISCOSITY (all types of snow)

Version
11.02
Parameters
Edata
Returns
Snow viscosity (Pa s)

◆ snowViscosityFudgeCALIBRATION()

double SnLaws::snowViscosityFudgeCALIBRATION ( const ElementData Edata,
const mio::Date &  date 
)
static

To calibrate the fudge factor for viscosity.

Note
This is the fudge playground for calibrations used in SnowViscosityCALIBRATION()
Version
10.11
Parameters
Edata
datecurrent date
Returns
Fudge factor for snow viscosity

◆ snowViscosityFudgeDEFAULT()

double SnLaws::snowViscosityFudgeDEFAULT ( const ElementData Edata)
static

Determines the fudge factor for viscosity
This fudge factor takes into account bond-ice imperfections and the effect of liquid water.

Version
11.06
Parameters
Edata
Returns
Fudge factor for snow viscosity

◆ snowViscosityKOJIMA()

double SnLaws::snowViscosityKOJIMA ( const ElementData Edata)
static

SNOW VISCOSITY according to formulation by Kojima.

Version
9.10
Parameters
Edata
Returns
Snow viscosity (Pa s)

◆ SnowViscosityMSU()

double SnLaws::SnowViscosityMSU ( const ElementData Edata)
static

MONTANA snow viscosity (non-dendritic snow, i.e. if dd=0.); From bb_lw_VS_Montana() (Laws.c r7.7). Bob Brown's MICRO-STRUCTURE law. Clearly the BEST law we have right now but also the most UNSTABLE: note that the viscosity is not only a function of the grain dimensions, but also a function of the overburden stress. A series of equations that collectively give the viscosity Vis = S/eDot. The microstructure parameters rb and rg are obtained through Edata pointer and are in mm. The secondary microstructure parameters L and rc are also calculated in mm. This means that the dimensions of rg, rb, L, & rc are in mm since they show up in the following equations as ratios to give dimensionless numbers. NOTE: The theory is NOT valid for NEW SNOW or WET SNOW. However, we try it for wet snow, since it makes sense for wet snow and is related to the wet snow bond growth (Pressure Sintering) If the snow is new, however, then use LEHNING's law, or any other law ... No, do not use any other law, they don't work.

Author
Guess!
Version
Antediluvienne
Parameters
*Edata
Returns
Snow viscosity Montana style

◆ snowViscosityTemperatureTerm()

double SnLaws::snowViscosityTemperatureTerm ( const double &  Te)
static

Computes the temperature term of viscosity The modifications for POLAR variant are described in: Steger CR, Reijmer CH, van den Broeke MR, Wever N, Forster RR, Koenig LS, Kuipers Munneke P, Lehning M, Lhermitte S, Ligtenberg SRM, Miège C and Noël BPY (2017) Firn Meltwater Retention on the Greenland Ice Sheet: A Model Comparison. Front. Earth Sci. 5:3. doi: 10.3389/feart.2017.00003: "To improve the agreement with observations, the tunable factors in the snow viscosity scheme (Groot Zwaaftink et al., 2013) for the activation energy of snow Qs and the critical exponent β are set to 16,080 J mol−1 and 0.3, respectively.".

Version
11.06
Parameters
TeElement temperature (K)
Returns
Temperature term of snow viscosity

◆ soilVaporDiffusivity()

double SnLaws::soilVaporDiffusivity ( const ElementData Edata)
static

Water vapor diffusion coefficient in soil. The formulation is based on Saito et al., 2006 "Numerical analysis of coupled water vapor and heat transport in the vadose zone", see eq. [14] and [15]. It is defined as the product of the tortuosity factor (-) as defined by Millington and Quirck (1961), the air-filled porosity (m3 m-3) and the diffusivity of water vapor in air (m2 s-1).

Author
Margaux Couttet
Parameters
Edataelement data
Returns
vapour diffusivity in soil (m2 s-1)

Member Data Documentation

◆ __init

const bool SnLaws::__init = SnLaws::setStaticData("DEFAULT", "BUCKET")
static

◆ alpha_por_tor

const double SnLaws::alpha_por_tor = 0.07
static

Ratio of porosity to tortuosity.

◆ alpha_por_tor_soil

const double SnLaws::alpha_por_tor_soil = 0.05
static

Ratio of porosity to tortuosity for Soil.

◆ current_variant

std::string SnLaws::current_variant = "DEFAULT"
static

◆ displacement_coef

const double SnLaws::displacement_coef = 0.7
static

Integral of snow density corresponding to the maximal displacement depth d_pump (m)

◆ event

SnLaws::EventType SnLaws::event = SnLaws::event_none
static

◆ event_wind_highlim

double SnLaws::event_wind_highlim = 0.0
static

◆ event_wind_lowlim

double SnLaws::event_wind_lowlim = 0.0
static

◆ field_capacity_soil

const double SnLaws::field_capacity_soil = 0.15
static

Field capacity of the soil (1). Above this levels, water begins to drain.

◆ jordy_new_snow

const bool SnLaws::jordy_new_snow = false
static

To use J. Hendrikx's parameterization for wind speeds > 2.9 m s-1.

◆ max_hn_density

double SnLaws::max_hn_density = 250.0
static

◆ min_hn_density

double SnLaws::min_hn_density = 30.
static

◆ montana_c_fudge

const double SnLaws::montana_c_fudge = 0.13
static

Factor controlling ice to ice conduction.

◆ montana_v_water_fudge

const double SnLaws::montana_v_water_fudge = 20.
static

Defines the formerly used MONTANA fudge factor for old wet snow settling (up to snowpack_r8.0).

◆ montana_vapor_fudge

const double SnLaws::montana_vapor_fudge = 2.5
static

Factor controlling increase in water vapor transport and thus energy transport in wet snow.

◆ pore_length_soil

const double SnLaws::pore_length_soil = 0.01
static

Pore length for surface soil for ventilation (m)

◆ relsatmin

const double SnLaws::relsatmin = 0.05
static

Minimum relative saturation of top soil layer.

◆ rsoilmin

const double SnLaws::rsoilmin = 50.0
static

Minimum soil surface resistance, 50 sm-1 (van den Hurk et al, 2000)

◆ setfix

bool SnLaws::setfix = false
static

◆ smallest_viscosity

const double SnLaws::smallest_viscosity = 1.0e6
static

Defines the smallest allowable viscosity (Pa s) that a viscosity law will return
Value is DAMM SMALL – smaller values than this are pretty unrealistic.

◆ sn_dt

double SnLaws::sn_dt
static

◆ swa_fb

std::vector< double > SnLaws::swa_fb
static

fudge_bohren

◆ swa_k

std::vector< double > SnLaws::swa_k
static

mean extinction coefficient for pure ice per band

◆ swa_nBands

size_t SnLaws::swa_nBands = 5
static

Snow extinction coefficients for absorbtion of SW radiation (swa)

Number of bands used

◆ swa_pc

std::vector< double > SnLaws::swa_pc
static

fraction of sun power spectrum per band

◆ t_term

◆ visc

◆ visc_ice_fudge

double SnLaws::visc_ice_fudge = 1.
static

◆ visc_sp_fudge

double SnLaws::visc_sp_fudge = 1.
static

◆ visc_water_fudge

double SnLaws::visc_water_fudge = 0.
static

◆ wind_ext_coef

const double SnLaws::wind_ext_coef = 0.1
static

Used to decribe advective flux attenuation.

◆ wind_pump

const bool SnLaws::wind_pump = true
static

Switch on or off wind pumping in snow.

◆ wind_pump_soil

const bool SnLaws::wind_pump_soil = true
static

Switch on or off wind pumping in soil.


The documentation for this class was generated from the following files: