Penman Monteith equation

meteorology
hydrology
soil science
crop science
Author

Dennis Walvoort

Published

January 16, 2023

Introduction

The Penman-Monteith equation is often used for modelling evapotranspiration of the soil-crop system. Penman-Monteith is a physics-based and comprehensive equation (Figure 1).

G T T T_v T v T->T_v Delta Δ T->Delta e0 e 0 T->e0 T_min T min e_a e a T_min->e_a R_nl R nl T_min->R_nl e_s e s T_min->e_s T_max T max T_max->e_a T_max->R_nl T_max->e_s T_0 T 0 P P T_0->P T_0->Delta T_0->e0 ETo ET o T_v->ETo e_a->T_v e_a->ETo e_a->R_nl R_s R s R_sr R sr R_s->R_sr R_ns R ns R_s->R_ns P->T_v gamma γ P->gamma P_0 P 0 P_0->P date date R_a R a date->R_a phi φ phi->R_a u u u2 u 2 u->u2 r_a r a u2->r_a z z z->P z->u2 R_so R so z->R_so g g g->P alpha α alpha->R_ns alpha_l α l alpha_l->P h h d d h->d z_om z om h->z_om LAI LAI h->LAI z_m z m z_m->u2 z_h z h z_h->r_a z_0 z 0 z_0->P RH_min RH min RH_min->e_a RH_max RH max RH_max->e_a Delta->ETo R_n R n R_n->ETo R_a->R_so R_so->R_sr R_sr->R_nl R_nl->R_n e_s->ETo gamma->ETo r_a->ETo r_s r s r_s->ETo rho_w ρ w rho_w->ETo lambda λ lambda->ETo lambda->gamma G G G->ETo G_sc G sc G_sc->R_a epsilon ε epsilon->T_v epsilon->ETo epsilon->gamma R_sp R sp R_sp->ETo R_sp->P R_ns->R_n R R R->R_sp M_v M v M_v->epsilon M_a M a M_a->epsilon M_a->R_sp N_a N a N_a->R k k k->R sigma σ k->sigma e0->Delta c_p c p c_p->gamma d->r_a z_om->r_a z_oh z oh z_om->z_oh karman κ karman->r_a r_l r l r_l->r_s c_1 c 1 c_1->LAI LAI_a LAI a LAI->LAI_a LAI_a->r_s sigma->R_nl ph ℎ rph ℏ ph->rph rph->sigma tau τ tau->rph c c c->sigma
Figure 1: Visualization of the Penman-Monteith equation. The notation follows Allen et al. (1998). Hover to see more details.

In this notebook, we will approximate the Penman-Monteith equation by a relatively simple Artificial Neural Network (ANN).

Data preparation

We use the Julia language for data analysis and modelling. First, we load some useful packages.

using PenmanMonteith
using Markdown
using CSV
using DataFrames
using Dates
using Unitful
using Flux
using Optimisers
using CairoMakie
using AlgebraOfGraphics
using Random

We initialize the pseudo random number generator to improve reproducibility.

Random.seed!(42);

Next, we calculate the reference crop evapotranspiration for meteorological station De Bilt (The Netherlands) by means of the Penman-Monteith equation. We use our own Julia package ‘PenmanMonteith’ for this. This package follows the implementation in Allen et al. (1998).

All daily observations needed to calculate reference crop evapotranspiration are downloaded from the KNMI-website.

d_meteo = CSV.read(
    "./data/etmgeg-260.csv", DataFrame,
    header=47, skipto=48,
    select=[:YYYYMMDD,:FG,:TN,:TX,:TG,:Q,:PG,:UX,:UN],
    normalizenames=true, missingstring=["     ", ""]);

Some conversions have to be done, before we can use the KNMI-data.

transform!(d_meteo,
    :YYYYMMDD => ByRow(x -> Date(string(x), dateformat"yyyymmdd")) => :date)

filter!(:date => x -> Date("1983-01-01") ≤ x ≤ Date("2022-12-31") , d_meteo)

