LE_to_ET(LE,Tair)
ET_to_LE(ET,Tair)

Convert evaporative water flux from mass (ET=evapotranspiration) to energy (LE=latent heat flux) units, or vice versa.

Arguments

• LE Latent heat flux (W m-2)
• ET Evapotranspiration (kg m-2 s-1)
• Tair Air temperature (deg C)

Details

The conversions are given by:

• $ET = LE/\lambda$
• $LE = \lambda ET$

where $\lambda$ is the latent heat of vaporization (J kg-1) as calculated by latent_heat_vaporization.

Examples

# LE of 200 Wm-2 and air temperature of 25degC
ET = LE_to_ET(200,25)
≈(ET, 8.19e-5, atol =1e-7)

Esat_slope(Tair; Esat_formula, constants) 
Esat_from_Tair(Tair; Esat_formula, constants) 
Esat_from_Tair_deriv(Tair; Esat_formula, constants)

Saturation Vapor Pressure (Esat) and Slope of the Esat Curve

Arguemtns

• Tair: Air temperature (deg C)
• Esat_formula=Val(:Sonntag_1990): Esatformula to be used. Either `Val(:Sonntag1990)(Default),Val(:Alduchov1996, orVal(:Allen1998)`.
• constants=bigleaf_constants(): Dictionary with entry :Pa2kPa

Details

Esat (kPa) is calculated using the Magnus equation:

Esat = a * exp((b * Tair) / (c + Tair)) / 1000 where the coefficients a, b, c take different values depending on the Esat_formula use The default values are from Sonntag 1990 (a=611.2, b=17.62, c=243.12). This versi of the Magnus equation is recommended by the WMO (WMO 2008; p1.4-29). Alternativel parameter values determined by Alduchov & Eskridge 1996 or Allen et al. 1998 can b used (see references). The slope of the Esat curve ($\Delta$) is calculated as the first derivative of the function:

$\Delta = {dEsat \over dTair}$

Value

• Esat_from_Tair: Saturation vapor pressure (kPa)
• Esat_from_Tair_deriv: Slope of the saturation vapor pressure curve (kPa K-1)
• Esat_slope: A tuple of the two above values

References

Sonntag D. 1990: Important new values of the physical constants of 1986, vapor pressure formulations based on the ITS-90 and psychrometric formulae. Zeitschrift fuer Meteorologie 70, 340-344.

World Meteorological Organization 2008: Guide to Meteorological Instruments and Methods of Observation (WMO-No.8). World Meteorological Organization, Geneva. 7th Edition,

Alduchov, O. A. & Eskridge, R. E., 1996: Improved Magnus form approximation of saturation vapor pressure. Journal of Applied Meteorology, 35, 601-609

Allen, R.G., Pereira, L.S., Raes, D., Smith, M., 1998: Crop evapotranspiration - Guidelines for computing crop water requirements - FAO irrigation and drainage paper 56, FAO, Rome.

Esat_from_Tair(20.0)          # Esat in kPa
Esat_from_Tair_deriv(20.0)    # its derivative to temperature in kPa K-1
Esat_slope(20.0)              # both as a tuple

Esat_slope(Tair; Esat_formula, constants) 
Esat_from_Tair(Tair; Esat_formula, constants) 
Esat_from_Tair_deriv(Tair; Esat_formula, constants)

Saturation Vapor Pressure (Esat) and Slope of the Esat Curve

Arguemtns

• Tair: Air temperature (deg C)
• Esat_formula=Val(:Sonntag_1990): Esatformula to be used. Either `Val(:Sonntag1990)(Default),Val(:Alduchov1996, orVal(:Allen1998)`.
• constants=bigleaf_constants(): Dictionary with entry :Pa2kPa

Details

Esat (kPa) is calculated using the Magnus equation:

Esat = a * exp((b * Tair) / (c + Tair)) / 1000 where the coefficients a, b, c take different values depending on the Esat_formula use The default values are from Sonntag 1990 (a=611.2, b=17.62, c=243.12). This versi of the Magnus equation is recommended by the WMO (WMO 2008; p1.4-29). Alternativel parameter values determined by Alduchov & Eskridge 1996 or Allen et al. 1998 can b used (see references). The slope of the Esat curve ($\Delta$) is calculated as the first derivative of the function:

$\Delta = {dEsat \over dTair}$

Value

• Esat_from_Tair: Saturation vapor pressure (kPa)
• Esat_from_Tair_deriv: Slope of the saturation vapor pressure curve (kPa K-1)
• Esat_slope: A tuple of the two above values

References

Sonntag D. 1990: Important new values of the physical constants of 1986, vapor pressure formulations based on the ITS-90 and psychrometric formulae. Zeitschrift fuer Meteorologie 70, 340-344.

World Meteorological Organization 2008: Guide to Meteorological Instruments and Methods of Observation (WMO-No.8). World Meteorological Organization, Geneva. 7th Edition,

Alduchov, O. A. & Eskridge, R. E., 1996: Improved Magnus form approximation of saturation vapor pressure. Journal of Applied Meteorology, 35, 601-609

Allen, R.G., Pereira, L.S., Raes, D., Smith, M., 1998: Crop evapotranspiration - Guidelines for computing crop water requirements - FAO irrigation and drainage paper 56, FAO, Rome.

Esat_from_Tair(20.0)          # Esat in kPa
Esat_from_Tair_deriv(20.0)    # its derivative to temperature in kPa K-1
Esat_slope(20.0)              # both as a tuple

Esat_slope(Tair; Esat_formula, constants) 
Esat_from_Tair(Tair; Esat_formula, constants) 
Esat_from_Tair_deriv(Tair; Esat_formula, constants)

Saturation Vapor Pressure (Esat) and Slope of the Esat Curve

Arguemtns

• Tair: Air temperature (deg C)
• Esat_formula=Val(:Sonntag_1990): Esatformula to be used. Either `Val(:Sonntag1990)(Default),Val(:Alduchov1996, orVal(:Allen1998)`.
• constants=bigleaf_constants(): Dictionary with entry :Pa2kPa

Details

Esat (kPa) is calculated using the Magnus equation:

Esat = a * exp((b * Tair) / (c + Tair)) / 1000 where the coefficients a, b, c take different values depending on the Esat_formula use The default values are from Sonntag 1990 (a=611.2, b=17.62, c=243.12). This versi of the Magnus equation is recommended by the WMO (WMO 2008; p1.4-29). Alternativel parameter values determined by Alduchov & Eskridge 1996 or Allen et al. 1998 can b used (see references). The slope of the Esat curve ($\Delta$) is calculated as the first derivative of the function:

$\Delta = {dEsat \over dTair}$

Value

• Esat_from_Tair: Saturation vapor pressure (kPa)
• Esat_from_Tair_deriv: Slope of the saturation vapor pressure curve (kPa K-1)
• Esat_slope: A tuple of the two above values

References

Sonntag D. 1990: Important new values of the physical constants of 1986, vapor pressure formulations based on the ITS-90 and psychrometric formulae. Zeitschrift fuer Meteorologie 70, 340-344.

World Meteorological Organization 2008: Guide to Meteorological Instruments and Methods of Observation (WMO-No.8). World Meteorological Organization, Geneva. 7th Edition,

Alduchov, O. A. & Eskridge, R. E., 1996: Improved Magnus form approximation of saturation vapor pressure. Journal of Applied Meteorology, 35, 601-609

Allen, R.G., Pereira, L.S., Raes, D., Smith, M., 1998: Crop evapotranspiration - Guidelines for computing crop water requirements - FAO irrigation and drainage paper 56, FAO, Rome.

Esat_from_Tair(20.0)          # Esat in kPa
Esat_from_Tair_deriv(20.0)    # its derivative to temperature in kPa K-1
Esat_slope(20.0)              # both as a tuple

Gb_Choudhury(; leafwidth, LAI, wind_zh, constants)
Gb_Choudhury!(df; leafwidth, LAI, wind_zh, constants)

Estimate the canopy boundary layer conductance for heat transfer according to Choudhury & Monteith 1988.

Arguments

• df : DataFrame where Gb_h is to be added/updated
• leafwidth : Leaf width (m)
• LAI : One-sided leaf area index
• wind_zh : Wind speed at canopy heihgt (m s-1), see wind_profile
• constants=bigleaf_constants()

Value

see compute_Gb!

Details

Boundary layer conductance according to Choudhury & Monteith 1988 is given by:

$Gb_h = LAI \left( 2a/\alpha \sqrt{u(z_h)/w} (1-e^{-\alpha/2})\right)$

where $\alpha$ is modeled as an empirical relation to LAI (McNaughton & van den Hurk 1995):

$\alpha = 4.39 - 3.97 e^{-0.258 \, LAI}$

$w$ is leafwidth and $u(zh)$ is the wind speed at the canopy surface.

It can be approximated from measured wind speed at sensor height zr and a wind extinction coefficient $\alpha$: $u(z_h) = u(z_r) / e^{\alpha(z_r/z_h -1)}$. However, here (if not explicitly given) it is estimated by wind_profile

References

• Choudhury, B. J., Monteith J_L., 1988: A four-layer model for the heat budget of homogeneous land surfaces. Q. J. R. Meteorol. Soc. 114, 373-398.
• McNaughton, K. G., Van den Hurk, BJJ_M., 1995: A 'Lagrangian' revision of the resistors in the two-layer model for calculating the energy budget of a plant canopy. Boundary-Layer Meteorology 74, 261-288.
• Hicks, BB., Baldocchi, DD., Meyers, TP., Hosker, JR., Matt, D_R., 1987: A preliminary multiple resistance routine for deriving dry deposition velocities from measured quantities. Water, Air, and Soil Pollution 36, 311-330.

Gb_Su(Tair,pressure,ustar; wind_zh, Dl, fc, N=2, Cd=0.2, hs=0.01, constants)
Gb_Su!(df; wind_zh, Dl, fc=nothing, N=2, Cd=0.2, hs=0.01, LAI, constants)

Estimate Boundary Layer Conductance to heat transfer using the physically based formulation according to Su et al. 2001.

Arguments

• Tair : Air temperature (degC)
• pressure : Atmospheric pressure (kPa)
• ustar : Friction velocity (m s-1)
• df : DataFrame or matrix containing the above variables
• Dl : Leaf characteristic dimension (m)
• fc : Fractional vegetation cover 0-1
• LAI : One-sided leaf area index (-) - alternative to fc.
• N : Number of leaf sides participating in heat exchange (defaults to 2)
• Cd : Foliage drag coefficient (-)
• hs : Roughness height of the soil (m)
• constants=bigleaf_constants()

Value

see compute_Gb!

Details

The formulation is based on the kB-1 model developed by Massman 1999. Su et al. 2001 derived the following approximation:

$k_{B1} = (k C_d f_c^2) / (4C_t u^*/u(z_h)) + k_{Bs-1}(1 - f_c)^2$

If $f_c$ (fractional vegetation cover) is missing, it is estimated from LAI: $f_c = 1 - e^{-LAI/2}$

The wind speed at the top of the canopy is calculated using function wind_profile.

Ct is the heat transfer coefficient of the leaf (Massman 1999):

