Garden with Insight v1.0 Help: Hydrology - Evapotranspiration

Potential Evapotranspiration

The model offers four options for estimating potential evaporation: Hargreaves and Samani (1985), Penman (1948), Priestly-Taylor (1972), and Penman-Monteith (Monteith, 1965). The Penman and Penman-Monteith methods require solar radiation, air temperatue, wind speed, and relative humidity as input. If wind speed, relative humidity, and solar radiation data are not available, the Hargreaves or Priestley-Taylor methods provide options that give realistic results in most cases.

The model computes evaporation from soils and plants separately, as described by Ritchie (1972). Potential soil water evaporation is estimated as a class function of potential evaporation and leaf area index (LAI, area of plant leaves relative to the soil surface area). Actual soil water evaporation is simulated as a linear function of potential evaporation and leaf area index.

Penman method

The Penman (1948) option for estimating potential evaporation is based on the equation [Equation 50] where E(0) is the potential evaporation in mm, delta is the slope of the saturation vapor pressure curvin kPa/degreesC, gamma is a psychrometer constant in kPa/degressC, h(0) is the net radiation in MJ/m2, G is the soil heat flux in MJ/m2, HV is the latent heat of vaporization in MJ/kg, f(V) is a wind speed function in mm/day*kPa, e(a) is the saturation vapor pressure at mean air temperature in kPa, and e(d) is the vapor pressure at mean air temperature in kPa.

Equation 50:

E(o) = (delta / (delta + gamma)) * (h(o) - G) / HV
+ (gamma / (delta + gamma)) * f(V) * (e(a) - e(d))
Code:
same
Variables:
E(o) = PotentialSoilEvapByPenman_mm
delta = slopeSaturVaporPressureCurve_kPaPdegC
gamma = psychromConstant_kPaPdegC
h(0) = netRadiationForDay_MJPm2
G = soilHeatFlux_MJPM2
HV = latentHeatOfVaporization_MJPkg
f(V) = windSpeedFunction_mmPdaykPa
e(a) - e(d) = vaporPressureDeficit_kPa

The latent heat of vaporization is estimated with the temperature function [Equation 51] where T is the mean daily temperature in degrees C.

Equation 51:

HV = 2.5 - 0.0022 * T
Code:
same
Variables:
HV = LatentHeatOfVaporization_MJPkg
T = meanTempForDay_degC

The saturation vapor pressure is also estimated as a function of temperature by using the equation [Equation 52].

Equation 52:

e(a) = 0.1 * exp(54.88 - 5.03 * ln(T + 273) - 6791 / (T + 273))
Code:
same
Variables:
e(a) = saturVaporPressure_kPa
T = meanTempForDay_degC

The vapor pressure is simulated as a function of the saturation value and the relative humidity [Equation 53] where RH is the relative humidity expressed as a fraction.

Equation 53:

e(d) = e(a) * RH
Code:
same
Variables:
e(d) = vaporPressureAtMeanTemp_kPa
e(a) = saturVaporPressureAtMeanTemp_kPa
RH = relHumForDay_frn

The slope of the saturation vapor pressure curve is estimated with the equation [Equation 54].

Equation 54:

delta = (e(a) / (T + 273)) * (6791 / (T + 273) -5.03)
Code:
same
Variables:
delta = SlopeSaturVaporPressureCurve_kPaPdegC
e(a) = saturVaporPressureAtMeanTemp_kPa
T = meanTempForDay_degC

The psychrometer constant is computed with the equation [Equation 55] where PB is the barometric pressure (kPa).

Equation 55:

gamma = 6.6e10-4 * PB
Code:
same
Variables:
gamma = psychromConstant_kPaPdegC
PB = barometricPressure_kPa

The barometric pressure is estimated as a function of elevation by using the equation [Equation 56] where ELEV is the elevation of the site in m.

Equation 56:

PB = 101 - 0.0115 * ELEV + 5.44e10-7 * sqr(ELEV)
Code:
same
Variables:
PB = barometricPressure_kPa
ELEV = siteElevation_m

The soil heat flux is estimated by using air temperature on the day of interest plus three days prior: [Equation 57] where T is the mean daily air temperature on day i in degrees C. Since G is usually relatively small, it is assumed zero in EPIC.