select!(d_meteo,
    :date,
    :date => ByRow(x -> Date("1983-01-01") ≤ x ≤ Date("2012-12-31") ? "training" : "testing") => :type,
    :TN => ByRow(x -> uconvert(u"K", 0.1x*u"°C")) => :T_min,
    :TG => ByRow(x -> uconvert(u"K", 0.1x*u"°C")) => :T_avg,
    :TX => ByRow(x -> uconvert(u"K", 0.1x*u"°C")) => :T_max,
    :UN => ByRow(x -> 1.0x) => :RH_min,
    :UX => ByRow(x -> 1.0x) => :RH_max,
    :Q  => ByRow(x -> uconvert(u"MJ/m^2", 1.0x*u"J/cm^2")) => :Rₛ, 
    :FG => ByRow(x -> 0.1x*u"m/s") => :u,
    :PG => ByRow(x -> 0.01x*u"kPa") => :P);

Finally, we can compute the reference evapotranspiration as follows:

d_meteo.ETₒ = reference_evapotranspiration.(
    d_meteo.T_avg, # mean temperature
    d_meteo.T_min, # minimum temperature
    d_meteo.T_max, # maximum temperature
    actual_vapour_pressure.(d_meteo.T_min, d_meteo.T_max, d_meteo.RH_min, d_meteo.RH_max),
    d_meteo.Rₛ,  # global radiation
    d_meteo.P, # Daily mean sea level pressure
    d_meteo.date, # date
    52.1u"°", # latitude
    d_meteo.u, # daily mean windspeed 
    z = 1.90u"m", # altitude
    zₘ = 10.0u"m") # altitude wind measurements;

The resulting data set is given in Table 1.