$C_t = P_r^{-2/3} R_{eh}^{-1/2} N$

where $P_r$ is the Prandtl number (set to 0.71), and $R_{eh}$ is the Reynolds number for leaves:

$R_{eh} = D_l \, wind(z_h) / v$

k{Bs-1}, the k{B-1} value for bare soil surface, is calculated according to Su et al. 2001:

$k_{Bs-1} = 2.46(Re)^{0.25} - ln(7.4)$

References

• Su, Z., Schmugge, T., Kustas, W. & Massman, W., 2001: An evaluation of two models for estimation of the roughness height for heat transfer between the land surface and the atmosphere. Journal of Applied Meteorology 40, 1933-1951.
• Massman, W., 1999: A model study of kB H- 1 for vegetated surfaces using localized near-field' Lagrangian theory. Journal of Hydrology 223, 27-43.
• Hicks, BB., Baldocchi, DD., Meyers, TP., Hosker, JR., Matt, D_R., 1987: A preliminary multiple resistance routine for deriving dry deposition velocities from measured quantities. Water, Air, and Soil Pollution 36, 311-330.

using DataFrames
df = DataFrame(Tair=25,pressure=100,wind=[3,4,5],ustar=[0.5,0.6,0.65],H=[200,230,250]) 
zh = 25; zr = 40
z0m = roughness_parameters(
  Val(:wind_profile), df.ustar, df.wind, df.Tair, df.pressure, df.H; zh, zr).z0m 
wind_zh = wind_profile(zh, df, 0.7*zh, z0m)
compute_Gb!(df,Val(:Su_2001); wind_zh, Dl=0.01, LAI=5)
# the same meteorological conditions, but larger leaves
compute_Gb!(df,Val(:Su_2001); wind_zh, Dl=0.1,LAI=5)
# same conditions, large leaves, and sparse canopy cover (LAI = 1.5)
compute_Gb!(df,Val(:Su_2001); wind_zh, Dl=0.1,LAI=1.5)
≈(df.Gb_h[1], 0.0638, rtol=1e-3)

Gb_Thom(ustar; constants)
compute_Gb!(df, Val{:Thom_1972})

Boundary Layer Conductance according to Thom 1972, an empirical formulation for the for heat transfer based on a simple ustar (friction velocity) dependency.

Arguments

• ustar : Friction velocity (m s-1)
• df : DataFrame with above variables
• constants=bigleaf_constants()

Details

The empirical equation for Rb suggested by Thom 1972 is:

$Rb = 6.2 {u^*}^{-0.67}$

Gb (=1/Rb) for water vapor and heat are assumed to be equal in this package.

Value

see compute_Gb!

References

• Thom, A., 1972: Momentum, mass and heat exchange of vegetation. Quarterly Journal of the Royal Meteorological Society 98, 124-134.
• Hicks, BB., Baldocchi, DD., Meyers, TP., Hosker, JR., Matt, D_R., 1987: A preliminary multiple resistance routine for deriving dry deposition velocities from measured quantities. Water, Air, and Soil Pollution 36, 311-330.

using DataFrames
df = DataFrame(ustar = SA[0.1,missing,0.3])
compute_Gb!(df, Val(:Thom_1972))
propertynames(df) == [:ustar, :Gb_h]

Gb_constant_kB1(ustar, kB_h; constants)
compute_Gb!(df, Val{:constant_kB1})

Boundary Layer Conductance using constant kB-1 value for heat transfer.

Arguments

• ustar : Friction velocity (m s-1)
• df : DataFrame with above variables
• kB_h : kB-1 value for heat transfer
• constants=bigleaf_constants()

Details

Rbh computed by ``kBh/(k * ustar)``, where k is the von Karman constant.

LE_to_ET(LE,Tair)
ET_to_LE(ET,Tair)

Convert evaporative water flux from mass (ET=evapotranspiration) to energy (LE=latent heat flux) units, or vice versa.

Arguments

• LE Latent heat flux (W m-2)
• ET Evapotranspiration (kg m-2 s-1)
• Tair Air temperature (deg C)

Details

The conversions are given by:

• $ET = LE/\lambda$
• $LE = \lambda ET$

where $\lambda$ is the latent heat of vaporization (J kg-1) as calculated by latent_heat_vaporization.

Examples

# LE of 200 Wm-2 and air temperature of 25degC
ET = LE_to_ET(200,25)
≈(ET, 8.19e-5, atol =1e-7)

Monin_Obukhov_length(Tair, pressure, ustar, H; constants)
Monin_Obukhov_length!(df;constants=bigleaf_constants())
Monin_Obukhov_length(df;constants=bigleaf_constants())

calculates the Monin-Obukhov length.

Arguments

• Tair : Air temperature (degC)
• pressure : Atmospheric pressure (kPa)
• ustar : Friction velocity (m s-1)
• H : Sensible heat flux (W m-2)
• df : DataFrame containing the above variables

optional

• constants=bigleaf_constants(): Dictionary with entries
  • Kelvin - conversion degree Celsius to Kelvin
  • cp - specific heat of air for constant pressure (J K-1 1)
  • k - von Karman constant (-)
  • g - gravitational acceleration (m s-2)

Details

The Monin-Obukhov length (L) is given by:

$L = - (\rho * cp * ustar^3 * Tair) / (k * g * H)$

where $\rho$ is air density (kg m-3).

Note

Note that L gets very small for very low ustar values with implications for subsequent functions using L as input. It is recommended to filter data and exclude low ustar values (ustar < ~0.2) beforehand.

Value

Monin-Obukhov length L (m). The non-mutating DataFrame variant returns a vector, the mutating variant add or modifies column :MOL.

References

Foken, T, 2008: Micrometeorology. Springer, Berlin, Germany.

See also

stability_parameter

Monin_Obukhov_length(
  Tair=25,pressure=100,
  ustar=seq(0.2,1,0.1),H=seq(40,200,20))

Rg_to_PPFD(Rg,J_to_mol=4.6,frac_PAR=0.5)
PPFD_to_Rg(PPFD,J_to_mol=4.6,frac_PAR=0.5)

Conversions between Global Radiation (W m-2) and Photosynthetic Photon Flux Density (umol m-2 s-1)

Arguments

• Rg: Global radiation = incoming short-wave radiation at the surface (W m-2)
• PPFD: Photosynthetic photon flux density (umol m-2 s-1)
• Jtomol: Conversion factor from J m-2 s-1 (= W m-2) to umol (quanta) m-2 s-1
• frac_PAR: Fraction of incoming solar irradiance that is photosynthetical- active radiation (PAR); defaults to 0.5

Details

The conversion is given by:

$PPFD = Rg * frac_PAR * J_to_mol$

by default, the combined conversion factor (frac_PAR * J_to_mol) is 2.3

Examples

# convert a measured incoming short-wave radiation of 500 Wm-2 to
# PPFD in umol m-2 s-1 and backwards
Rg_to_PPFD(500)
PPFD_to_Rg(1150)

Reynolds_Number(Tair,pressure,ustar,z0m; constants)

calculates the Roughness Reynolds Number.

Arguments

• Tair : Air temperature (deg C)
• pressure : Atmospheric pressure (kPa)
• ustar : Friction velocity (m s-1)
• z0m : Roughness length (m)

optional

• constants=bigleaf_constants()

Details

The Roughness Reynolds Number is calculated as in Massman 1999a: $Re = z0m * ustar / v$, where v is the kinematic viscosity (m2 s-1).

Value

Roughness Reynolds Number (-)

References

Massman, W_J., 1999a: A model study of kB H- 1 for vegetated surfaces using 'localized near-field' Lagrangian theory. Journal of Hydrology 223, 27-43.

Tair,pressure,ustar,z0m = 25,100,0.5,0.5
R = Reynolds_Number(Tair,pressure,ustar,z0m)                             
≈(R, 15870, rtol=1e-3)

Rg_to_PPFD(Rg,J_to_mol=4.6,frac_PAR=0.5)
PPFD_to_Rg(PPFD,J_to_mol=4.6,frac_PAR=0.5)

Conversions between Global Radiation (W m-2) and Photosynthetic Photon Flux Density (umol m-2 s-1)

Arguments

• Rg: Global radiation = incoming short-wave radiation at the surface (W m-2)
• PPFD: Photosynthetic photon flux density (umol m-2 s-1)
• Jtomol: Conversion factor from J m-2 s-1 (= W m-2) to umol (quanta) m-2 s-1
• frac_PAR: Fraction of incoming solar irradiance that is photosynthetical- active radiation (PAR); defaults to 0.5

Details

The conversion is given by:

$PPFD = Rg * frac_PAR * J_to_mol$

by default, the combined conversion factor (frac_PAR * J_to_mol) is 2.3

Examples

# convert a measured incoming short-wave radiation of 500 Wm-2 to
# PPFD in umol m-2 s-1 and backwards
Rg_to_PPFD(500)
PPFD_to_Rg(1150)

VPD_to_rH(VPD,Tair; ...)
H_to_VPD(rH,Tair; ...)
e_to_rH(e,Tair; ...)
VPD_to_e(VPD,Tair; ...)
e_to_VPD(e,Tair; ...)
e_to_q(e,pressure; ...)
q_to_e(q,pressure; ...)
q_to_VPD(q,Tair,pressure; ...)
VPD_to_q(VPD,Tair,pressure; ...)

Conversion between vapor pressure (e), vapor pressure deficit (VPD), specific humidity (q), and relative humidity (rH).

Arguments

• Tair: Air temperature (deg C)
• pressure: Atmospheric pressure (kPa)
• e: Vapor pressure (kPa)
• q: Specific humidity (kg kg-1)
• VPD: Vapor pressure deficit (kPa)
• rH: Relative humidity (-)

All functions accept the optional arguemtns:

• Esat_formula: formula used in Esat_from_Tair
• constants: dictionary from bigleaf_constants with entries eps and Pa2kPa

Rreferences

Foken, T, 2008: Micrometeorology_ Springer, Berlin, Germany.

VPD_to_rH(VPD,Tair; ...)
H_to_VPD(rH,Tair; ...)
e_to_rH(e,Tair; ...)
VPD_to_e(VPD,Tair; ...)
e_to_VPD(e,Tair; ...)
e_to_q(e,pressure; ...)
q_to_e(q,pressure; ...)
q_to_VPD(q,Tair,pressure; ...)
VPD_to_q(VPD,Tair,pressure; ...)

Conversion between vapor pressure (e), vapor pressure deficit (VPD), specific humidity (q), and relative humidity (rH).

Arguments

• Tair: Air temperature (deg C)
• pressure: Atmospheric pressure (kPa)
• e: Vapor pressure (kPa)
• q: Specific humidity (kg kg-1)
• VPD: Vapor pressure deficit (kPa)
• rH: Relative humidity (-)

All functions accept the optional arguemtns:

• Esat_formula: formula used in Esat_from_Tair
• constants: dictionary from bigleaf_constants with entries eps and Pa2kPa

Rreferences

Foken, T, 2008: Micrometeorology_ Springer, Berlin, Germany.

