Calculating Daily Penman-Monteith Reference ET

This is a step-by-step example of how reference evapotranspiration (ET_o)is calculated in the ETo package following the FAO ¹ guidelines and the steps outlined in Zotarelli et al (2010)².

ET_o Definition

Reference evapotranspiration (ET_o) is defined as the potential water lost to evaporation and plant transpiration from a reference surface:

The reference surface is a hypothetical grass reference crop with an assumed crop height of 0.12 m, a fixed surface resistance of 70 s m^-1 and an albedo of 0.23. The reference surface closely resembles an extensive surface of green, well-watered grass of uniform height, actively growing and completely shading the ground. The fixed surface resistance of 70 s m^-1 implies a moderately dry soil surface resulting from about a weekly irrigation frequency.³

To calculate ET_o, we use the Food and Agriculture Organization’s Penman-Monteith ET_o algorithm. This method allows ET_o to be calculated with using daily meteorological observations. With this method, daily ET_o is calculated as:

\[ ET_o = \frac{0.408 \Delta (R_n - G) + \gamma \frac{900}{T + 273} u_2 (e_s - e_a)}{\Delta + \gamma (1 + 0.34 u_2)} \]

where:

\(ET_o\) is the reference surface evapotranspiration (mm/day)
\(R_n\) is the net radiation at the crop surface (MJ/m²/day)
\(G\) is soil heat flux (MJ/m²/day)
\(T\) is the average temperature (°C)
\(u_2\) is the wind speed at 2 meters height (m/s)
\(e_s\) is the saturation vapor pressure (kPa)
\(e_a\) is the actual vapor pressure (kPa)
\(\Delta\) is the slope of the vapor pressure curve (kPa/°C)
\(\gamma\) is the psychrometric constant (kPa/°C)

Calculation of Daily ET_o

The steps to derive all the intermediate variables necessary to calculate ET_o are outlined in the following sections. The following input data is used in all of the examples:

Example meteorological data as measured on 6 July in Uccle (Brussels, Belgium)⁴.
Maximum Air Temperature (T_max)	21.5 [°C]
Minimum air temperature (T_min)	12.3 [°C]
Maximum relative humidity (RH_max)	84 [%]
Minimum relative humidity (RH_min)	63 [%]
Average 10 meter wind speed (u₂)	2.7778 [m/s]
Average incoming shortwave radiation (R_s)	255.4398 [W m^-2]
Elevation (z)	100 [m]
Latitude	50.8 [degrees]

Step 1: Calculate Daily Average Temperature

The daily average temperature is derived from T_min and T_max:

\[ T_{mean} = \frac{T_{max} + T_{min}}{2} \] Plugging in the numbers from our example, we get:

\[ T_{mean} = \frac{21.5 + 12.3}{2} = 16.9°C \]

Step 2: Convert Solar Radiation to Units of MJ m^-2 day^-1

One watt per square meter is equal to one joule per meter squared per second. Therefore, we convert daily incoming shortwave radiation by multiplying it by the number of seconds in a day (3600 * 24), then converting from joules to megajoules (multiply by 1e-6): \[ R_s [MJ m^{-2} day^{-1}] = R_s [W m^{2}] * 3600 * 24 * 1e-6 \] Using our example, we get: \[ R_s [MJ m^{-2} day^{-1}] = 255.4398 * 3600 * 24 * 1e-6 = 22.07 [MJ m^{-2} day^{-1}] \] ETo provides the wm2_to_mj() function, which does this calculation for the user.

Step 3: Adjust Wind Speed Based on Height of Observation

The Penman-Monteith method requires wind speed measured at 2 meters above the ground surface for ET_o calculation. Because wind speeds measured from taller heights tend to be larger than 2-meter wind speeds, we need to adjust the observations using the following equation:

\[ u_2 = u_2 * \frac{4.87}{log(67.8 * observation\_height - 5.42)} \] Using our example: \[ u_2 = 2.7778 * \frac{4.87}{log(67.8 * 10 - 5.42)} = 2.0777 \] ETo provides the adjust_wind_speed() function, which does this calculation for the user.

Step 4: Calculate the Slope of the Saturation Vapor Pressure Curve (\(\Delta\))

The saturation vapor pressure of the air changes with temperature. The slope of this relationship represents how effectively moisture can be wicked from the plant or soil surface into the air. This in turn affects the rate of ET_o. \(\Delta\) is defined as follows:

\[ \Delta = \frac{4098 * [0.6108 * exp(\frac{17.27 * T_{mean}}{T_{mean} + 237.3})]}{(T_{mean} + 237.3)^2} \] Using our example, we get: \[ \Delta = \frac{4098 * [0.6108 * exp(\frac{17.27 * 16.9}{16.9 + 237.3})]}{(16.9 + 237.3)^2} = 0.122 \] ETo provides the calc_svp_slope() function, which does the calculation for the user.

Step 5: Calculate Atmospheric Pressure (\(P\))

Atmospheric pressure has an effect on how much ET can occur. At lower pressures ET increases slightly, due to less force being exerted on the earth’s surface⁵. In the absence of measured pressure data, it can be estimated using elevation:

\[ P = 101.3 * [\frac{293 - 0.0065*z}{293}]^{5.26} \] Using our example, we get: \[ P = 101.3 * [\frac{293 - 0.0065*100}{293}]^{5.26} = 100.124 \] ETo provides the calc_pressure() function, which does the calculation for the user.

Step 6: Calculate the Psychrometric Constant (\(\gamma\))

Zotarelli et al (2010) explain \(\gamma\) as relating “the partial pressure of water in air to the air temperature so that vapor pressure can be estimated using paired dry and wet thermometer bulb temperature readings. Another way to describe the psychrometric constant is the ratio of specific heat of moist air at constant pressure to latent heat of vaporization. It is calculated as:

\[ \gamma = \frac{C_p*P}{\epsilon * \lambda} = 0.000665 * P \] where:

\(\gamma\) is the psychrometric constant, kPa °C^-1.
\(C_p\) is the specific heat at a constant pressure, 1.013 10^-3, MJ kg^-1 °C^-1.
\(P\) is the atmospheric pressure, kPa.
\(\epsilon\) is the ratio of the molecular weight of water vapor to dry air, 0.622.
\(\lambda\) is the latent heat of vaporization, 2.45, MJ kg^-1

Using our example:

\[ \gamma = \frac{C_p*100.1235083}{\epsilon * \lambda} = 0.067 \] ETo provides the calc_psychrometric_constant() function, which does the calculation for the user.

Steps 7, 8 and 9 outlined in Zotarelli et al (2010) are skipped here, as they are all just sub-components of the main ET_o equation.

Step 10: Calculate the Mean Saturation Vapor Pressure

As described in step 4, the saturation vapor pressure is related to the air temperature. For a given temperature, the saturation vapor pressure can be calculated as:

\[ e_{(t)} = 0.6108 * exp[\frac{17.27 * T}{T + 237.3}] \]

Where:

\(e_{(t)}\) is the saturation vapor pressure for a given temperature (kPa).
\(T\) is the air temperature (°C).

If daily minimum and maximum temperature are available, the daily average saturation vapor pressure should be calculated as the mean of the \(e_{(T_{min})}\) and \(e_{(T_{max})}\). Otherwise, \(e_{(T_{avg})}\) can be used. Using our example, this means \(e_s\) is calculated as: \[ e_s = \frac{e_{(T_{min})} + e_{(T_{max})}}{2} \] Using our example, we get: \[ e_s = \frac{1.431 + 2.564}{2} = 1.997 \] ETo provides the calc_sat_vapor_pressure() function, which does the calculation for the user.

Step 11: Calculate the Actual Vapor Pressure (\(e_a\))

The actual vapor pressure can be calculated using \(e_s\) and daily \(RH_{min}\) and \(RH_{max}\) values: \[ e_a = \frac{e_{(T_{min})} * [\frac{RH_{max}}{100}] + e_{(T_{max})} * [\frac{RH_{min}}{100}]}{2} \] In the absence of daily \(RH_{min}\) and \(RH_{max}\) observations, \(e_a\) can be calculated using daily mean \(RH\): \[ e_a = \frac{RH_{mean}}{100} * e_s \] Since we have both \(RH_{min}\) and \(RH_{max}\), we can calculate \(e_a\) as: \[ e_a = \frac{1.431 * 0.84 + 3 * 0.63}{2} = 1.409 \] ETo provides the calc_act_vapor_pressure() function, which does the calculation for the user.

Step 12: Calculate the Inverse Relative Earth-Sun Distance and the Solar Declination

These variables, along with those calculated in the next few steps, are required to model the solar energy reaching earth’s atmosphere. The inverse relative distance between the earth and sun (\(d_r\)) and the solar declination (the angle of the sun in the sky; \(\delta\)) are given by: \[ d_r = 1 + 0.033 * cos[\frac{2 * \pi}{365} * J] \] \[ \delta = 0.409 * sin[\frac{2*\pi}{365} * J - 1.39] \] Where:

J is the julian day (day of the year) ranging between 1 and 365 (366 for leap years).

Using our example, we calculate \(d_r\) and \(\delta\) as: \[ d_r = 1 + 0.033 * cos[\frac{2 * \pi}{365} * 187] = 0.967099 \] \[ \delta = 0.409 * sin[\frac{2*\pi}{365} * 187 - 1.39] = 0.3954358 \] ETo provides the calc_solar_declination() and calc_inverse_relative_distance() functions, which do the calculations for the user.

Step 13: Convert Latitude (\(\phi\)) from Decimal Degrees to Radians

To be compatible with the variables we just created, the latitude must be in units of radians: \[ \phi [rad] = \frac{\pi}{180} * \phi [decimal\ degrees] \] In our example, this gives us: \[ \phi [rad] = \frac{\pi}{180} * 50.8 = 0.887 \] ETo provides the lat_to_radians() function, which does the calculation for the user.

Step 14: Calculate the Sunset Hour Angle (\(\omega_s\))