Equation 57:

G = 0.12 * (T(i) - (T(i-1) + T(i-2) + T(i-3)) / 3)
Code:
G is assumed to be zero in the code
We are implementing it because in a microclimate the heat flux might be more important.
Variables:
G = SoilHeatFlux_MJPM2
T(i) = meanTempForDay_degC
T(i-1) = dailyMeanTempMinus1Day_degC
T(i-2) = dailyMeanTempMinus2Days_degC
T(i-3) = dailyMeanTempMinus3Days_degC

Solar radiation is adjusted to obtain net radiation by using the equation [Equation 58] where RA is the solar radiation in MJ/m2, AB is albedo, RAB is the net outgoing long wave radiation in MJ/m2 for clear days, and RAMX is the maximum solar radiation possible in MJ/m2 for the location on day i.

Equation 58:

h(o) = RA * (1.0 - AB) - RAB * ((0.9 * RA) / RAMX + 0.1)
Code:
same
Variables:
h(o) = NetRadiationForDay_MJPm2
RA = radiationForDay_MJPm2
AB = albedo_frn
RAB = netOutgoingRadiationIfDayIsClear_MJPm2
RAMX = maxPossibleRadiation_MJPm2

The RAB value is estimated with the equation [Equation 59].

Equation 59:

RAB = 4.9x10-9 * (0.34 - 0.14 * sqrt(e(d))) * power(T + 273, 4)
Code:
same
Variables:
RAB = netOutgoingRadiationIfDayIsClear_MJPm2
e(d) = vaporPressureAtMeanTemp_kPa
T = meanTempForDay_degC

The maximum possible solar radiation is computed with the equations [Equation 60] and [Equation 61] where LAT is the latitude of the site in degrees, SD is the sun's declination angle (radians), and i is the day of the year.

Equation 60:

RAMX = 30 * (1.0 + 0.0335 * sin((2pi/365) * (i + 88.2)) * (XT * sin((2pi/360) * LAT)
* sin(SD) + cos((2pi/360) * LAT) * cos(SD) * sin(XT)))
Code:
same except in the code this is calculated only once per month
we will calculate it daily
Variables:
RAMX = MaxPossibleRadiation_MJPm2
i = julianDay
XT = maxPossibleRadiationVariable
LAT = siteLatitude_rad
SD = declinationAngleOfSun_rad

Equation 61:

XT = acos(-tan((2pi/360) * LAT) * tan(SD)), 0 <= XT <= pi
Code:
same
Variables:
XT = MaxPossibleRadiationVariable
LAT = siteLatitude_rad
SD = declinationAngleOfSun_rad

Finally, the wind function of the Penman equation is approximated with the relationship [Equation 63] where V is the mean daily wind speed (at a 10 meter height) in meters/second.

Equation 63:

f(V) = 2.7 + 1.63 * V
Code:
same
Variables:
f(V) = PenmanWindSpeedFunction_mmPdaykPa
V = meanWindSpeedForDay_mPsec

Penman-Monteith method

The Penman-Monteith method (Monteith, 1965) was recently added to EPIC to provide a means for estimating the effects of CO2 changes (Stockle et al., 1992). The Penman-Monteith equation is expressed as [Equation 64] where AD is the air density in g/m3 and AR is the aerodynamic resistance for heat and vapor transfer in s/m.

Equation 64:

E(o) = (delta * (h(o) - G) + 86.7 * AD * (e(a) - e(d)) / AR)
/ (HV * (delta + gamma))
Code:
same
Variables:
E(o) = PotentialSoilEvapByPenmanMonteith_mm
delta = slopeSaturVaporPressureCurve_kPaPdegC
h(0) = netRadiationForDay_MJPm2
G = soilHeatFlux_MJPM2
AD = airDensity_gPm3
e(a) - e(d) = vaporPressureDeficit_kPa
AR = aeroResistForHeatAndVaporTransfer_secPm
HV = latentHeatOfVaporization_MJPkg
gamma = psychromConstant_kPaPdegC