VPD_to_rH(VPD,Tair; ...)
H_to_VPD(rH,Tair; ...)
e_to_rH(e,Tair; ...)
VPD_to_e(VPD,Tair; ...)
e_to_VPD(e,Tair; ...)
e_to_q(e,pressure; ...)
q_to_e(q,pressure; ...)
q_to_VPD(q,Tair,pressure; ...)
VPD_to_q(VPD,Tair,pressure; ...)

Conversion between vapor pressure (e), vapor pressure deficit (VPD), specific humidity (q), and relative humidity (rH).

Arguments

• Tair: Air temperature (deg C)
• pressure: Atmospheric pressure (kPa)
• e: Vapor pressure (kPa)
• q: Specific humidity (kg kg-1)
• VPD: Vapor pressure deficit (kPa)
• rH: Relative humidity (-)

All functions accept the optional arguemtns:

• Esat_formula: formula used in Esat_from_Tair
• constants: dictionary from bigleaf_constants with entries eps and Pa2kPa

Rreferences

Foken, T, 2008: Micrometeorology_ Springer, Berlin, Germany.

add_Ga(Gb_h, Ga_m, Sc::Vararg{Pair,N}; constants)
add_Ga!(df::AbstractDataFrame, Sc; Gb_h = df.Gb_h, Ga_m = df.Ga.m, kwargs...)

compute additional aerodynamic conductance quantities for given Schmidt-numbers

Arguments

• Gb_h : Boundary layer conductance for heat transfer (m s-1)
• Ga_m : Aerodynamic conductance for momentum (m s-1)
• Sc : several Pair{Symbol,Number} Output name and Schmidt number of additional conductances to be calculated
• df : DataFrame to add output columns

optional

• constants=bigleaf_constants(): Dictionary with entries
  • Pr - Prandtl number

Details

Aerodynamic conductance is calculated as

$G_{a_x} = 1/(1/G_{a_m} + 1/G_{b_x})$

where Gb_x is the Boundary layer conductance for other quantities x is calculated based on boundary layer for heat transfer, Schmidt-Number, and Prantl number, as documented in add_Gb!.

Value

a NameTuple or df with keys Ga_x where x are the keys in Sc and corresponding aerodynamic conductances (m s-1).

Examples

using DataFrames
df = DataFrame(Gb_h=[0.02, missing, 0.055], Ga_m = 1 ./ [0.03, 0.03, 0.03])
add_Ga!(df, :O2 => 0.84, :CH4 => 0.99)
propertynames(df)[3:4] == [:Ga_O2, :Ga_CH4]

add_Gb(Gb_h::Union{Missing,Number}, Sc::Vararg{Pair,N}; constants)
add_Gb!(df::AbstractDataFrame, Sc::Vararg{Pair,N}; Gb_h = df.Gb_h, kwargs...)

compute boundary layer conductance for additional quantities for given Schmidt-numbers

Arguments

• Gb_h : Boundary layer conductance for heat transfer (m s-1)
• Sc : several Pair{Symbol,Number} Output name and Schmidt number of additional conductances to be calculated
• df : DataFrame to add output columns

optional

• constants=bigleaf_constants(): Dictionary with entries
  • Pr - Prandtl number

Details

Boundary layer conductance for other quantities x is calculated based on boundary layer for heat transfer as (Hicks et al. 1987):

$Gb_x = Gb_h / (Sc_x / Pr)^{0.67}$

where Sc_x is the Schmidt number of quantity x, and Pr is the Prandtl number (0.71).

Value

a NameTuple or df with keys Gb_x where x are the keys in Sc and corresponding boundary layer conductances (m s-1).

Examples

using DataFrames
df = DataFrame(Gb_h=[0.02, missing, 0.055])
add_Gb!(df, :O2 => 0.84, :CH4 => 0.99)
propertynames(df)[2:3] == [:Gb_O2, :Gb_CH4]

aerodynamic_conductance!(df; 
  Gb_model = Val(:Thom_1972), Ram_model = Val(:wind_zr),
  zr=nothing,zh=nothing, d = isnothing(zh) ? nothing : 0.7*zh,
  ...
  )

Bulk aerodynamic conductance, including options for the boundary layer conductance formulation and stability correction functions.

Arguments

• df: DataFrame with columns
  • ustar : Friction velocity (m s-1)
  • wind : Wind speed at sensor height (m s-1)
• Gb_model : model for computing boundary layer conductance (see compute_Gb!)
• Ram_model : model for computing aerodynamic resistance (see compute_Ram)
• zh : canopy height (m)
• zr : Instrument (reference) height (m)

Further required columns of df and keyword argument depend on Gb_model (see compute_Gb!) and Ram_model (see compute_Ram).

If only columns ustar and wind are available, use default models (Val(:Thom_1972) and Val(:wind_zr)).

Details

Aerodynamic conductance for heat (Ga_h) is calculated as:

$Ga_h = 1 / (Ra_m + Rb_h)$

where $Ra_m$ is the aerodynamic resistance for momentum and $Rb_h = 1/Gb_h$ the (quasi-laminar) canopy boundary layer resistance ('excess resistance') for heat.

Ra_m is computed and described with compute_Ram using model Ram_model.

Rb_h is computed and described with 1/compute_Gb! using a given Gb_model.

Value

combined results of compute_Gb! and

• Ra_m: Aerodynamic resistance for momentum transfer (s m-1)
• Ga_m: Aerodynamic conductance for momentum transfer (m s-1)
• Ga_h: Aerodynamic conductance for heat transfer (m s-1)
• Ra_h: Aerodynamic resistance for heat transfer (s m-1)
• Ga_CO2: Aerodynamic conductance for CO2 transfer (m s-1)

Note

The roughness length for water and heat (z0h) can be computed by roughness_z0h.

TODO check Input variables such as LAI, Dl, or zh can be either constants, or vary with time, i.e. are vectors of the same length as df.

Note that boundary layer conductance to water vapor transfer (Gb_w) is often assumed to equal Gb_h. This assumption is also made in Bigleaf.jl, for example in the function surface_conductance.

If the roughness length for momentum (z0m) is not provided as input, it is estimated using roughness_parameters, which estimates a single z0m value for the entire time period. If a varying z0m value (e.g. across seasons or years) is required, z0m should be provided as input argument.

Examples

using DataFrames
df = DataFrame(Tair=25,pressure=100,wind=[3,4,5],
  ustar=[0.5,0.6,0.65],H=[200,230,250])   
# simple calculation of Ga  
aerodynamic_conductance!(df;Gb_model=Val(:Thom_1972)) 
# calculation of Ram using a model derived from the logarithmic wind profile
aerodynamic_conductance!(df;Gb_model=Val(:Thom_1972),Ram_model = Val(:wind_profile), 
  zr=40,zh=25,d=17.5,z0m=2) 
# simple calculation of Ga, but a physically based canopy boundary layer model
aerodynamic_conductance!(df,Gb_model=Val(:Su_2001),
  zr=40,zh=25,d=17.5,Dl=0.05,N=2,fc=0.8)
all(isfinite.(df.psi_h))

air_density(Tair,pressure; ...)

Air density of moist air from air temperature and pressure_

Arguments

• Tair: Air temperature (deg C)
• pressure: Atmospheric pressure (kPa)

optional

• constants=bigleaf_constants(): Dictionary with entries Kelvin, kPa2Pa

Details

Air density $\rho$ is calculated as:

$\rho = {pressure \over Rd * Tair}$

Value

air density (kg m-3)

Examples

# air density at 25degC and standard pressure (101.325kPa)
air_density(25,101.325)

References

Foken, T, 2008: Micrometeorology. Springer, Berlin, Germany.

bigleaf_constants(;...)

Constants Used in the Bigleaf.jl Package

Return a Dictionary Dict{Symbol, Real}. Default values can be overridden by the named arguments:

• cp Specific heat of air for constant pressure (J K-1 kg-1)
• Rgas Universal gas constant (J mol-1 K-1)
• Rv Gas constant of water vapor (J kg-1 K-1) (Stull 1988 p.641)
• Rd Gas constant of dry air (J kg-1 K-1) (Foken p. 245)
• Md Molar mass of dry air (kg mol-1)
• Mw Molar mass of water vapor (kg mol-1)
• eps Ratio of the molecular weight of water vapor to dry air (=Mw/Md)
• g Gravitational acceleration (m s-2)
• solar_constant Solar constant (W m-2)
• pressure0 Reference atmospheric pressure at sea level (Pa)
• Tair0 Reference air temperature (K)
• k von Karman constant
• Cmol Molar mass of carbon (kg mol-1)
• Omol Molar mass of oxygen (kg mol-1)
• H2Omol Molar mass of water (kg mol-1)
• sigma Stefan-Boltzmann constant (W m-2 K-4)
• Pr Prandtl number
• Sc_CO2 Schmidt number for CO2
• Kelvin Conversion degree Celsius to Kelvin
• DwDc Ratio of the molecular diffusivities for water vapor and CO2
• days2seconds Seconds per day
• kPa2Pa Conversion kilopascal (kPa) to pascal (Pa)
• Pa2kPa Conversion pascal (Pa) to kilopascal (kPa)
• umol2mol Conversion micromole (umol) to mole (mol)
• mol2umol Conversion mole (mol) to micromole (umol)
• kg2g Conversion kilogram (kg) to gram (g)
• g2kg Conversion gram (g) to kilogram (kg)
• kJ2J Conversion kilojoule (kJ) to joule (J)
• J2kJ Conversion joule (J) to kilojoule (kJ)
• se_median Conversion standard error (SE) of the mean to SE of the median
• frac2percent Conversion between fraction and percent

calc_sun_position_MOD(JD::Number)

Compute the Sun position at the Julian Day JD.

Results are represented in the Mean Equinox of Date (MOD), i.e. accounting for precession but not for nutation and smaller pertubation of the polar axes, in spherical ecliptic and equatorial coordinates. The algorithm was adapted from [Vallado 2013, p. 277-279].

Arguments:

• JD: time given as Julian Day .

Value

NamedTuple: sun position where

• Ecliptic coordinates (1:3)
  • λ: Ecliptic longitude of the Sun [rad].
  • β: Ecliptic latitude of the Sun [rad] is assumed 0.
  • r: Distance of the Sun from Earth [m].
• Equatorial coordinate (4:5)
  • α: ascention [rad]
  • δ: declination [rad]
• ϵ: Obliquity of the ecliptic [rad].

References

• Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications. Microcosm Press, Hawthorne, CA.

calc_sun_position_hor(datetime, lat, long)

Compute the Sun position at given time and observer coordinates in horizontal coordinates.

Arguments:

• datetime: time: Either a ZonedDateTime, or DateTime assumed in UTC
• lat, long: latitude and longitude in degree

Value

NamedTuple: sun position with entries

• altitude: angle above the horizon [rad].
• azimuth: angle ange the horizon plane eastwards of north [rad]
• hourangle: [rad] as output by AstroLib.eq2hor Seems to represent time [day/2pi] after solar noon. Value at local timezone noon provdes (local time - solar time).

calc_sun_position_hor(datetime, lat, long)

Compute the Sun position at given time and observer coordinates in horizontal coordinates.

Arguments:

• datetime: time: Either a ZonedDateTime, or DateTime assumed in UTC
• lat, long: latitude and longitude in degree

Value

NamedTuple: sun position with entries

• altitude: angle above the horizon [rad].
• azimuth: angle ange the horizon plane eastwards of north [rad]
• hourangle: [rad] as output by AstroLib.eq2hor Seems to represent time [day/2pi] after solar noon. Value at local timezone noon provdes (local time - solar time).

compute_Gb!(df::AbstractDataFrame, approach; kwargs...)

Estimate boundary layer conductance.

Arguments

• df : DataFrame with required variables (depend on approach)
• approach : one of
  • Val(:Thom_1972): see Gb_Thom
  • Val(:Choudhury_1988): see Gb_Choudhury
  • Val(:Su_2001): see Gb_Su
  • Val(:constant_kB1): see Gb_constant_kB1

The different approaches require different variables to be present in df and different keyword arguments.

Value

updated DataFrame df with the following columns:

• Gb_h: Boundary layer conductance for heat transfer (m s-1)

To subsequently derived quantities see

• compute_Gb_quantities for Resistance, kB-1 constant, and CO2 conductance
• add_Gb! for conductances of other species given their Schmidt numbers.

See also

Gb_Thom, Gb_Choudhury, Gb_Su, Gb_constant_kB1, aerodynamic_conductance!

using DataFrames
df = DataFrame(wind=[3,4,5], ustar=[0.5,0.6,0.65]) 
compute_Gb!(df, Val(:Thom_1972))
≈(df.Gb_h[1], 0.102, rtol=1e-2)

compute_Gb_quantities(Gb_h)
compute_Gb_quantities!(df:::AbstractDataFrame)

Based on boundary layer conductance for heat, compute derived quantities.

Arguments

• Gb_h : Boundary layer conductance for heat transfer (m s-1)
• df : DataFrame with above columns
• constants=bigleaf_constants(): entries Sc_CO2 and Pr

Value

NamedTuple with entries

• Gb_h: Boundary layer conductance for heat transfer (m s-1)
• Rb_h: Boundary layer resistance for heat transfer (s m-1)
• kB_h: kB-1 parameter for heat transfer
• Gb_CO2: Boundary layer conductance for CO2 (m s-1).

compute_Ram(::Val{:wind_profile}, ustar; 
  zr, d, z0m, psi_h, constants=bigleaf_constants())
compute_Ram!(df, method::Val{:wind_profile};  
  zr, d, z0m, psi_h = df.psi_h, kwargs...)

compute_Ram(::Val{:wind_zr}, ustar, wind)
compute_Ram!(df, method::Val{:wind_zr}; kwargs...)

Estimate bulk aerodynamic conductance.

Arguments

• ustar : Friction velocity (m s-1)
• df : DataFrame with above columns
• zr : Instrument (reference) height (m)
• d : Zero-plane displacement height (-), can be estimated using roughness_parameters
• z0m : Roughness length for momentum (m). Can be estimated using from roughness_parameters
• psi_h : the value of the stability function for heat and water vapor (-) see stability_correction

Details

The aerodynamic resistance for momentum $R_{a_m}$ is given by (Ram_method = Val(:wind_zr)):

$R_{a_m} = u/{u^*}^2$

Where u is the horizontal wind velocity. Note that this formulation accounts for changes in atmospheric stability, and does not require an additional stability correction function.

An alternative method to calculate $Ra_m$ is provided (Ram_method = Val(:wind_profile)):

$R_{a_m} = (ln((z_r - d)/z_{0m}) - \psi_h) / (k \, u^*)$

If the roughness parameters z0m and d are unknown, they can be estimated using roughness_parameters. The argument stab_formulation determines the stability correction function used to account for the effect of atmospheric stability on Ra_m (Ra_m is lower for unstable and higher for stable stratification). Stratification is based on a stability parameter zeta $\zeta=(z-d/L)$, where z is the height, d the zero-plane displacement height, and L the Monin-Obukhov length, calculated with Monin_Obukhov_length The stability correction function is chosen by the argument stab_formulation. Options are Val(:Dyer_1970) and Val(:Businger_1971) and Val(:no_stability_correction).

Note

For adding aerodynamic conductance for other species see add_Ga!.

Value

Aerodynamic resistance for momentum transfer (s m-1) ($Ra_m$)

References

• Verma, S., 1989: Aerodynamic resistances to transfers of heat, mass and momentum. In: Estimation of areal evapotranspiration, IAHS Pub, 177, 13-20.
• Verhoef, A., De Bruin, H., Van Den Hurk, B., 1997: Some practical notes on the parameter kB-1 for sparse vegetation. Journal of Applied Meteorology, 36, 560-572.
• Hicks, BB., Baldocchi, DD., Meyers, TP., Hosker, JR., Matt, D_R., 1987: A preliminary multiple resistance routine for deriving dry deposition velocities from measured quantities. Water, Air, and Soil Pollution 36, 311-330.
• Monteith, JL., Unsworth, MH., 2008: Principles of environmental physics. Third Edition. Elsevier Academic Press, Burlington, USA.

See also

aerodynamic_conductance!, add_Ga!

TODO; implement decoupling.

This stub is there to satisfy links im Help-pages.

dew_point(Tair,VPD; ...)
dew_point_from_e(ea; ...)

Calculate the dew point, the temperature to which air must be cooled to become saturated (ie e = Esat(Td))

Arguments

• Tair: Air temperature (degC)
• VPD: Vapor pressure deficit (kPa)
• ea: actual water vapor pressure (kPa)

optional

• accuracy = 1e-03: Accuracy of the result (deg C)
• Esat_formula: formula used in Esat_from_Tair
• constants=bigleaf_constants(): Dictionary with entries Pa2kPa

Details

Dew point temperature (Td) is defined by the temperature for which saturaed vapour pressure equals current vapour pressure, i.e. below which water would start to condensate.

$e = Esat(Td)$

where e is vapor pressure of the air and Esat is the vapor pressure deficit. This equation is solved for Td by optimization.

Value

dew point temperature (degC)

References

Monteith J.L., Unsworth M.H., 2008: Principles of Environmental Physics. 3rd edition. Academic Press, London.

Examples

Tair = 20.0
VPD = 1.5
Td = dew_point(Tair, VPD; accuracy = 1e-2)                
(e = VPD_to_e(VPD,Tair), esat_Td = Esat_from_Tair(Td))

dew_point(Tair,VPD; ...)
dew_point_from_e(ea; ...)

Calculate the dew point, the temperature to which air must be cooled to become saturated (ie e = Esat(Td))

Arguments

• Tair: Air temperature (degC)
• VPD: Vapor pressure deficit (kPa)
• ea: actual water vapor pressure (kPa)

optional

• accuracy = 1e-03: Accuracy of the result (deg C)
• Esat_formula: formula used in Esat_from_Tair
• constants=bigleaf_constants(): Dictionary with entries Pa2kPa

Details

Dew point temperature (Td) is defined by the temperature for which saturaed vapour pressure equals current vapour pressure, i.e. below which water would start to condensate.

$e = Esat(Td)$

where e is vapor pressure of the air and Esat is the vapor pressure deficit. This equation is solved for Td by optimization.

Value

dew point temperature (degC)

References

Monteith J.L., Unsworth M.H., 2008: Principles of Environmental Physics. 3rd edition. Academic Press, London.

Examples

Tair = 20.0
VPD = 1.5
Td = dew_point(Tair, VPD; accuracy = 1e-2)                
(e = VPD_to_e(VPD,Tair), esat_Td = Esat_from_Tair(Td))

VPD_to_rH(VPD,Tair; ...)
H_to_VPD(rH,Tair; ...)
e_to_rH(e,Tair; ...)
VPD_to_e(VPD,Tair; ...)
e_to_VPD(e,Tair; ...)
e_to_q(e,pressure; ...)
q_to_e(q,pressure; ...)
q_to_VPD(q,Tair,pressure; ...)
VPD_to_q(VPD,Tair,pressure; ...)

Conversion between vapor pressure (e), vapor pressure deficit (VPD), specific humidity (q), and relative humidity (rH).

Arguments

• Tair: Air temperature (deg C)
• pressure: Atmospheric pressure (kPa)
• e: Vapor pressure (kPa)
• q: Specific humidity (kg kg-1)
• VPD: Vapor pressure deficit (kPa)
• rH: Relative humidity (-)

All functions accept the optional arguemtns:

• Esat_formula: formula used in Esat_from_Tair
• constants: dictionary from bigleaf_constants with entries eps and Pa2kPa

Rreferences

Foken, T, 2008: Micrometeorology_ Springer, Berlin, Germany.

VPD_to_rH(VPD,Tair; ...)
H_to_VPD(rH,Tair; ...)
e_to_rH(e,Tair; ...)
VPD_to_e(VPD,Tair; ...)
e_to_VPD(e,Tair; ...)
e_to_q(e,pressure; ...)
q_to_e(q,pressure; ...)
q_to_VPD(q,Tair,pressure; ...)
VPD_to_q(VPD,Tair,pressure; ...)

Conversion between vapor pressure (e), vapor pressure deficit (VPD), specific humidity (q), and relative humidity (rH).

Arguments

• Tair: Air temperature (deg C)
• pressure: Atmospheric pressure (kPa)
• e: Vapor pressure (kPa)
• q: Specific humidity (kg kg-1)
• VPD: Vapor pressure deficit (kPa)
• rH: Relative humidity (-)

All functions accept the optional arguemtns:

• Esat_formula: formula used in Esat_from_Tair
• constants: dictionary from bigleaf_constants with entries eps and Pa2kPa

Rreferences

Foken, T, 2008: Micrometeorology_ Springer, Berlin, Germany.

VPD_to_rH(VPD,Tair; ...)
H_to_VPD(rH,Tair; ...)
e_to_rH(e,Tair; ...)
VPD_to_e(VPD,Tair; ...)
e_to_VPD(e,Tair; ...)
e_to_q(e,pressure; ...)
q_to_e(q,pressure; ...)
q_to_VPD(q,Tair,pressure; ...)
VPD_to_q(VPD,Tair,pressure; ...)

Conversion between vapor pressure (e), vapor pressure deficit (VPD), specific humidity (q), and relative humidity (rH).

Arguments

• Tair: Air temperature (deg C)
• pressure: Atmospheric pressure (kPa)
• e: Vapor pressure (kPa)
• q: Specific humidity (kg kg-1)
• VPD: Vapor pressure deficit (kPa)
• rH: Relative humidity (-)

All functions accept the optional arguemtns:

• Esat_formula: formula used in Esat_from_Tair
• constants: dictionary from bigleaf_constants with entries eps and Pa2kPa

Rreferences

Foken, T, 2008: Micrometeorology_ Springer, Berlin, Germany.

equilibrium_imposed_ET(Tair,pressure,VPD,Gs, Rn; ...)
equilibrium_imposed_ET!(df; ...)

Evapotranspiration (ET) split up into imposed ET and equilibrium ET.

Argumens