Table 1: Meteorological data for meteorological station de Bilt (The Netherlands).
14610×11 DataFrame
14585 rows omitted
Row date type T_min T_avg T_max RH_min RH_max Rₛ u P ETₒ
Date String Quantity… Quantity… Quantity… Float64 Float64 Quantity… Quantity… Quantity… Quantity…
1 1983-01-01 training 271.25 K 273.95 K 276.25 K 76.0 99.0 2.74 MJ m^-2 5.1 m s^-1 102.3 kPa 0.275534 kg kPa^-1 s^-2
2 1983-01-02 training 275.85 K 276.95 K 277.75 K 91.0 100.0 1.45 MJ m^-2 3.6 m s^-1 102.55 kPa 0.209148 kg kPa^-1 s^-2
3 1983-01-03 training 276.55 K 281.85 K 284.15 K 89.0 99.0 0.27 MJ m^-2 6.7 m s^-1 101.59 kPa 0.527503 kg kPa^-1 s^-2
4 1983-01-04 training 277.05 K 281.75 K 284.85 K 65.0 97.0 0.89 MJ m^-2 7.2 m s^-1 101.1 kPa 1.11914 kg kPa^-1 s^-2
5 1983-01-05 training 277.95 K 282.95 K 286.05 K 82.0 97.0 0.96 MJ m^-2 6.2 m s^-1 101.24 kPa 0.6717 kg kPa^-1 s^-2
6 1983-01-06 training 278.15 K 283.45 K 286.05 K 66.0 93.0 1.99 MJ m^-2 8.2 m s^-1 101.44 kPa 1.14383 kg kPa^-1 s^-2
7 1983-01-07 training 276.25 K 279.35 K 280.95 K 71.0 87.0 2.69 MJ m^-2 5.7 m s^-1 102.24 kPa 0.65378 kg kPa^-1 s^-2
8 1983-01-08 training 275.25 K 278.35 K 280.75 K 69.0 95.0 3.48 MJ m^-2 4.1 m s^-1 102.91 kPa 0.370785 kg kPa^-1 s^-2
9 1983-01-09 training 276.55 K 279.35 K 282.15 K 75.0 96.0 0.69 MJ m^-2 5.7 m s^-1 102.84 kPa 0.751661 kg kPa^-1 s^-2
10 1983-01-10 training 275.65 K 279.95 K 282.35 K 86.0 97.0 2.26 MJ m^-2 3.6 m s^-1 103.17 kPa 0.277468 kg kPa^-1 s^-2
11 1983-01-11 training 279.65 K 281.35 K 283.15 K 82.0 96.0 2.04 MJ m^-2 3.6 m s^-1 103.16 kPa 0.395708 kg kPa^-1 s^-2
12 1983-01-12 training 277.25 K 279.05 K 280.65 K 80.0 89.0 0.67 MJ m^-2 5.1 m s^-1 102.7 kPa 0.697853 kg kPa^-1 s^-2
13 1983-01-13 training 273.75 K 277.35 K 280.05 K 74.0 95.0 0.57 MJ m^-2 5.7 m s^-1 101.49 kPa 0.730486 kg kPa^-1 s^-2
⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮
14599 2022-12-20 testing 281.45 K 283.15 K 284.35 K 75.0 97.0 0.55 MJ m^-2 4.6 m s^-1 101.0 kPa 0.743504 kg kPa^-1 s^-2
14600 2022-12-21 testing 279.65 K 281.05 K 282.35 K 85.0 98.0 1.37 MJ m^-2 3.3 m s^-1 100.83 kPa 0.336518 kg kPa^-1 s^-2
14601 2022-12-22 testing 280.05 K 281.85 K 283.25 K 83.0 99.0 2.14 MJ m^-2 2.9 m s^-1 100.34 kPa 0.20909 kg kPa^-1 s^-2
14602 2022-12-23 testing 281.55 K 282.65 K 284.25 K 81.0 98.0 0.43 MJ m^-2 3.1 m s^-1 100.27 kPa 0.571801 kg kPa^-1 s^-2
14603 2022-12-24 testing 280.05 K 282.25 K 284.15 K 78.0 95.0 2.39 MJ m^-2 3.9 m s^-1 101.12 kPa 0.360101 kg kPa^-1 s^-2
14604 2022-12-25 testing 280.75 K 281.75 K 282.45 K 94.0 98.0 0.55 MJ m^-2 2.1 m s^-1 101.07 kPa 0.355915 kg kPa^-1 s^-2
14605 2022-12-26 testing 277.05 K 280.15 K 282.95 K 66.0 96.0 1.96 MJ m^-2 4.6 m s^-1 101.41 kPa 0.654877 kg kPa^-1 s^-2
14606 2022-12-27 testing 277.05 K 278.85 K 280.45 K 78.0 89.0 2.64 MJ m^-2 4.7 m s^-1 102.11 kPa 0.393862 kg kPa^-1 s^-2
14607 2022-12-28 testing 279.15 K 282.35 K 284.55 K 81.0 93.0 0.65 MJ m^-2 7.1 m s^-1 100.63 kPa 0.789226 kg kPa^-1 s^-2
14608 2022-12-29 testing 279.65 K 282.35 K 284.85 K 75.0 88.0 1.28 MJ m^-2 6.9 m s^-1 100.22 kPa 0.939953 kg kPa^-1 s^-2
14609 2022-12-30 testing 277.75 K 281.65 K 286.05 K 78.0 95.0 0.95 MJ m^-2 6.2 m s^-1 100.29 kPa 0.796519 kg kPa^-1 s^-2
14610 2022-12-31 testing 283.55 K 287.15 K 289.05 K 58.0 97.0 0.73 MJ m^-2 6.6 m s^-1 100.4 kPa 1.46106 kg kPa^-1 s^-2

Modelling

We use the Flux-package to train an Artificial Neural Network (ANN). We first divide the data into a training set and a test set. Training data are from 1983 to 2012, test data from 2013 to 2022.

d_train = subset(d_meteo, :type => ByRow(x -> x == "training"))
d_test  = subset(d_meteo, :type => ByRow(x -> x == "testing"))

x_train = ustrip.(select(d_train, :T_min, :T_max, :RH_min, :RH_max, :Rₛ, :u))
x_test  = ustrip.(select(d_test,  :T_min, :T_max, :RH_min, :RH_max, :Rₛ, :u))
x_train = transpose(Matrix(x_train))
x_test  = transpose(Matrix(x_test))

y_train = ustrip.(d_train.ETₒ)
y_test  = ustrip.(d_test.ETₒ)
y_train .= max.(0.0, y_train)
y_test  .= max.(0.0, y_test);