The Penman-Monteith method also estimates plant water evaporation (all other methods use a common equation which will come later), which is [Equation 65] where CR is the canopy resistance for vapor transfer in s/m.

Equation 65:

E(p) = (delta * (h(o) - G) + 86.7 * AD * (e(a) - e(d)) / AR)
/ (HV * (delta + gamma * (1.0 + CR / AR)))
Code:
same
Variables:
E(p) = PotentialPlantEvapByPenmanMonteith_mm
delta = slopeSaturVaporPressureCurve_kPaPdegC
h(0) = netRadiationForDay_MJPm2
G = soilHeatFlux_MJPM2
AD = airDensity_gPm3
e(a) - e(d) = vaporPressureDeficit_kPa
AR = aeroResistForHeatAndVaporTransfer_secPm
HV = latentHeatOfVaporization_MJPkg
gamma = psychromConstant_kPaPdegC
CR = canopyResisForVaporTransfer_secPm

Air density is estimated with the equation [Equation 66].

Equation 66:

AD = 0.01276 * PB / (1.0 + 0.0367 * T)
Code:
same
Variables:
AD = AirDensity_gPm3
PB = barometricPressure_kPa
T = meanTempForDay_degC

The aerodynamic resistance is computed with the equation [Equation 67] where ZD is the displacement height of the crop in m, ZO is the surface roughness parameter in m, and V is the daily mean wind speed in m/sec.

Equation 67:

AR = 6.25 * sqr(ln((10.0 - ZD) / ZO)) / V
Code:
same, but if cropHeight_m > 8.0,
the 10.0 is replaced with cropHeight_m + 2
and V is replaced with V * ln((cropHeight_m + 2) / 0.0005) / 9.9035
Variables:
AR = AeroResistForHeatAndVaporTransfer_secPm
ZD = displacementHeightOfCrop_m
ZO = surfaceRoughnessParam_m
V = meanWindSpeedForDay_mPsec

The surface roughness is estimated with the equation [Equation 68] and the crop displacement height is estimated with the equation [Equation 69] where CHT is crop height in m.

Equation 68:

ZO = 0.131 * power(CHT, 0.997)
Code:
ZO = power(10.0, 0.979 * log10(CHT) - 0.883)
Variables:
ZO = SurfaceRoughnessParam_m
CHT = cropHeight_m

Equation 69:

ZD = 0.702 * power(CHT, 0.979)
Code:
ZD = power(10.0, 0.979 * log10(CHT) - 0.154)
Variables:
ZD = DisplacementHeightOfCrop_m
CHT = cropHeight_m

When no crop is growing the aerodynamic resistance is estimated with the equation [Equation 70].

Equation 70:

AR = 350.0 / V
Code:
same
Variables:
AR = AeroResistForHeatAndVaporTransferNoCrop_secPm
V = meanWindSpeedForDay_mPsec

The canopy resistance is computed with the equation [Equation 71] where p(1) is a parameter ranging from 1.0 to 2.0, LAI is the leaf area index of the crop, g*(o) is the leaf conductance in m/sec, and CO2 is the carbon dioxide level in the atmosphere in ppm.

Equation 71:

CR = p(1) / (LAI * g*(o) * (1.4 - 0.00121 * CO2))
Code:
same (0.4 / 330.0 = 0.00121)
Variables:
CR = CanopyResistForVaporTransfer_secPm
p(1) = canopyResistParam
LAI = leafAreaIndex
g*(o) = leafConductance_mPsec
CO2 = carbonDioxideInAtmosphere_ppm

Leaf conductance is estimated from the crop input rate adjusted for vapor pressure deficit (VPD) [Equation 72] where g(o) is the crop's leaf resistance when the VPD is less than the crop's threshold VPD, and FV is the VPD correction factor.

Equation 72:

g*(o) = g(o) * FV
Code:
same
Variables:
g*(o) = LeafConductance_mPsec
g(o) = leafResistIfVaporPressureDeficitBelowThresholdForCrop_mPsec
FV = vaporPressureDeficitCorrectionFactor

FV is calculated by [Equation 73] where b(v) is a crop coefficient and VPD(t) is threshold VPD for the crop.

Equation 73:

FV = 1.0 - b(v) * (VPD - VPD(t)), >= 0.1
Code:
same, except if VPD < VPD(t), FV = 1.0
Variables:
FV = VaporPressureDeficitCorrectionFactor
b(v) = fractionOfMaxLeafConductForHighVPD_frn
VPD = vaporPressureDeficit_kPa
VPD(t) = thresholdVaporPressureDeficit_kPa

Priestley-Taylor method

The Priestley-Taylor (1972) method provides estimates of potential evaporation without wind and relative humidity inputs. The simplified equation based only on temperature and radiation is [Equation 74].

Equation 74:

E(o) = 1.28 * (h(o) / HV) * (delta / (delta + gamma))
Code:
same
Variables:
E(o) = PotentialSoilEvapByPriestleyTaylor_mm
h(o) = netRadiationByPriestleyTaylor_MJPm2
HV = latentHeatOfVaporization_MJPkg
delta = slopeSaturVaporPressureCurve_kPaPdegC
gamma = psychromConstant_kPaPdegC

The net radiation (for the Priestley-Taylor method) is estimated with the equation [Equation 75] instead of equation 58 which is used in the Penman method.

Equation 75:

h(o) = RA * (1.0 - AB)
Code:
same
Variables:
h(o) = NetRadiationByPriestleyTaylor_MJPm2
RA = radiationForDay_MJPm2
AB = albedo_frn

Hargreaves method

The Hargreaves method (Hargreaves and Samani, 1985) estimates potential evapotranspiration as a function of extraterrestrial radiation and air temperature. Hargreaves' method was modified for use in EPIC by increasing the temperature difference exponent from 0.5 to 0.6. Also, extraterrestrial radiation is replaced by RAMX (maximum possible solar radiation at the earth's surface) and the coefficient is adjusted from 0.0023 to 0.0032 for proper conversion. The modified equation is [Equation 76] where T(mx) and T(mn) are the daily maximum and minimum air temperatures in degrees C.

Equation 76:

E(o) = 0.0032 * (RAMX / HV) * (T + 17.8) * power(T(mx) - T(mn), 0.6)
Code:
same
Variables:
E(o) = PotentialSoilEvapByHargreaves_mm
RAMX = maxPossibleRadiation_MJPm2
HV = latentHeatOfVaporization_MJPkg
T = meanTempForDay_degC
T(mx) = maxTempForDay_degC
T(mn) = minTempForDay_degC

All four methods estimate albedo by considering the soil, crop, and snow cover. If a snow cover exists with 5 mm or greater water content, the value of albedo is set to 0.6. If the snow cover is less than 5 mm and no crop is growing, the soil albedo is the appropriate value. When crops are growing, albedo is determined by using the equation [Equation 77] where 0.23 is the albedo for plants, AB(s) is the soil albedo, and EA is a soil cover index.

Equation 77:

if snowWaterContent_mm > 5 AB = 0.6
if snowWaterContent_mm < 5 and no crop is growing, AB = AB(s)
if crops are growing AB = 0.23 * (1.0 - EA) + AB(s) * EA
Code:
same
Variables:
AB = Albedo_frn
EA = soilCoverIndex_frn
AB(s) = soilAlbedo_frn

The value of EA ranges from 0.0 to 1.0 according to the equation [Equation 78] where CV is the sum of the above ground biomass and crop residue in t/ha.

Equation 78:

EA = exp(-0.05 * CV)
Code:
same except if snowWaterContent_mm > 5 EA = 1.0
Variables:
EA = SoilCoverIndex_frn
CV = aboveGroundBiomassAndResidue_tPha

Actual Evapotranspiration

The model computes evaporation from soils and plants separately by an approach similar to that of Ritchie (1972). For all methods except Penman-Monteith, potential plant water evaporation is computed with the equations [Equation 79] and [Equation 80] where E(p) is the predicted plant water evaporation rate in mm/day. If soil water is limited, plant water evaporation will be reduced as described in the plant growth section of this section.

Equations 79, 80:

E(p) = (E(o) * LAI) / 3.0 if 0.0 < LAI < 3.0
E(p) = E(o) if LAI > 3.0
Code:
E(p) = min((E(o) * LAI) / 3.0, E(o))
this comes out the same because if LAI > 3, E(o) * LAI / 3 > E(o)
Variables:
E(p) = potPlantEvap_mm
E(o) = potentialSoilEvap_mm
LAI = leafAreaIndex

Potential soil water evaporation is simulated by considering soil cover according to the following equation [Equation 81] where E(s) is the potential soil water evaporation rate in mm/day.

Equation 81:

E(s) = E(o) * EA
Code:
same
Variables:
E(s) = PotentialSoilEvapAdjByCover_mm
E(o) = potentialSoilEvap_mm
EA = soilCoverIndex_frn

Potential soil water evaporation is reduced during periods of high plant water use with the equation [Equation 82]. When E(p) is low E*(s) approaches E(s), but as E(p) approaches E(o), E*(s) approaches E(s)/(1.0 + EA).

Equation 82:

E*(s) = min(ES, E(s) * E(o) / (E(s) + E(p)))
Code:
same
Variables:
E*(s) = PotentialSoilEvapAdjByCoverAndPotPlantEvap_mm
E(s) = potentialSoilEvapAdjByCover_mm
E(o) = potentialSoilEvap_mm
E(p) = potPlantEvap_mm

Actual soil water evaporation is estimated on the basis of the top 0.2 m of soil and snow cover, if any. If 5 mm or more (water content) of snow is present albedo is set to 0.6 and EA to 0.5 for estimating E(o) and snow is evaporated at that rate. When all snow is evaporated, soil water evaporation begins. Such evaporation is governed by soil depth and water content according to the equation [Equation 83] where EV is the total soil water evaporation in mm from soil of depth Z in mm. The coefficients of equation 83 are set to give EV = 0.5 E*(s) when Z = 10 mm and EV = 0.95 E*(s) when Z = 100 mm.

Equation 83:

EV(Z) = E*(s) * (Z / (Z + exp(2.374 - 0.00713 * Z)))
Code:
same
Variables:
EV(Z) = PotentialSoilEvapForDepth_mm
E*(s) = potentialSoilEvapAdjByCoverAndPotPlantEvap_mm
Z = depth_mm

Potential soil water evaporation for a layer is estimated by taking the difference between EV's at the layer boundaries [Equation 84] where SEV is the potential soil evaporation for layer l in mm.

Equation 84:

SEV = EV(Z(l)) - EV(Z(l-1))
Code:
same
Variables:
SEV = PotentialSoilEvapForLayer_mm

The depth distributed estimate of soil water evaporation may be reduced according to the following equation if soil water is limited in a layer [Equation 85] and [Equation 86] where SEV*(l) is the adjusted soil water evaporation
estimate in mm.

Equations 85, 86:

SEV* = SEV * exp(2.5 * (SW - FC) / (FC - WP)) if SW < FC
SEV* = SEV if SW >= FC
Code:
same
Variables:
SEV* = PotentialSoilEvapForLayerAdjForDryness_mm
SEV = potentialSoilEvapForLayer_mm
SW = waterContent_mm
FC = fieldCapacity_mm
WP = wiltingPoint_mm

The final step in adjusting the evaporation estimate is to assure that the soil water supply is adequate to meet the demand [Equation 87] where b(w) may range from 0.0 to 1.0 in the top 0.5 m of soil and is set to 1.0 below 0.5 m. Thus, EPIC can be adjusted to allow the top 0.5 m to dry down to any fraction of wilting point.

Equation 87:

SEV* = min(SEV*, SW - b(w) * WP)
Code:
same
code has
if (SW - SEV* < b(w) * WP) SEV* = SW - b(W) * WP otherwise SEV* = SEV*
if SW - SEV* - b(w) * WP < 0
if SW - b(w) * WP - SEV* < 0
if SEV* > SW - b(w) * WP
which is the same as min(SEV*, SW - b(w) * WP)
Variables:
SEV* = SoilEvaporationForLayer_mm
SEV* = potentialSoilEvapForLayer_mm
SW = waterContent_mm
b(w) = lowerLimitWaterContentInTopP5MAsFractOfWP_frn
WP = wiltingPoint_mm