• Tair : Air temperature (deg C)
• pressure : Atmospheric pressure (kPa)
• VPD : Air vapor pressure deficit (kPa)
• Gs : surface conductance to water vapor (m s-1)
• Rn : Net radiation (W m-2)

optional :

• G=0 : Ground heat flux (W m-2)
• S=0 : Sum of all storage fluxes (W m-2)
• Esat_formula=Val(:Sonntag_1990): formula used in Esat_from_Tair
• constants=bigleaf_constants(): pysical constants (cp, eps)

Details

Total evapotranspiration can be written in the form (Jarvis & McNaughton 6):

$ET = \Omega \mathit{ET}_{eq} + (1 - \Omega) \mathit{ET}_{imp}$

where $\Omega$ is the decoupling coefficient as calculated from decoupling. ET_eq is the equilibrium evapotranspiration i.e., the ET rate that would occur under uncoupled conditions, where the budget is dominated by radiation (when Ga -> 0):

$ET_{eq} = (\Delta \, (R_n - G - S) \, \lambda) / ( \Delta \gamma)$

where $\Delta$ is the slope of the saturation vapor pressur(kPa K-1), $\lambda$ is the latent heat of vaporization (J kg-1), and $\gamma$ is the psychrometric constant (kPa K-1). ET_imp is the imposed evapotranspiration rate, the ET rate that would occur under fully coupled conditions (when Ga -> inf):

$ET_{imp} = (\rho \, c_p \, \mathit{VPD} ~ G_s \, \lambda) / \gamma$

where $\rho$ is the air density (kg m-3).

Value

A NamedTuple or DataFrame with the following columns:

• ET_eq: Equilibrium ET (kg m-2 s-1)
• ET_imp: Imposed ET (kg m-2 s-1)
• LE_eq: Equilibrium LE (W m-2)
• LE_imp: Imposed LE (W m-2)

References

• Jarvis, PG., McNaughton, KG., 1986: Stomatal control of transpiration: scaling up from leaf to region. Advances in Ecological Rese1-49.
• Monteith, JL., Unsworth, MH., 2008: Principles of ironmPhysics. 3rd edition. Academic Press, London.

Examples

Tair,pressure,Rn, VPD, Gs = 20.0,100.0,50.0, 0.5, 0.01
ET_eq, ET_imp, LE_eq, LE_imp = equilibrium_imposed_ET(Tair,pressure,VPD,Gs, Rn)    
≈(ET_eq, 1.399424e-05; rtol = 1e-5)

extraterrestrial_radiation(doy::Number; ...)
extraterrestrial_radiation(datetime::TimeType; ...)

Compute the extraterrestrial solar radiation with the eccentricity correction. Computation follows Lanini, 2010 (Master thesis, Bern University).

Arguments

• doy: integer vector with day of year (DoY)

optional

• constants=bigleaf_constants(): Dictionary with entries
  • solar_constant
• year=2030: year to create timestamps. Due to precession results slightly change across decades.

Value

numeric vector of extraterrestrial radiation (W_m-2)

Examples

ex_rad = extraterrestrial_radiation(1)
≈(ex_rad, 1414, atol = 1)

extraterrestrial_radiation(doy::Number; ...)
extraterrestrial_radiation(datetime::TimeType; ...)

Compute the extraterrestrial solar radiation with the eccentricity correction. Computation follows Lanini, 2010 (Master thesis, Bern University).

Arguments

• doy: integer vector with day of year (DoY)

optional

• constants=bigleaf_constants(): Dictionary with entries
  • solar_constant
• year=2030: year to create timestamps. Due to precession results slightly change across decades.

Value

numeric vector of extraterrestrial radiation (W_m-2)

Examples

ex_rad = extraterrestrial_radiation(1)
≈(ex_rad, 1414, atol = 1)

frac_hour(float::Number)
frac_hour(period::Type{<:Period}, float::Number)

Create a period in given type (defaults to Nanosecond) from fractional hours.

using Dates
frac_hour(1+1/60) == Hour(1) + Minute(1)

umolCO2_to_gC(CO2_flux; constants=bigleaf_constants())
gC_to_umolCO2(C_flux; constants=bigleaf_constants())

Convert CO2 quantities from (umol CO2 m-2 s-1) to (g C m-2 d-1) and vice versa.

Arguments

• CO2_flux CO2 flux (umol CO2 m-2 s-1)
• C_flux Carbon (C) flux (gC m-2 d-1)
• constants: dictionary from bigleaf_constants with entries: Cmol, umol2mol, mol2umol, kg2g, g2kg, says2seconds

Examples

umolCO2_to_gC(20)  # gC m-2 d-1

get_datetime_for_doy_hour(doy, hour=12; year = 2030)

Create DateTime for given dayofyear and hour. Hour defaults to noon and year to 2030, a near future where earth axis precession does not differ too much from year where the function is called. Fractional hours can be provided.

Examples

get_datetime_for_doy_hour(2; year = 2021)
# output
2021-01-02T12:00:00

get_growingseason(GPPd, tGPP; ws=15, min_int=5, warngap=true)

Filters annual time series for growing season based on smoothed daily GPP data.

Arguments

• GPPd: daily GPP (any unit)
• tGPP: GPP threshold (fraction of 95th percentile of the GPP time series). Takes values between 0 and 1.

optional

• ws: window size used for GPP time series smoothing
• min_int: minimum time interval in days for a given state of growing season
• warngap: set to false to suppress warning on too few non-missing data

Details

The basic idea behind the growing season filter is that vegetation is considered to be active when its carbon uptake (GPP) is above a specified threshold, which is defined relative to the peak GPP (95th percentile) observed in the year. The GPP-threshold is calculated as:

$GPP_{threshold} = quantile(GPPd,0.95)*tGPP$

GPPd time series are smoothed with a moving average to avoid fluctuations in the delineation of the growing season. The window size defaults to 15 days, but depending on the ecosystem, other values can be appropriate.

The argument min_int serves to avoid short fluctuations in the status growing season vs. no growing season by defining a minimum length of the status. If a time interval shorter than min_int is labeled as growing season or non-growing season, it is changed to the status of the neighboring values, i.e its opposite.

Value

A NamedTuple with entries

• is_growingseason: a BitVector of the same length as the input GPPd in which false indicate no growing season (dormant season) and true indicate growing season.
• GPPd_smoothed: smoothed GPPd

get_nonoverlapping_periods(dfo::DataFrame)

Get non-overlapping periods.

Arguments

• dfo: DataFrame with columns p_start and p_end sorted by increasing p_start

Value

• DataFrame with columns p_start and p_end with p_start[i] > p_end[i-1]

kg_to_mol(mass, molarMass=bigleaf_constants()[:H2Omol])

Conversion between Mass (kg) and Molar Units (mol).

kinematic_viscosity(Tair,pressure; ...)

Calculate the kinematic viscosity of air.

Parameters

• Tair Air temperature (deg C)
• pressure Atmospheric pressure (kPa)

optional

• constants=bigleaf_constants(): Dictionary with entries Kelvin, pressure0, Tair0, kPa2Pa

Details

Eq where v is the kinematic viscosity of the air (m2 s-1), given by (Massman 1999b):

$v = 1.327 * 10^-5(pressure0/pressure)(Tair/Tair0)^{1.81}$

Value

kinematic viscosity of air (m2 s-1)

References

Massman, W.J., 1999b: Molecular diffusivities of Hg vapor in air, O2 and N2 near STP and the kinematic viscosity and thermal diffusivity of air near STP. Atmospheric Environment 33, 453-457.

Examples

Tair,pressure = 25,100
vis = kinematic_viscosity(Tair,pressure)
≈(vis, 1.58e-5, atol =1e-7)

latent_heat_vaporization(Tair)

Latent heat of vaporization as a function of air temperature using

$\lambda = (2.501 - 0.00237 \, Tair) 10^6$.

Arguments:

• Tair: Air temperature (deg C)

Value

$\lambda$: Latent heat of vaporization (J kg-1)

References

• Stull, B, 1988: An Introduction to Boundary Layer Meteorology (p641) Kluwer Academic Publishers, Dordrecht, Netherlands
• Foken, T, 2008: Micrometeorology_ Springer, Berlin, Germany

latent_heat_vaporization.(5:5:45)        

ms_to_mol(G_ms,Tair,pressure; constants=bigleaf_constants())
mol_to_ms(G_mol,Tair,pressure; constants=bigleaf_constants())

Converts conductances from mass (m s-1) to molar units (mol m-2 s-1), or vice versa

Details

The conversions are given by

• $G_{mol} = G_{ms} \, pressure / (Rgas Tair)$
• $G_{ms} = G_{mol} \, (Rgas Tair) / pressure$

where Tair is in Kelvin and pressure in Pa (converted from kPa internally).

References

Jones, HG 1992_ Plants and microclimate: a quantitative approach to environmental plant physiology_ 2nd Edition, Cambridge University Press, Cambridge 428

Examples

G_ms,Tair,pressure = 0.005,25,100
rmol = ms_to_mol(G_ms,Tair,pressure)
≈(rmol, 0.2017, atol =1e-4)

moving_average(vs,n; nmin = n ÷ 2)

Compute the moving average over a vector allowing for missings.

Arguments

• vs: numeric vector to average over
• n: window size: number of items to average
• nmin: minimum number of non-missing records

Values outside the edges are assumed missing.

If the number of non-missing records within a window is smaller than nmin then the averaged value is assumed missing. This avoids average at edges of periods with many missings to be very sensitive to the edge values.

ms_to_mol(G_ms,Tair,pressure; constants=bigleaf_constants())
mol_to_ms(G_mol,Tair,pressure; constants=bigleaf_constants())

Converts conductances from mass (m s-1) to molar units (mol m-2 s-1), or vice versa

Details

The conversions are given by

• $G_{mol} = G_{ms} \, pressure / (Rgas Tair)$
• $G_{ms} = G_{mol} \, (Rgas Tair) / pressure$

where Tair is in Kelvin and pressure in Pa (converted from kPa internally).

References

Jones, HG 1992_ Plants and microclimate: a quantitative approach to environmental plant physiology_ 2nd Edition, Cambridge University Press, Cambridge 428

Examples

G_ms,Tair,pressure = 0.005,25,100
rmol = ms_to_mol(G_ms,Tair,pressure)
≈(rmol, 0.2017, atol =1e-4)

potential_ET(::Val{:PriestleyTaylor}, Tair, pressure, Rn; G=0.0, S=0.0, ...)
potential_ET(::Val{:PriestleyTaylor}, Tair, pressure, Rn, G, S; ...)

potential_ET(::Val{:PenmanMonteith}, Tair, pressure, Rn, VPD, Ga_h;
  G=zero(Tair),S=zero(Tair), ...)
fpotential_ET(::Val{:PenmanMonteith}, Tair, pressure, Rn, VPD, Ga_h, G, S; ...)

potential_ET!(df, approach; ...)

Potential evapotranspiration according to Priestley & Taylor 1972 or the Penman-Monteith equation with a prescribed surface conductance.

Arguments