x_train = Float32.(x_train)
y_train = Float32.(y_train)
x_test = Float32.(x_test)
y_test = Float32.(y_test);

We keep our ANN-model as simple as possible to prevent overfitting and to reduce computational demand. We only discern an input layer of six neurons, a single hidden layer of four neurons, and an output layer of one neuron. We also add a nonnegativity constraint to the output layer of our ANN to guarantee that evapotranspiration is always positive or zero. In Flux, this ANN is given by:

model = Chain(
    Dense(6 => 4),
    Dense(4 => 1),
    x -> abs2.(x)
)
Chain(
  Dense(6 => 4),                        # 28 parameters
  Dense(4 => 1),                        # 5 parameters
  var"#27#28"(),
)                   # Total: 4 arrays, 33 parameters, 388 bytes.

The ANN is visualized in Figure 2.

G cluster_1 input layer cluster_2 hidden layer cluster_3 output layer i1 h1 i1->h1 h2 i1->h2 h3 i1->h3 h4 i1->h4 i2 i2->h1 i2->h2 i2->h3 i2->h4 i3 i3->h1 i3->h2 i3->h3 i3->h4 i4 i4->h1 i4->h2 i4->h3 i4->h4 i5 i5->h1 i5->h2 i5->h3 i5->h4 i6 i6->h1 i6->h2 i6->h3 i6->h4 o1 h1->o1 h2->o1 h3->o1 h4->o1
Figure 2: visualization of the artificial neural network.

Next, we will define the loss-function (in our case, the Mean Square Error). The loss-function quantifies the difference between observations and predictions.

loss(ŷ, y) = Flux.Losses.mse(vec(ŷ), vec(y));

Next, we select and initialize an optimizer (in our case, ADAM). The aim of the optimizer is to reduce the value of the loss-function. In other words, by minimizing the loss-function the difference between the observations and the predictions will be reduced leading to a better approximation.

optimizer = Flux.setup(ADAM(), model);

Now we are ready to train our model. In Flux this can be achieved by running the following code:

for epoch ∈ 1:100_000
    Flux.train!(model, [(x_train,  y_train)], optimizer) do m, x, y
        loss(m(x), y)
    end
end

Results

The resulting value of the loss function for the training period is 0.064 (= 0.25 mm/d). More interestingly however, the value of the loss function for the testing period is 0.067 (= 0.26 mm/d).

Figure 3 and Figure 4 give scatter plots of the reference evapotranspiration computed by means of Penman Monteith and by the ANN for the training period and testing period. We see good agreement.

Figure 3: Reference evapotranspiration by means of Penman Monteith and by the ANN for the training period (1983-2012).
Figure 4: Reference evapotranspiration by means of Penman Monteith and by the ANN for the testing period (2013-2022).

Figure 5 gives time-series for the reference evapotranspiration estimated by means of Penman Monteith (PM) and by means of the ANN for the training period. Figure 6 gives these results for the testing period.

Figure 5: Evapotranspiration by means of Penman Monteith (PM) and by means of the ANN for the training period (1983-2012).
Figure 6: Evapotranspiration by means of Penman Monteith (PM) and by means of the ANN for the testing period (2013-2022).

Conclusion

It is possible to approximate the relatively complicated physics-based Penman-Monteith equation by a relatively simple and fast artificial neural network (ANN). Although the FAO (Allen et al. 1998) recommends the Penman-Monteith equation, and it is physics-based, it is certainly not without errors (Widmoser 2009). It is therefore not necessary and even undesirable to aim for a ‘perfect’ fit. The ANN may be used as top boundary condition of Richardson-Richards’ equation.

References

Allen, R. G., L. S. Pereira, D. Raes, and M. Smith. 1998. Crop Evapotranspiration - Guidelines for Computing Crop Water Requirements. FAO Irrigation and Drainage Paper 56. Rome, Italy: Food; Agriculture Organization of the United Nations. https://www.fao.org/3/X0490E/x0490e00.htm.
Widmoser, P. 2009. “A Discussion on and Alternative to the Penman–Monteith Equation.” Agricultural Water Management 96 (4): 711–21. https://doi.org/https://doi.org/10.1016/j.agwat.2008.10.003.

Reuse