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
|