• Tair: Air temperature (degC)
• pressure: Atmospheric pressure (kPa)
• Rn: Net radiation (W m-2)
• VPD: Vapor pressure deficit (kPa)
• Ga: Aerodynamic conductance to heat/water vapor (m s-1)
• df: DataFrame with the above variables
• approach: Approach used: Either Val(:PriestleyTaylor) or Val(:PenmanMonteith).

optional:

• G=0.0: Ground heat flux (W m-2). Defaults to zero.
• S=0.0: Sum of all storage fluxes (W m-2) . Defaults to zero.
• constants=bigleaf_constants(): physical constants (cp, eps, Rd, Rgas)

for PriestleyTaylor:

• alpha = 1.26: Priestley-Taylor coefficient

for PenmanMonteith:

• Gs_pot = 0.6: Potential/maximum surface conductance (mol m-2 s-1);
• Esat_formula: formula used in Esat_from_Tair

Details

Potential evapotranspiration is calculated according to Priestley & Taylor, 1972 (approach = Val(:PriestleyTaylor):

$\mathit{LE}_{pot} = (\alpha \, \Delta \, (Rn - G - S)) / (\Delta + \gamma)$

$\alpha$ is the Priestley-Taylor coefficient, $\Delta$ is the slope of the saturation vapor pressure curve (kPa K-1), and $\gamma$ is the psychrometric constant (kPa K-1).

If approach = Val(:PenmanMonteith), potential evapotranspiration is calculated according to the Penman-Monteith equation:

$\mathit{LE}_{pot} = (\Delta \, (R_n - G - S) + \rho \, c_p \, \mathit{VPD} \; G_a) / (\Delta + \gamma \, (1 + G_a/G_{s pot})$

where $\Delta$ is the slope of the saturation vapor pressure curve (kPa K-1), $\rho$ is the air density (kg m-3), and $\gamma$ is the psychrometric constant (kPa K-1). The value of $G_{s pot}$ is typically a maximum value of $G_s$ observed at the site, e.g. the $90^{th}$ percentile of $G_s$ within the growing season.

Ground heat flux and storage heat flux G or S are provided as optional arguments. In the input-explicit variants, they default to zero. In the data-frame arguments, they default to missing, which results in assuming them to be zero which is displayed in a log-message. Note that in difference ot the bigleaf R package, you explitly need to care for missing values (see examples).

Value

NamedTuple with the following entries:

• ET_pot: Potential evapotranspiration (kg m-2 s-1)
• LE_pot: Potential latent heat flux (W m-2)

References

• Priestley, CHB., Taylor, R_J., 1972: On the assessment of surface heat flux and evaporation using large-scale parameters. Monthly Weather Review 100, 81-92.
• Allen, RG., Pereira LS., Raes D., Smith M., 1998: Crop evapotranspiration - Guidelines for computing crop water requirements - FAO Irrigation and drainage paper 56.
• Novick, K_A., et al. 2016: The increasing importance of atmospheric demand for ecosystem water and carbon fluxes. Nature Climate Change 6, 1023 - 1027.

See also

surface_conductance

Examples

Calculate potential ET of a surface that receives a net radiation of 500 Wm-2 using Priestley-Taylor:

Tair, pressure, Rn = 30.0,100.0,500.0
ET_pot, LE_pot = potential_ET(Val(:PriestleyTaylor), Tair, pressure, Rn)    
≈(ET_pot, 0.000204; rtol = 1e-2)

Calculate potential ET for a surface with known Gs (0.5 mol m-2 s-1) and Ga (0.1 m s-1) using Penman-Monteith:

Tair, pressure, Rn = 30.0,100.0,500.0
VPD, Ga_h, Gs_pot = 2.0, 0.1, 0.5
ET_pot, LE_pot = potential_ET(
  Val(:PenmanMonteith), Tair,pressure,Rn,VPD, Ga_h; Gs_pot)    
# now cross-check with the inverted equation
Gs_ms, Gs_mol = surface_conductance(
  Val(:PenmanMonteith), Tair,pressure,VPD,LE_pot,Rn,Ga_h)
Gs_mol ≈ Gs_pot

DataFrame variant with explicitly replacing missings using coalesce.:

using DataFrames
df = DataFrame(
  Tair = 20.0:1.0:30.0,pressure = 100.0, Rn = 500.0, G = 105.0, VPD = 2.0, 
  Ga_h = 0.1) 
allowmissing!(df, Cols(:G)); df.G[1] = missing
#
# need to provide G explicitly
df_ET = potential_ET!(copy(df), Val(:PriestleyTaylor); G = df.G)    
ismissing(df_ET.ET_pot[1])
#
# use coalesce to replace missing values by zero
df_ET = potential_ET!(
  copy(df), Val(:PriestleyTaylor); G = coalesce.(df.G, zero(df.G)))    
!ismissing(df_ET.ET_pot[1])

potential_radiation(datetime, lat, long)
potential_radiation(doy, hour, lat, long; timezone,year)

Compute potential radiation for given geolocation and time.

Because potential radiation does not change across years (for a given location and daytime and day of year), time can be specified alternatively by doy and hour.

Arguments

• datetime: UTC timestamp.
• lat: Latitude (decimal degrees)
• long: Longitude (decimal degrees)
• doy: day of year (integer starting from 1)
• hour: hour within the day (0.0 .. 24.0 float)

optional

• timezone: Timezone for doy and hour, defaults to "GMT+x" nearest to given longitude.
• year: specific year for doy and hour

Value

vector of potential radiation (W m-2)

Examples

# assume hours in the GMT+x that is closest to given longitude
potrad = potential_radiation(160, 10.5, 51.0, 11.5)

potential_radiation(datetime, lat, long)
potential_radiation(doy, hour, lat, long; timezone,year)

Compute potential radiation for given geolocation and time.

Because potential radiation does not change across years (for a given location and daytime and day of year), time can be specified alternatively by doy and hour.

Arguments

• datetime: UTC timestamp.
• lat: Latitude (decimal degrees)
• long: Longitude (decimal degrees)
• doy: day of year (integer starting from 1)
• hour: hour within the day (0.0 .. 24.0 float)

optional

• timezone: Timezone for doy and hour, defaults to "GMT+x" nearest to given longitude.
• year: specific year for doy and hour

Value

vector of potential radiation (W m-2)

Examples

# assume hours in the GMT+x that is closest to given longitude
potrad = potential_radiation(160, 10.5, 51.0, 11.5)

pressure_from_elevation(elev,Tair,VPD=missing;...)

An estimate of mean pressure at a given elevation as predicted by the hypsometric equation.

Arguments

• elev: Elevation asl_ (m)
• Tair: Air temperature (deg C)
• VPD: Vapor pressure deficit (kPa); optional

optional

• constants=bigleaf_constants(): Dictionary with entries Kelvin, pressure0, Rd, g, Pa2kPa

Details

Atmospheric pressure is approximated by the hypsometric equation:

$pressure = pressure_0 / (exp(g * elevation / (Rd Temp)))$

The hypsometric equation gives an estimate of the standard pressure at a given altitude. If VPD is provided, humidity correction is applied and the virtual temperature instead of air temperature is used_ VPD is internally converted to specific humidity.

Value

Atmospheric pressure (kPa)

References

Stull B_, 1988: An Introduction to Boundary Layer Meteorology. Kluwer Academic Publishers, Dordrecht, Netherlands.

Examples

# mean pressure at 500m altitude at 25 deg C and VPD of 1 kPa
pressure_from_elevation(500,25,1)

psychrometric_constant(Tair,pressure; ...)

Computes the psychrometric 'constant'.

Argumens

• Tair: Air temperature (deg C)
• pressure: Atmospheric pressure (kPa)

optional

• constants=bigleaf_constants(): Dictionary with entries cp, eps

Details

The psychrometric constant ($\gamma$) is given as:

$\gamma = cp * pressure / (eps * \lambda)$

where $\lambda$ is the latent heat of vaporization (J kg-1), as calculated from latent_heat_vaporization.

Value

the psychrometric constant (kPa K-1)

References

Monteith J.L., Unsworth M.H., 2008: Principles of Environmental Physics. 3rd Edition. Academic Press, London.

Examples

psychrometric_constant.(5.0:5.0:25.0, 100)

VPD_to_rH(VPD,Tair; ...)
H_to_VPD(rH,Tair; ...)
e_to_rH(e,Tair; ...)
VPD_to_e(VPD,Tair; ...)
e_to_VPD(e,Tair; ...)
e_to_q(e,pressure; ...)
q_to_e(q,pressure; ...)
q_to_VPD(q,Tair,pressure; ...)
VPD_to_q(VPD,Tair,pressure; ...)

Conversion between vapor pressure (e), vapor pressure deficit (VPD), specific humidity (q), and relative humidity (rH).

Arguments

• Tair: Air temperature (deg C)
• pressure: Atmospheric pressure (kPa)
• e: Vapor pressure (kPa)
• q: Specific humidity (kg kg-1)
• VPD: Vapor pressure deficit (kPa)
• rH: Relative humidity (-)

All functions accept the optional arguemtns:

• Esat_formula: formula used in Esat_from_Tair
• constants: dictionary from bigleaf_constants with entries eps and Pa2kPa

Rreferences

Foken, T, 2008: Micrometeorology_ Springer, Berlin, Germany.

VPD_to_rH(VPD,Tair; ...)
H_to_VPD(rH,Tair; ...)
e_to_rH(e,Tair; ...)
VPD_to_e(VPD,Tair; ...)
e_to_VPD(e,Tair; ...)
e_to_q(e,pressure; ...)
q_to_e(q,pressure; ...)
q_to_VPD(q,Tair,pressure; ...)
VPD_to_q(VPD,Tair,pressure; ...)

Conversion between vapor pressure (e), vapor pressure deficit (VPD), specific humidity (q), and relative humidity (rH).

Arguments

• Tair: Air temperature (deg C)
• pressure: Atmospheric pressure (kPa)
• e: Vapor pressure (kPa)
• q: Specific humidity (kg kg-1)
• VPD: Vapor pressure deficit (kPa)
• rH: Relative humidity (-)

All functions accept the optional arguemtns:

• Esat_formula: formula used in Esat_from_Tair
• constants: dictionary from bigleaf_constants with entries eps and Pa2kPa

Rreferences

Foken, T, 2008: Micrometeorology_ Springer, Berlin, Germany.

VPD_to_rH(VPD,Tair; ...)
H_to_VPD(rH,Tair; ...)
e_to_rH(e,Tair; ...)
VPD_to_e(VPD,Tair; ...)
e_to_VPD(e,Tair; ...)
e_to_q(e,pressure; ...)
q_to_e(q,pressure; ...)
q_to_VPD(q,Tair,pressure; ...)
VPD_to_q(VPD,Tair,pressure; ...)

Conversion between vapor pressure (e), vapor pressure deficit (VPD), specific humidity (q), and relative humidity (rH).

Arguments

• Tair: Air temperature (deg C)
• pressure: Atmospheric pressure (kPa)
• e: Vapor pressure (kPa)
• q: Specific humidity (kg kg-1)
• VPD: Vapor pressure deficit (kPa)
• rH: Relative humidity (-)

All functions accept the optional arguemtns:

• Esat_formula: formula used in Esat_from_Tair
• constants: dictionary from bigleaf_constants with entries eps and Pa2kPa

Rreferences

Foken, T, 2008: Micrometeorology_ Springer, Berlin, Germany.

roughness_parameters(::Val{:canopy_height}    , zh; frac_d=0.7, frac_z0m=0.1)
roughness_parameters(::Val{:canopy_height_LAI}, zh, LAI; cd=0.2, hs=0.01)
roughness_parameters(::Val{:wind_profile}     , ustar, wind, psi_m; 
  zh, zr, d = 0.7*zh, constants)

roughness_parameters(method::Val{:wind_profile}, utar, wind, Tair, pressure, H; 
  zh, zr, d = 0.7*zh, stab_formulation=Val(:Dyer_1970), constants)
roughness_parameters(method::Val{:wind_profile}, df::DFTable; ...)

A approximations of the two roughness parameters displacement height (d) and roughness length for momentum (z0m).

Arguments

• zh : Vegetation height (m)
• constants=bigleaf_constants():

By canopy height:

• frac_d : Fraction of displacement height on canopy height (-)
• frac_z0m : Fraction of roughness length on canopy height (-)

By canopy height and LAI

• LAI : Leaf area index (-)
• cd : Mean drag coefficient for individual leaves. Defaults to 0.2.
• hs : roughness length of the soil surface (m).

By wind profile

• wind : Wind speed at height zr (m s-1)
• ustar : Friction velocity (m s-1)
• zr : Instrument (reference) height (m)
• d : Zero-plane displacement height (-)
• psi_m : value of the stability function for heat, see stability_correction

Another variant estimates of psi_m by stability_correction, which requires further input arguments. For convenience, these arguments can be provided using a DataFrame, however, the result then is not type stable.

Details

The two main roughness parameters, the displacement height (d) and the roughness length for momentum (z0m) can be estimated from simple empirical relationships with canopy height (zh). If method = Val(:canopy_height), the following formulas are used:

$d = frac_d * zh$

$z0m = frac_{z0m} * zh$

where $frac_d$ defaults to 0.7 and $frac_{z0m}$ to 0.1.

Alternatively, d and z0m can be estimated from both canopy height and LAI (If method = Val(:canopy_height_LAI)). Based on data from Shaw & Pereira 1982, Choudhury & Monteith 1988 proposed the following semi-empirical relations:

$X = cd * LAI$

$d = 1.1 * zh * ln(1 + X^{1/4})$

$z0m = hs + 0.3 * zh * X^{1/2}$ for $0 <= X <= 0.2$

$z0m = hs * zh * (1 - d/zh)$ for $0.2 < X$

If method = Val(:wind_profile), z0m is estimated by solving the wind speed profile for z0m:

$z0m = median((zr - d) * exp(-k*wind / ustar - psi_m)$

By default, d in this equation is fixed to 0.7*zh, but can be set to any other value.

Value

a NamedTuple with the following components:

• d: Zero-plane displacement height (m)
• z0m: Roughness length for momentum (m)
• z0m_se: Only if method = wind_profile: Standard Error of the median for z0m (m)

References

• Choudhury, B. J., Monteith J_L., 1988: A four-layer model for the heat budget of homogeneous land surfaces. Q. J. R. Meteorol. Soc. 114, 373-398.
• Shaw, R. H., Pereira, A., 1982: Aerodynamic roughness of a plant canopy: a numerical experiment. Agricultural Meteorology, 26, 51-65.

See also

wind_profile

using DataFrames
# estimate d and z0m from canopy height for a dense (LAI=5) and open (LAI=2) canopy
zh = 25.0
roughness_parameters(Val(:canopy_height_LAI),zh,5)
roughness_parameters(Val(:canopy_height_LAI),zh,2)   
   
# fix d to 0.7*zh and estimate z0m from the wind profile
df = DataFrame(Tair=[25,25,25],pressure=100,wind=[3,4,5],ustar=[0.5,0.6,0.65],H=200)
roughness_parameters(Val(:wind_profile),df;zh,zr=40,d=0.7*zh)

# assume d = 0.8*zh
rp = roughness_parameters(Val(:wind_profile),df;zh,zr=40,d=0.8*zh)
≈(rp.z0m, 0.55, rtol=0.1)

roughness_z0h(z0m, kB_h)

Arguments

• z0m : Roughness length for momentum (m). Can be calculated by roughness_parameters.
• kB_h : kB-1 parameter, Output of aerodynamic_conductance!

Details

The roughness length for water and heat (z0h) is calculated from the relationship (e.g. Verma 1989):

${k_B}_h = ln(z_{0m}/z_{0h})$

it follows:

$z_{0h} = z_{0m} / e^{k_{B_h}}$

References

Verma, S., 1989: Aerodynamic resistances to transfers of heat, mass and momentum. In: Estimation of areal evapotranspiration, IAHS Pub, 177, 13-20.

set_datetime_ydh!(df, timezone = tz"UTC+0"; 
    year = df.year, doy = df.doy, hour = df.hour)

Update column DateTime column datetime according to year, dayOfYear, and fractional hour.

Value

df with updated or added column :datetime moved to the first position.

setinvalid_afterprecip!(df; min_precip = 0.01, hours_after = 24.0)

Set records after precipitation to false in :valid column.

Arguments

• df: DataFrame with columns :datetime and :precip sorted by :datetime in increasing order.

optional:

• min_precip (in mm per timestep): minimum precip to be considered effective precipitation.
• hours_after: time after the precipitation event to be considered invalid

Value

df with modified column :valid.

setinvalid_nongrowingseason!(df, tGPP; kwargs...)

Set non-growseason to false in :valid column.

Arguments

• df: DataFrame with columns :GPP and :datetime
• tGPP: scalar threshold of daily GPP (see get_growingseason)

optional:

• update_GPPd: set to true additionally update :GPPd_smoothed column to results from get_growingseason
• and others passed to get_growingseason

Value

df with modified columns :valid and if :GPPd_smoothed, where all non-growing season records are set to false.

setinvalid_qualityflag!(df; 
  vars=["LE","H","NEE","Tair","VPD","wind"],
  qc_suffix="_qc",
  good_quality_threshold = 1.0,
  missing_qc_as_bad = true,
  setvalmissing = true, 
)

Set records with quality flags indicating problems to false in :valid column.

Arguments

• df: DataFrame with column :GPP

optional

• vars=["LE","H","NEE","Tair","VPD","wind"]: columns to theck for quality
• qc_suffix="_qc": naming of the corresponding quality-flag column
• good_quality_threshold = 1.0: threshold in quality flag up to which data is considered good quality
• missing_qc_as_bad = true: set to false to not mark records with missing quality flag as invalid
• setvalmissing = true: set to false to prevent replacing values in value column corresponding to problematic quality flag to missing.

Value

df with modified :valid and value columns.

Example

using DataFrames
df = DataFrame(
  NEE = 1:3, NEE_qc = [1,1,2],
  GPP = 10:10:30, GPP_qc = [1,missing,1])
setinvalid_qualityflag!(df; vars = ["NEE", "GPP"])
df.valid == [true, false, false]
ismissing(df.GPP[2]) && ismissing(df.NEE[3])

setinvalid_range!(df, var_ranges...; setvalmissing = true, ...)

Set records with values outside specified ranges to false in :valid column.

If their is no limit towards small or large values, supply -Inf or Inf as the minimum or maximum respectively. If there were false values in the :value column before, they are kept false. In addition, set values outside ranges to missing.

Arguments

• df: DataFrame with column :GPP
• var_ranges: Pair Varname_symbol => (min,max): closed valid interval for respective column

optional

• setvalmissing: set to false to prevent replacing values in value column outside ranges to missing.

Value

df with modified :valid and value columns.

using DataFrames
df = DataFrame(NEE = [missing, 2,3], GPP = 10:10:30)
setinvalid_range!(df, :NEE => (-2.0,4.0), :GPP => (8.0,28.0))
df.valid == [false, true, false]
ismissing(df.GPP[3])

stability_correction(zeta; 
  stab_formulation=Val(:Dyer_1970))
stability_correction(z,d, Tair,pressure,ustar,H; constants,
  stab_formulation=Val(:Dyer_1970))
stability_correction!(df; zeta, z, d; 
  stab_formulation=Val(:Dyer_1970), constants = bigleaf_constants())

Integrated Stability Correction Functions for Heat and Momentum

Arguments

• zeta : Stability parameter zeta (-)
• Tair,pressure,ustar,H : see Monin_Obukhov_length
• z,d : see stability_parameter
• df : DataFrame containting the variables required by Monin_Obukhov_length
• stab_formulation : Formulation for the stability function. Either Val(:Dyer_1970), or Val(:Businger_1971) or Val(:no_stability_correction)

In the second and third form computes zeta by stability_parameter and Monin_Obukhov_length and requires respective arguments.

Details

These dimensionless values are needed to correct deviations from the exponential wind profile under non-neutral conditions. The functions give the integrated form of the universal functions. They depend on the value of the stability parameter $\zeta$, a function of heigh z, which can be calculated from the function stability_parameter. The integration of the universal functions is:

$\psi = -x * zeta$

for stable atmospheric conditions ($\zeta$ >= 0), and

$\psi = 2 * log( (1 + y(zeta)) / 2)$

for unstable atmospheric conditions ($\zeta$ < 0).

The different formulations differ in their value of x and y.

Value

a NamedTuple with the following columns:

• psi_h: the value of the stability function for heat and water vapor (-)
• psi_m: the value of the stability function for momentum (-)

References

• Dyer, A_J., 1974: A review of flux-profile relationships. Boundary-Layer Meteorology 7, 363-372.
• Dyer, A. J., Hicks, B_B., 1970: Flux-Gradient relationships in the constant flux layer. Quart. J. R. Meteorol. Soc. 96, 715-721.
• Businger, J_A., Wyngaard, J. C., Izumi, I., Bradley, E. F., 1971: Flux-Profile relationships in the atmospheric surface layer. J. Atmospheric Sci. 28, 181-189.
• Paulson, C_A., 1970: The mathematical representation of wind speed and temperature profiles in the unstable atmospheric surface layer. Journal of Applied Meteorology 9, 857-861. Foken, T, 2008: Micrometeorology. Springer, Berlin, Germany.

Examples

using DataFrames
zeta = -2:0.5:0.5
df2 = DataFrame(stability_correction.(zeta; stab_formulation=Val(:Businger_1971)))                         
propertynames(df2) == [:psi_h, :psi_m]

stability_parameter(z,d,MOL)
stability_parameter(z,d,Tair, pressure, ustar, H; constants)
stability_parameter!(df::AbstractDataFrame; z,d, MOL=nothing, constants)

calculates stability parameter "zeta", a parameter characterizing stratification in the lower atmosphere.

Arguments

• z : height (m)
• d : Zero-plane displacement height (m)
• MOL : Monin-Obukhov-length L (m)
• df : DataFrame containting the variables required by Monin_Obukhov_length

optional

• constants=bigleaf_constants()

In the second variant and if MOL=nothing in the DataFrame variants, MOL is computed by Monin_Obukhov_length.

Details

The stability parameter $\zeta$ is given by:

$\zeta = (z - d) / L$

where L is the Monin-Obukhov length (m), calculated by Monin_Obukhov_length. The displacement height can be estimated from the function roughness_parameters.

Value

$\zeta$: stability parameter (-). The nonmutainting DataFrame variant returns a vector. The mutating variant modifies or adds column [:zeta].

using DataFrames
df = DataFrame(Tair=25, pressure=100, ustar=0.2:0.1:1.0, H=40:20:200)
z=40;d=15
zeta = stability_parameter.(z,d, df.Tair, df.pressure, df.ustar, df.H)
all(zeta .< 0)

  surface_conductance(::Val{:FluxGradient},   Tair,pressure,VPD,LE; constants)
  surface_conductance(::Val{:PenmanMonteith}, Tair,pressure,VPD,LE,Rn,Ga_h; 
        G = 0.0, S = 0.0, Esat_formula=Val(:Sonntag_1990), constants)

  surface_conductance!(df, method::Val{:FluxGradient}; kwargs...)
  surface_conductance!(df, method::Val{:PenmanMonteith}; G = 0.0, S = 0.0, kwargs...)

Calculate surface conductance to water vapor from the inverted Penman-Monteith equation or from a simple flux-gradient approach.

Arguments

• Tair : Air temperature (deg C)
• pressure : Atmospheric pressure (kPa)
• Rn : Net radiation (W m-2)
• df : DataFrame with above variables
• constants=bigleaf_constants(): Dictionary with physical constants

additional for PenmanMonteith

• LE : Latent heat flux (W m-2)
• VPD : Vapor pressure deficit (kPa)
• Ga_h : Aerodynamic conductance towater vapor (m s-1), assumed equal to that of heat
• G=0.0 : Ground heat flux (W m-2), defaults to zero
• S=0.0 : Sum of all storage fluxes (W m-2), defaults to zero
• Esat_formula=Val(:Sonntag_1990): formula used in Esat_from_Tair

Details

For Val(:PenmanMonteith), surface conductance (Gs) in m s-1 is calculated from the inverted Penman-Monteith equation:

$G_s = ( LE \. G_a \, \gamma ) / ( \Delta \, A + \rho \, c_p \, G_a \, VPD - LE \, (\Delta + \gamma ))$

Where $\gamma$ is the psychrometric constant (kPa K-1), $\Delta$ is the slope of the saturation vapor pressure curve (kPa K^-1), and $\rho$ is air density (kg m-3). Available energy (A) is defined as $A = R_n - G - S$.

Ground heat flux and total storage flux can be provided as scalars or vectors of the lenght of the DataFrame in the DataFrame variant. While the bigleaf R package by default converts any missings in G and S to 0, in Bigleaf.jl the caller must take care, e.g. by using G = coalesce(myGvector, 0.0).

For Val(:FluxGradient), Gs (in mol m-2 s-1) is calculated from VPD and ET only:

$Gs = ET/p \, VPD$

where ET is in mol m-2 s-1 and p is pressure. Note that this formulation assumes fully coupled conditions (i.e. Ga = inf). This formulation is equivalent to the inverted form of Eq.6 in McNaughton & Black 1973:

$Gs = LE \, \gamma / (\rho \, c_p \, VPD)$

which gives Gs in m s-1. Note that Gs > Gc (canopy conductance) under conditions when a significant fraction of ET comes from interception or soil evaporation.

If pressure is not available, it can be approximated by elevation using the function pressure_from_elevation

Value

NamedTuple with entries:

• Gs_ms: Surface conductance in m s-1
• Gs_mol: Surface conductance in mol m-2 s-1

References

• Monteith, J., 1965: Evaporation and environment. In Fogg, G. E. (Ed.), The state and movement of water in living organisms (pp.205-234). 19th Symp. Soc. Exp. Biol., Cambridge University Press, Cambridge
• McNaughton, KG., Black, TA., 1973: A study of evapotranspiration from a Douglas Fir forest using the energy balance approach. Water Resources Research 9, 1579-1590.

Examples

Tair,pressure,VPD,LE,Rn,Ga_h,G = (14.8, 97.7, 1.08, 183.0, 778.0, 0.116, 15.6)
Gs = surface_conductance(Val(:PenmanMonteith), Tair,pressure,VPD,LE,Rn,Ga_h;G)
isapprox(Gs.Gs_mol, 0.28, atol=0.1)

umolCO2_to_gC(CO2_flux; constants=bigleaf_constants())
gC_to_umolCO2(C_flux; constants=bigleaf_constants())

Convert CO2 quantities from (umol CO2 m-2 s-1) to (g C m-2 d-1) and vice versa.

Arguments

• CO2_flux CO2 flux (umol CO2 m-2 s-1)
• C_flux Carbon (C) flux (gC m-2 d-1)
• constants: dictionary from bigleaf_constants with entries: Cmol, umol2mol, mol2umol, kg2g, g2kg, says2seconds

Examples

umolCO2_to_gC(20)  # gC m-2 d-1

virtual_temp(Tair,pressure,VPD; ...)

Virtual temperature, defined as the temperature at which dry air would have the same density as moist air at its actual temperature.

Arguments

• Tair: Air temperature (deg C)
• pressure: Atmospheric pressure (kPa)
• VPD: Vapor pressure deficit (kPa)
• Esat_formula: formula used in Esat_from_Tair

optional

• constants=bigleaf_constants(): Dictionary with entries
  • :Kelvin - conversion degree Celsius to Kelvin
  • :eps - ratio of the molecular weight of water vapor to dry air (-)

Details

The virtual temperature is given by:

Tv = Tair / (1 - (1 - eps) e/pressure)

where Tair is in Kelvin (converted internally)_ Likewise, VPD is converted to actual vapor pressure (e in kPa) with VPD_to_e internally.

Value

virtual temperature (deg C)

References

• Monteith J.L., Unsworth M.H., 2008: Principles of Environmental Physics. 3rd edition. Academic Press, London.

Examples

Tair,pressure,VPD = 25,100,1.5
vt = virtual_temp(Tair,pressure,VPD)  
≈(vt, 26.9, atol =0.1)

wetbulb_temp(Tair, pressure, VPD; ...)

Calculate the wet bulb temperature, ie the temperature that the air would have if it was saturated

Arguments

• Tair: Air temperature (deg C)
• pressure: Atmospheric pressure (kPa)
• VPD: Vapor pressure deficit (kPa)

optional

• accuracy: Accuracy of the result (deg C)
• Esat_formula: formula used in Esat_from_Tair
• constants=bigleaf_constants(): Dictionary with entries cp, eps,Pa2kPa

Details

Wet-bulb temperature (Tw) is calculated from the following expression:

$e = Esat(Tw) - gamma* (Tair - Tw)$

The equation is optimized for Tw. Actual vapor pressure e (kPa) is calculated from VPD using VPD_to_e. The psychrometric constant gamma (kPa K-1) is calculated using psychrometric_constant.

Value

wet-bulb temperature (degC)

References

Monteith J.L., Unsworth M.H., 2008: Principles of Environmental Physics. 3rd edition. Academic Press, London.

Examples

wetbulb_temp.([20,25.0], 100, [1,1.6])             

wetbulb_temp(Tair, pressure, VPD; ...)

Calculate the wet bulb temperature, ie the temperature that the air would have if it was saturated

Arguments

• Tair: Air temperature (deg C)
• pressure: Atmospheric pressure (kPa)
• VPD: Vapor pressure deficit (kPa)

optional

• accuracy: Accuracy of the result (deg C)
• Esat_formula: formula used in Esat_from_Tair
• constants=bigleaf_constants(): Dictionary with entries cp, eps,Pa2kPa

Details

Wet-bulb temperature (Tw) is calculated from the following expression:

$e = Esat(Tw) - gamma* (Tair - Tw)$

The equation is optimized for Tw. Actual vapor pressure e (kPa) is calculated from VPD using VPD_to_e. The psychrometric constant gamma (kPa K-1) is calculated using psychrometric_constant.

Value

wet-bulb temperature (degC)

References

Monteith J.L., Unsworth M.H., 2008: Principles of Environmental Physics. 3rd edition. Academic Press, London.

Examples

wetbulb_temp.([20,25.0], 100, [1,1.6])             

wind_profile(z::Number, ustar, d, z0m, psi_m = zero(z); constants)
wind_profile(z, df::AbstractDataFrame, d, z0m, psi_m::AbstractVector; constants)
wind_profile(z, df::DFTable, d, z0m; psi_m = nothing,
  stab_formulation = Val(:Dyer_1970), constants)

Wind speed at a given height above the canopy estimated from single-level measurements of wind speed.

Arguments

• z : Height above ground for which wind speed is calculated.
• ustar : Friction velocity (m s-1)
• d : Zero-plane displacement height (-)
• z0m : Roughness length (m)
• constants=bigleaf_constants()

For DataFrame variant with supplying stability_parameter

• df: : DataFrame with columns
  • ustar : Friction velocity (m s-1)
• psi_m : value of the stability function for heat, see stability_correction Pass psi_m = 0.0 to neglect stability correction.

For DataFrame varinat where psi_m is to be estimated

• additional columns in df:
  • Tair, pressure, H : see stability_correction

Details

The underlying assumption is the existence of a logarithmic wind profile above the height d + z0m (the height at which wind speed mathematically reaches zero according to the Monin-Obhukov similarity theory). In this case, the wind speed at a given height z is given by:

$u(z) = (ustar/k) * (ln((z - d) / z0m) - \psi_m$

The roughness parameters zero-plane displacement height (d) and roughness length (z0m) can be approximated from roughness_parameters. $\psi_m$ is the stability correction. Set it to zero (not recommended) to neglect statbility correction. By default it is estimated from wind profile using stability_correction

Note

Note that this equation is only valid for z >= d + z0m, and it is not meaningful to calculate values closely above d + z0m. All values in heights smaller than d + z0m will return 0.

Value

wind speed at given height z.

References

• Monteith, JL., Unsworth, MH., 2008: Principles of Environmental Physics. 3rd edition. Academic Press, London.
• Newman, JF., Klein, PM., 2014: The impacts of atmospheric stability on the accuracy of wind speed extrapolation methods. Resources 3, 81-105.

See also

roughness_parameters

using DataFrames
heights = 18:2:40  # heights above ground for which to calculate wind speed
df = DataFrame(Tair=25,pressure=100,wind=[3,4,5],ustar=[0.5,0.6,0.65],H=[200,230,250]) 
zr=40;zh=25;d=16
# z0m and MOL are independent of height, compute before
MOL = Monin_Obukhov_length.(df.Tair, df.pressure, df.ustar, df.H)
z0m = roughness_parameters(
  Val(:wind_profile), df.ustar, df.wind, df.Tair, df.pressure, df.H; zh, zr).z0m 
ws = map(heights) do z
  wind_profile(z,df,d,z0m; MOL)
end
using Plots # plot wind profiles for the three rows in df
plot(first.(ws), heights, ylab = "height (m)", xlab = "wind speed (m/s)", legend=:topleft)
plot!(getindex.(ws, 2), heights)
plot!(getindex.(ws, 3), heights)
nothing
# output