The sunset hour angle is calculated using the latitude and the declination: \[ \omega_s = arccos[-tan(\phi) * tan(\delta)] \] In our example: \[ \omega_s = arccos[-tan(0.887) * tan(0.395)] = 2.108 \] ETo provides the calc_sunset_hour_angle() function, which does the calculation for the user.

Step 15: Calculate Extraterrestrial Radiation (\(R_a\))

The extraterrestrial radiation is the amount of radiation that hits the earth’s atmosphere and varies with the time of the year and latitude. It is calculated as: \[ R_a = \frac{24 * 60}{\pi}*G_{sc}*d_r*[(\omega_s * sin(\phi)*sin(\delta)) + (cos(\phi) * cos(\delta) * sin(\omega_s))] \] Where:

\(R_a\) is the extraterrestrial radiation (MJ M^-2 day^-1)
\(G_{sc}\) is the solar constant of 0.082 MJ M^-2 min^-1

Using our example, this gives: \[ \begin{multline} R_a = \\ \frac{24 * 60}{\pi}*G_{sc}*0.967*[(2.108 * sin(0.887)*sin(0.395)) + (cos(0.887) * cos(0.395) * sin(2.108))] = \\ 41.088 \end{multline} \] ETo provides the calc_extraterrestrial_rad() function, which does the calculation for the user.

Step 16: Calculate Clear Sky Solar Radiation (\(R_{so}\))

The extraterrestrial radiation calculated in the previous step can be scaled to radiation received at the Earth’s surface under a clear sky:

\[ R_{so} = (0.75 + 2e10^{-5}*z)*R_a \] Where z is the elevation above sea level in meters. Using our example: \[ R_{so} = (0.75 + 2e10^{-5}*100)*41.088 = 30.898 \] ETo provides the calc_clear_sky_radiation() function, which does the calculation for the user.

Step 17: Calculate Net Shortwave Solar Radiation (\(R_{ns}\))

The net solar radiation is calculated using daily incoming shortwave radiation calculated in step 2: \[ R_{ns} = (1-\alpha)*R_s \] Where:

\(R_{ns}\) is the net incoming shortwave radiation (MJ M^-2 day^-1).
\(\alpha\) is the crop surface albedo. This defaults to 0.23, but can be adjusted with the reference argument.
\(R_s\) is the incoming shortwave solar radiation (MJ M^-2 day^-1)

In our example:

\[ R_{ns} = (1-0.23)*22.07 \]

ETo provides the calc_shortwave_raidation() function, which does the calculation for the user.

Step 18: Calculate Net Outgoing Longwave Solar Radiation (\(R_{nl}\))

The Stefan-Boltzmann law states that radiation leaving earth’s surface is proportional to temperature raised to the fourth power. As such, we define \(R_{nl}\) as follows: \[ R_{nl} = \sigma*[\frac{(T_{max}+273.16)^4 + (T_{min}+273.16)^4}{2}] * (0.34 - 0.14*\sqrt{e_a}) * (1.35 * \frac{R_s}{R_{so}}-0.35)) \] Where:

\(R_{nl}\) is the net outgoing longwave radiation (MJ M^-2 day^-1).
\(\sigma\) is the Stefan-Boltzmann constant (4.903 * 10^-9; MJ K^-4 M^-2 day^-1)

Using our example, we get:

\[ R_{nl} = \sigma*[\frac{(21.5+273.16)^4 + (12.3+273.16)^4}{2}] * (0.34 - 0.14*\sqrt{1.409}) * (1.35 * \frac{16.994}{30.898}-0.35)) \] \[ = 3.712 \]

Note that if only \(T_{mean}\) temperature data are available, \((T_{mean} + 273.16)^4\) can be used in place of the \(T_{min}\) and \(T_{max}\) in the bracketed portion of the equation. ETo provides the calc_longwave_radiation() function, which does the calculation for the user.

Step 19: Calculate Net Radiation (\(R_n\))

Net radiation is simply the difference betwen incoming shortwave radiation and outgoing longwave radiation:

\[ R_n = R_{ns} - R_{nl} \] where \(R_n\) has units of MJ M^-2 day^-1. Using our example: \[ R_n = 16.994 - 3.712 = 13.282 \] ETo provides the calc_net_radiation() function, which does the calculation for the user.

Step 20: Calculate ET_o

Now we have all the components needed to calculate ET_o. Recall that our equation for calculating ET_o is \[ ET_o = \frac{0.408 \Delta (R_n - G) + \gamma \frac{900}{T + 273} u_2 (e_s - e_a)}{\Delta + \gamma (1 + 0.34 u_2)} \] We have not calculated the soil heat flux (\(G\)), as it is assumed to be zero at a daily time step. Finally, plugging in the values from our example: \[ ET_o = \frac{0.408 * 0.122 * (13.282 - 0) + 0.067 * \frac{900}{16.9 + 273} * 2.078 * (1.997 - 1.409)}{0.122 + 0.067 * (1 + 0.34 * 2.078)} = 3.88 \] which gives \(ET_o\) in millimeters. Luckily, you don’t have to do all of these steps if you want to calculate \(ET_o\)! ETo provides the etr_penman_monteith() function that does all of these calculations under the hood.

ETo Definition

Calculation of Daily ETo