PyEt v1.3.1: a Python package for the estimation of potential evapotranspiration

Evapotranspiration (ET) is a crucial flux of the hydrological water balance, commonly estimated using (semi-)empirical formulas. The estimated flux may strongly depend on the formula used, adding uncertainty to the outcomes of environmental studies using ET. Climate change may cause additional uncertainty, as the ET estimated by each formula may respond differently to changes in meteorological input data. To include the effects of model uncertainty and climate change and facilitate the use of these formulas in a consistent, tested, and reproducible workflow, we present PyEt. PyEt is an open-source Python package for the estimation of daily potential evapotranspiration (PET) using available meteorological data. It allows the application of 20 different PET methods on both time series and gridded datasets. The majority of the implemented methods are benchmarked against literature values and tested with continuous integration to ensure the correctness of the implementation. This article provides an overview of PyEt's capabilities, including the estimation of PET with 20 PET methods for station and gridded data, a simple procedure for calibrating the empirical coefficients in the alternative PET methods, and estimation of PET under warming and elevated atmospheric CO₂ concentration. Further discussion on the advantages of using PyEt estimates as input for hydrological models, sensitivity and uncertainty analyses, and hindcasting and forecasting studies (especially in data-scarce regions) is provided.

1 Introduction

Evaporation – the process by which water is converted from its liquid to vapor phase – is a central component of the global hydrological cycle (Katul et al., 2012). Evaporation has far-reaching impacts on both human societies and ecosystems (Oki and Kanae, 2006; Fisher et al., 2011). In the remainder of this paper, the term evapotranspiration (ET) is used to refer to the total evaporation flux from soil and waterbodies (evaporation) and vegetated surfaces (transpiration; Allen et al., 1998; Dingman, 2015). Information about the magnitude of the ET flux is important across different geoscience disciplines: it assists in predicting irrigation demands and crop water requirements in agriculture, supports efficient water resource management, guides operational strategies in hydropower and meteorological studies, and plays a crucial role in ecological research and climate-change impact assessments. Given that climate change – through warming and elevated CO₂ concentrations – is set to alter evapotranspiration, the need for its accurate estimation is of paramount importance, as it affects our understanding and assessment of past, present, and potential future impacts on ecosystem functioning (Milly and Dunne, 2016; Yang et al., 2019; Caretta et al., 2022).

Evapotranspiration can hardly be measured directly (Wang and Dickinson, 2012; Jensen and Allen, 2016) and is therefore commonly estimated using (semi-)empirical formulas from other, more easily obtained meteorological variables such as temperature, wind speed, and radiation. Over time, dozens of methods have been proposed and applied. Each of these methods generally results in slightly different estimates of evapotranspiration, depending on the methods and data used (Oudin et al., 2005; McMahon et al., 2013; Xu and Singh, 2000, 2001; Lemaitre-Basset et al., 2022). Most of these formulas estimate either the reference crop evapotranspiration (ET₀), which is ET from a reference surface or crop that is not short of water (Allen et al., 1998), or the potential evapotranspiration (PET), which is the maximum rate of ET that would occur given a sufficient water supply (Xiang et al., 2020). Potential evapotranspiration is determined by meteorological conditions, whereas water availability determines if actual evapotranspiration occurs at its potential rate (Jensen and Allen, 2016). Differences in the potential evapotranspiration estimate may cascade through a modeling chain and ultimately impact the results of a study. For example, Prudhomme and Williamson (2013), Lemaitre-Basset et al. (2022), and Bormann (2010) showed that the method used affects the results from hydrological climate change impact studies. Similarly, the estimation of water demand for efficient crop and irrigation management depends on potential evapotranspiration, and it may thus be impacted by the methods used (Kumar et al., 2012).

To account for the structural uncertainty of the different PET models, it has been recommended to use multiple methods (Seiller and Anctil, 2016; Beven and Freer, 2001; Velázquez et al., 2013). Such an approach can help improve the understanding of the effect of model uncertainty on PET estimates in, for example, historical climate studies (Zhou et al., 2020; Dakhlaoui et al., 2020; Yang et al., 2019) and climate change impact studies (Bormann, 2010; Seiller and Anctil, 2016; Gharbia et al., 2018; Shi et al., 2020). Climate change impact studies often rely on climate projection data or global observational datasets, which are generally available in a gridded format (i.e., netCDF, GRIB), thus requiring tools that can efficiently process such data. Some of these datasets also contain only a limited set of observed or projected meteorological variables, requiring PET estimation methods that use fewer inputs. It may also be necessary to account for environmental variables that change over time and impact the evapotranspiration, such as vegetation changes and increases in atmospheric CO₂ concentrations (Fatichi et al., 2016; Ainsworth and Rogers, 2007; Vremec et al., 2022). Studies like those mentioned above require software programs that (1) include multiple PET estimation methods, (2) are flexible in adjusting input parameters (e.g., empirical coefficients, crop data, and meteorological inputs), and (3) are applicable to both time series and gridded data given the spatial nature of many of these studies.

Existing tools for calculating evapotranspiration, such as “Evapotranspiration” in R (Guo et al., 2016) and “PyETo” (Richards, 2019) and “pyfao56” in Python (Thorp, 2022), are primarily designed for station-based time series data. This limits their applicability with gridded datasets. While Peterson et al. (2020) extended “Evapotranspiration” to “AWAPer” to process gridded data, its use is limited to Australia. For the large community of geoscientists working with Python, the number of available PET methods from existing packages is limited (three for PyETo and one for pyfao56) compared to the 21 methods in the R package. This highlights a gap in the availability of a software for the estimation of multiple PET methods for both time series and gridded data, with the input parameter flexibility required for advanced studies on PET. Given the increasing need to understand and predict environmental changes accurately across the globe, the availability of such software is of paramount importance for the geoscience community.

Opportunities also exist to further align these tools with the FAIR standards of findability, accessibility, interoperability, and reusability for research software, which is of crucial importance for the credibility and reproducibility of scientific studies (Barker et al., 2022). This involves improving methodological testing through continuous integration, inclusion of additional alternative PET methods, and enabling more flexibility in adjusting internal empirical coefficients. Such enhancements not only adhere to best practices in software development but also broaden the scope and applicability of these tools in diverse geoscientific contexts. The refinement and development of evapotranspiration estimation tools that fully embrace the FAIR principles are therefore crucial steps toward advancing the field, ensuring more reliable and comprehensive research outcomes in the face of evolving scientific needs (Wood et al., 1998; DeJonge and Thorp, 2017).

In this paper we introduce PyEt, an open-source Python package for the estimation of potential evapotranspiration. The aim of PyEt is to provide researchers and practitioners with a wide variety of tested, documented, and flexible Python functions that support multiple PET methods for both station and gridded data. All methods have a common application programming interface, allowing users to easily test different PET models for their application and, if desired, address structural uncertainty and changing conditions. The majority of the implemented methods are benchmarked against literature values and tested with continuous integration to ensure the correctness of the implementation. Allowing different types of input data, PyEt is also applicable in regions with sparsely distributed measurement stations where standard meteorological data (e.g., wind, relative humidity) are often unavailable. The software is available under MIT license from the Python Package Index (PyPI) (Vremec and Collenteur, 2022) and developed as a community project on GitHub (http://www.github.com/pyet-org/PyEt, last access: 17 September 2024).

The remainder of this paper is structured as follows. In the next section, the software design, capabilities, and benchmarking tests are described. Section 3 introduces the software through four examples, showing potential future users how to apply PyEt in real-world applications. These examples concentrate on addressing practical problems commonly faced by geoscientists in their daily work. Section 4 discusses future potential applications of PyEt, and we detail how we think it can help the scientific community improve the estimation of potential evapotranspiration. Conclusions and future plans are outlined in Sect. 5.

2 PyEt Python package

2.1 Software design

The basic design principle for PyEt was to build a software that is intuitive and easy to use by novice users with little programming experience yet flexible enough to allow advanced users to perform more complex analyses. The software uses a modular design, with formulas shared by different PET methods implemented as a single function. This reduces the amount of code and makes it easier to maintain the software and implement new methods. All of the PET methods are intended to work with the minimum input data required by the PET models (e.g., radiation, temperature) but also allow more user input if such data are available and allowed by the PET method (e.g., humidity, surface resistance in the Penman–Monteith model). Utility functions are available to the user or are called internally to compute unavailable variables (e.g., solar radiation from latitude value). Moreover, the constants in the empirical PET formulas (e.g., the Stefan–Boltzmann constant) are function arguments with default values from the literature, which may also be changed by the user to adapt the empirical relationship to another region. Finally, the available methods should work for both station (1D) and gridded data (2D/3D).

PyEt is part of the wider Python ecosystem and depends on three widely used and well-developed Python packages from the scientific Python stack: NumPy (Harris et al., 2020), Pandas (McKinney, 201), and Xarray (Hoyer and Hamman, 2017). The input and output data of PyEt are formatted as time series data in Pandas.Series or Xarray.DataArrays, which allows using all the Pandas or Xarray functions on the data (Harris et al., 2020; McKinney, 201; Hoyer and Hamman, 2017). These functions include gap-filling and selection functions for interpolation, resampling, clustering, and many more. Being part of a wider ecosystem, users can leverage other Python packages for visualization (e.g., Matplotlib (Hunter, 2007), MetPy (May et al., 2022)) and optimization and uncertainty analyses (SciPy (Virtanen et al., 2020), SpotPy (Houska et al., 2015)).

The software is hosted and developed on the GitHub platform and distributed under MIT license through the Python Packaging Index (PyPI). Documentation and example applications are available on a dedicated ReadTheDocs website (http://pyet.readthedocs.io, last access: 17 September 2024). The documentation for individual methods is also directly available in Python from the documentation strings. Each release of PyEt is automatically stored in the Zenodo repository and assigned a Digital Object Identifier (DOI). As such, PyEt complies with many of the recommendations for good research software development as given in, for example, Hutton et al. (2016) and the FAIR4RS (FAIR for Research Software) principles (Barker et al., 2022). The scripts or the Jupyter notebooks used to apply PyEt improve the reproducibility and provide a transparent report of the entire calculation process (Kluyver et al., 2016).

Figure 1 Scatter plots showing estimated PET with PyEt against PET values estimated with the R package “Evapotranspiration” from Guo et al. (2016) and literature values from Allen et al. (1998), McMahon et al. (2013), and Schrödter (1985).

2.2 Implemented methods and benchmarking

There are 20 methods currently implemented in PyEt for the estimation of daily potential evapotranspiration. Apart from the Penman–Monteith method, which is considered the standard by the Food and Agriculture Organization (FAO) (Allen et al., 1998) and the World Meteorological Organization (WMO), multiple alternative methods are also available in PyEt. An overview of these methods and the required input data is provided in Table 1. Depending on the method, different (amounts of) input data are required to compute the potential evapotranspiration. It is often also possible to provide different input data to the same method (e.g., the average or the minimum and maximum daily temperatures) or even that some input data is optional (as described in the footnotes of Table 1). In the case of optional input data, utility functions are used internally to estimate that data. In the example of the Penman–Monteith method, solar radiation does not necessarily need to be provided by the user and can be estimated from the latitude and actual duration of sunshine hours instead.

Monteith (1965)Walter et al. (2000)Allen et al. (1998)Penman (1948)Priestley and Taylor (1972)Wright (1982)Thom and Oliver (1977)Blaney and others (1952)McMahon et al. (2013)Xu and Singh (2001)Schrödter (1985)Hamon (1963)Oudin et al. (2005)Ansorge and Beran (2019)Romanenko (1961)Xu and Singh (2001)Linacre (1977)Xu and Singh (2001)Haude (1955)Schiff (1975)Turc (1961)Xu and Singh (2000)Jensen and Haise (1963)Oudin et al. (2005)Jensen and Allen (2016)McGuinness and Bordne (1972)Xu and Singh (2000)Hargreaves and Samani (1982)McMahon et al. (2013)Jensen and Allen (2016)Jensen et al. (1990)Abtew (1996)Xu and Singh (2000)Makkink (1957)McMahon et al. (2013)Oudin et al. (2005)

Table 1 Data requirements for different PET or ET₀ models, the corresponding PyEt function, and if benchmarking of the method was performed. The references include both the original publications of the models and the papers from which the equations were taken: (1) McMahon et al. (2013), (2) Oudin et al. (2005), (3) Xu and Singh (2001), (4) Ansorge and Beran (2019), Rosenberry et al. (2004), (5) Schrödter (1985), (6) Schiff (1975), (7) Jensen and Allen (2016), (8) Xu and Singh (2000).

^a The corresponding literature for each method is provided in the Appendix B22.
^b T_max and T_min can also be provided.
^c RH_max and RH_min can also be provided.
^d If actual vapor pressure is provided, RH is not needed.
^e Input for radiation can be (1) net radiation, (2) solar radiation, or (3) sunshine hours. If it is (1), then latitude is not needed. If it is (1, 3) latitude and elevation are needed.
^f One must provide either the atmospheric pressure or elevation.
^g The PM method can be used to estimate potential crop evapotranspiration if leaf area index or crop height data are available.
^h The effect of CO₂ on stomatal resistance can be included using the formulation of Yang et al. (2019).
ⁱ If net radiation is provided, RH and latitude are not needed.
^j If method 2 is used, u₂, RH_min and sunshine hours are required.
^k Additional input of T_max and T_min or T_dew.
^l Input can be RH or actual vapor pressure.
^m If method 1 is used, latitude is needed instead of R_s.
ⁿ T_max and T_min are also needed.

The PyEt project is intended to be used by a wide community, and any errors in the code may have consequences for other studies applying PyEt to obtain PET estimates. Special attention was therefore paid to benchmark the available methods to published literature values and data from well-known research and meteorological institutes (Allen et al., 1998; McMahon et al., 2013; Schrödter, 1985; Walter et al., 2000). These benchmarks are also implemented in the continuous integration and tested using the unittest testing framework (unittest, 2022). This ensures that the benchmarks are satisfied each time the software is updated in the future. New methods added to PyEt will be required to be accompanied by the appropriate benchmark data and tests. Figure 1 shows the results for each benchmarked method, indicating that the PET estimates from all these methods are equal to the benchmark values (i.e., all values are on the 1 : 1 line). Despite our best efforts, we acknowledge here that four methods have not (yet) been benchmarked due to a lack of appropriate data.

In various sections of the paper, the performance of PET estimation methods are evaluated using three key performance metrics: model bias (mm d⁻¹), the coefficient of determination (–), and the Kling–Gupta Efficiency (KGE, –) (Gupta et al., 2009). These metrics enabled comparisons between benchmark PET values from literature and those estimated with PyEt, as well as between PET values derived from alternative models and the Penman–Monteith method, both before and after calibration. The Python implementations of the package SpotPy (Houska et al., 2015) were used for each performance metric, and the formulas of these implementations are detailed in Appendix B22.

2.3 Performance

The computational efficiency of various PET methods was assessed by examining the computation times in relation to the time series length and the number of cells in a Xarray dataset. Computation time was evaluated by running all models on a benchmark configuration (with time series of varying lengths) using a 12th Gen Intel Core i7-1255U processor with 10 cores and 12 logical processors. All Xarray.DataArrays cover a period of 1 year, whereas the spatial resolution changes. This comparison highlights the trade-offs between computational complexity and data size but also demonstrates the performance of the methods.

Figure 2 demonstrates that all PyEt methods maintain computational times below 1 s for time series data with lengths ranging from 10 to 10 000 d. For multidimensional data, the computation time does not exceed 10 s for larger Xarray.DataArrays up to 100 000 cells. Moreover, the results in Fig. 2 show that the problem scales well and proportionally does not take more time for larger datasets. Notably, the Penman–Monteith and Priestley–Taylor methods exhibited the largest processing times, whereas methods like Jensen–Haise, Turc, Makkink, and Romanenko were faster. Future improvements will aim to increase this efficiency, in particular to support faster calculations in large-scale global studies.

Figure 2 Computational efficiency of different PET methods. This figure shows the comparative processing times of different PET estimation methods regarding the length of the time series (a) and the size of the Xarray data (b).

3 Example use cases

Below, four use cases of PyEt are presented to illustrate how the software can be used. The first example shows how to efficiently compute different potential evapotranspiration estimates using 20 various methods for station data. This example also illustrates how to use PyEt in general. The second example illustrates how to provide 3D estimates of PET using three different methods and gridded Xarray data. The third example shows how to calibrate different PET methods to local conditions and use the calibrated formula for hindcasting. The fourth and final example illustrates a workflow to account for the effects of warming and elevated CO₂ in climate change impact studies. The source code for these and other examples can be found in a Zenodo repository related to this paper (Vremec and Collenteur, 2024a).

3.1 Example 1: estimation of PET from station data

In this example, potential evapotranspiration is estimated for the town of De Bilt in The Netherlands using data provided by the Royal Netherlands Meteorological Institute (KNMI). The reference method used by the KNMI for the estimation of potential evapotranspiration is the Makkink method, which is also implemented in PyEt. The PET computed with the Makkink method is compared to the PET values from all other methods in PyEt. Several steps are taken in a Python script to estimate PET. The code implementing these steps is shown in the code example below. PyEt provides a convenient way to compute the PET with all available methods, pyet.calculate_all().

1. Import the necessary Python packages.

   import pandas as pd
   import pyet

2. Load the meteorological data.

   meteo = pd.read_csv("meteo.csv",
     index_col=0, parse_dates=True)

3. Determine the necessary input data for the PET model.

   tmean, tmax, tmin, rh, rs, wind, \
   pet_knmi = (meteo[col] for col in
               meteo.columns)
   lat = 0.91  # latitude
   elev = 4  # elevation

4. Estimate the potential evapotranspiration with all methods or the method of choice.

   pet_df = pyet.calculate_all(tmean,
     wind, rs, elev, lat, tmax, tmin, rh)
   pet_mak = pyet.makkink(tmean, rs,
     elevation=elev)

5. Visualize and analyze the results.

   pet_df.plot()
   pet_df.boxplot()
   pet_df.cumsum().plot()

Figure 3 Potential evapotranspiration estimates for the year 2018 computed with all available PET methods plotted as (a) time series, (b) box plots, and (c) cumulative PET.

The results from this analysis are shown in Fig. 3. From these visualizations, it is clear that the potential evapotranspiration depends on the chosen method. This can accumulate up to a 35 % deviation of the estimated annual flux from the mean in this example. Such substantial differences between the estimated fluxes motivate the use of multiple methods (ensemble modeling) (Beven and Freer, 2001; Krueger et al., 2010; Shi et al., 2020; Oudin et al., 2005). This example showed how PyEt can be used efficiently for this task.

3.2 Example 2: estimate PET for gridded data

Gridded three-dimensional data (x, y, and t) obtained from satellites, radar imagery, or post-processed products are rapidly becoming widely available. More and more public datasets exist with global PET estimates at 0.1° resolution (e.g., Martens et al., 2017; Xie et al., 2022). PyEt also supports such gridded data, as illustrated here for the E-OBS gridded dataset (Cornes et al., 2018) for Europe. The application of PyEt on gridded datasets is illustrated for the FAO-56, Makkink, and Hargreaves method. Xarray.DataArrays are used as input data instead of Pandas.Series. PyEt methods will return the same data type, a Xarray.DataArray. The workflow is comparable to that in the first example, except that now the individual PET methods are used.

Figure 4 Daily PET estimates for Europe from 6 to 8 June 2018 using meteorological data obtained from the E-OBS dataset (Cornes et al., 2018).

The results for the three methods and three time steps are shown in Fig. 4. These again show that results may also differ spatially depending on the PET method. Looking more closely at Fig. 4, we can observe that the FAO-56 and Makkink methods do not compute PET in eastern parts of Europe. The data do not include relative humidity and solar radiation for these areas, and thus PET cannot be computed using the FAO-56 or Makkink method. If NaN (not-a-number) values are present in the required input data for a PyEt method, the method also returns a NaN value. The Hargreaves method, on the other hand, does not require solar radiation or relative humidity data. It can therefore be used to compute PET in the eastern parts of Europe. This example showed how PyEt can be applied to estimate PET using gridded data and demonstrated the benefits of using alternative PET methods when data such as radiation or relative humidity are missing.

3.3 Example 3: calibration of PET models

The available input data often does not suffice to compute potential evapotranspiration with the Penman–Monteith equation. This can be the case in data-scarce regions or time periods or when using historical data or data from climate models. In such cases, alternative PET methods can be calibrated to the estimates obtained from the Penman–Monteith equation for a period when sufficient data are available. The calibrated method can then be used to estimate PET in periods of data scarcity. As concluded by several authors (Jensen and Allen, 2016; Valipour, 2015; Yang et al., 2021; Dlouhá et al., 2021), calibration of alternative models is often crucial to ensure that the model fits the regional climate. In this example, it is shown how the calibration of temperature-based PET models affects the model uncertainty for studies focusing on current and past climates.

The approach is illustrated for the town of Graz, Austria, where the input data required for Penman–Monteith are only available from 2000 to 2021. Imagine, however, that for our study we also need potential evapotranspiration data for the period 1961 to 2021, but only temperature data are available (e.g., from the Spartacus temperature dataset (Hiebl and Frei, 2016)). Several steps are taken to calibrate the following five temperature-based methods: Oudin, Hargreaves, McGuinness–Bordne, Hamon, and Blaney–Criddle. First, the PET for the period 2000–2021 is computed using the Penman–Monteith equation. In the second step, the coefficients of the temperature-based PET equations are estimated by calibrating the estimated PET from temperature-based methods to the Penman–Monteith PET. Calibration is done by minimizing the sum of the squared residuals between these two PET estimates using SciPy's (Virtanen et al., 2020) least_squares method. In the third and final step, these calibrated coefficients are used to estimate the PET for the period 1961–2021.

Figure 5 shows the computed PET with the default (row 1) and calibrated coefficients (row 2). The model bias (mm d⁻¹) and the Kling–Gupta efficiency between simulated and observed (Penman–Monteith) PET show an improved model fit for all methods after calibration. The use of calibrated methods reduces the model bias, which is visually illustrated by the annual PET flux (composed of daily values) in the last column of Fig. 5. Using the Spartacus temperature dataset (Hiebl and Frei, 2016), PET can now be estimated back to 1961 using the calibrated alternative PET methods.

Figure 5 Density scatter plots comparing simulated and observed (FAO-56) PET for uncalibrated (row 1) and calibrated models (row 2). The performance of the calibrated models is evaluated using the model bias in mm d⁻¹ (Bias), the coefficient of determination (R²), and the Kling–Gupta efficiency (KGE). The last column shows the annual PET sums for the period 1961–2021 using uncalibrated (row 1) and calibrated models (row 2).

3.4 Example 4: the effect of CO₂ on future PET estimates

In this example, it is shown how to account for changing environmental conditions affecting the PET flux when modeling the effects of climate change. Under a warmer and CO₂-richer future (Caretta et al., 2022), potential evapotranspiration tends to increase with increasing temperature (and vapor pressure deficit). A reduction in PET is expected under elevated CO₂ due to an increased stomatal resistance (Field et al., 1995; Ainsworth and Rogers, 2007). The increase in CO₂ is still commonly ignored in PET models employed for climate change studies, although excluding its stomatal effect may lead to an overestimation of PET (Kingston et al., 2009; Milly and Dunne, 2016; Vremec et al., 2022; Riedel et al., 2023). The effect of temperature increases on PET can be easily modeled with all available PET methods, as temperature is an input for all methods. The CO₂ stomatal effect, however, can only be directly accounted for with the Penman–Monteith method (Liu et al., 2022). Using a CO₂-dependent stomatal resistance model implemented in PyEt (Yang et al., 2019), the effect of elevated CO₂ on stomatal resistance can be considered (see Eq. B2). When calculating PET with alternative methods, Kruijt et al. (2008) and Trnka et al. (2014) argued that an adjustment factor for the atmospheric CO₂ concentration ($f_{{\mathrm{CO}}_{\mathrm{2}}}$) can be used to account for the effect of elevated CO₂ concentrations on PET. The scaling factor can be obtained from literature values (Kruijt et al., 2008; Trnka et al., 2014). Alternatively, the factor can be calibrated using the Penman–Monteith equation together with the CO₂-dependent stomatal resistance model (Eq. B1) to match the local climate and vegetation:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\mathrm{PET}}_{{\mathrm{CO}}_{\mathrm{2}}}& =f_{{\mathrm{CO}}_{\mathrm{2}}}{\mathrm{PET}}_{\mathrm{300}}\\ & =(\mathrm{1}+S_{{\mathrm{PET}}_{{\mathrm{CO}}_{\mathrm{2}}}}({\mathrm{CO}}_{\mathrm{2}}-\mathrm{300}\left)\right){\mathrm{PET}}_{\mathrm{300}}\end{array},\end{array}$$

where $S_{{\mathrm{PET}}_{{\mathrm{CO}}_{\mathrm{2}}}}$ is the relative sensitivity of PET to CO₂, PET₃₀₀ is the computed Penman–Monteith estimate at 300 ppm [CO₂] (preindustrial concentration), and PET${}_{{\mathrm{CO}}_{\mathrm{2}}}$ is the computed Penman–Monteith estimate under elevated CO₂ concentration (Yang et al., 2019). Such relationships can be easily implemented in PyEt, and $f_{{\mathrm{CO}}_{\mathrm{2}}}$ can be obtained by calculating PET₃₀₀ and PET${}_{{\mathrm{CO}}_{\mathrm{2}}}$ with the Penman–Monteith equation (Eq. B1) at ambient and elevated CO₂ concentration, respectively.

Building on the previous example, the Graz study area served as a practical example to demonstrate the application of the calibrated models in assessing the impact of warming and elevated CO₂ concentration on PET based on the projected increase in temperature and CO₂ concentration from the representative concentration pathways (RCPs) (Van Vuuren et al., 2011). Daily PET was calculated for each RCP scenario (2.6, 4.5, 6.0, and 8.5) by adding the projected increase in temperature and CO₂ concentration to the existing data for 2020–2021. Figure 6 shows the increase in the average annual PET (aggregated from daily values) under warming and elevated CO₂ concentrations according to the RCP scenarios. In Fig. 6c, the effects of elevated CO₂ concentration on PET were neglected, and only increases in temperature were considered. Similar to Milly and Dunne (2016), Yang et al. (2019), and Vremec et al. (2022), this example shows that neglecting the effect of elevated CO₂ on PET (Fig. 6c) can lead to overestimation of PET under future conditions.

Figure 6 Projected increase in temperature (a) and atmospheric CO₂ concentration (b) under the RCP scenarios, calculated increase in the average annual PET with warming (c), and PET with warming and elevated CO₂ concentration (d). The uncertainty bounds represent the 5th–95th percentile of the PET model ensemble.

4 Discussion

4.1 Improved handling of PET in scientific studies

Evapotranspiration data from lysimeter or eddy correlation measurements (Pastorello et al., 2020) are rare and, if available at all, only locally available for relatively short time periods. Thus, there is a widespread need to estimate evapotranspiration from more readily available meteorological data using (semi-)empirical approaches. In general, these approaches follow three steps, as outlined for example by Allen et al. (1998). Firstly, the potential evapotranspiration of a reference surface (hence reference evapotranspiration) is estimated using meteorological data. Secondly, a crop coefficient may be applied to transform the reference evapotranspiration into the potential crop evapotranspiration. Thirdly, a soil–water balance approach is used to account for reduced actual evapotranspiration if the soil–water storage is depleted. PyEt is designed to perform the first two steps. It can be easily complemented by soil–water balance approaches to calculate actual evapotranspiration. Hydrological models, however, often use PET directly as input.

Rainfall–runoff models represent one type of hydrological model where PET is commonly used as an input, either as gridded data in distributive models or as a spatially aggregated values in lumped-parameter models. Some studies (e.g., Andréassian et al., 2004; Oudin et al., 2005; Sperna Weiland et al., 2012) found that PET had little impact on the performance of such models, and they thus advocated for the use of simplistic PET models. However, Jayathilake and Smith (2021) found that model performance was clearly sensitive to PET at sites with water-limited evapotranspiration. More importantly, the choice of the PET model has been shown to affect the results of hydrological projections in climate change impact assessments (Kay and Davies, 2008; Seiller and Anctil, 2016; Dallaire et al., 2021; Lemaitre-Basset et al., 2022). PET is expected to be even more influential in the assessment of groundwater recharge (e.g., Bakundukize et al., 2011) and crop water demands (e.g., Webber et al., 2016), which – compared to runoff – are more directly linked to evapotranspiration. Thus, the selection of appropriate PET models needs to account for the research context and variable of interest.

Bormann (2010) found that PET models that are based on the same or similar climate variables exhibit different sensitivity to observed climate change. This finding suggests that appropriate PET models need to be specifically selected for the given region of interest. Guo et al. (2017) provides pointers to examining which variables are likely to be the most important for a particular location. For more detailed insights into the application of individual PET methods across various climates, refer to Table A2 and the studies by Allen et al. (1998), McMahon et al. (2013), Jensen and Allen (2016), Yang et al. (2021), and Pimentel et al. (2023). The comparison of PET estimates for Europe shown in Fig. 4 illustrates the spatial variability of differences in PET estimates obtained from different methods; as can be seen, the magnitude and pattern of PET estimates are similar in some regions (e.g., Scandinavia) but differ more strongly in others (e.g., southeastern Europe).

As indicated above, the performance of PET models may vary depending on the region considered. Approaches that were found to be applicable in one region may perform less well in other regions. In this case, PET models can be calibrated to a reference dataset by adjustment of the coefficients in the model equation, as shown in the third example. The reference dataset can either be observed evapotranspiration (e.g., from lysimeters) or PET obtained from a model considered to be reliable. This has been illustrated by Example 3, where the coefficients of temperature-based models were adjusted to achieve the best fit to the Penman–Monteith model. This approach can also be used to obtain consistent spatial distributions of PET. As shown in Example 2 (Fig. 4), the limited data availability for eastern Europe did not allow the application of the FAO-56 or Makkink method, while sufficient data were available for the Hargreaves method. Thus, one may consider calibrating the latter to one of the former methods where these are applicable and only then applying it to obtain estimates for the entire region. For a more advanced calibration procedure, see, for example, Haslinger and Bartsch (2016).

Often the range of PET models that can potentially be employed is pre-determined by data availability. This may be the case if historical records of climate data are to be used for the PET estimation, for example, as many weather stations do not measure all climate variables included in the Penman–Monteith equation. However, this is also often the case in assessments of hydrological impacts of climate change if projected climate variables have high uncertainty. Lai et al. (2022), for example, concluded that the high uncertainty of wind speed projected in complex terrain may increase the uncertainty in PET, whereas air temperature and solar radiation have low uncertainty and thus should be the parameters preferred in the PET model. Given the climate variables for which data is available, Table 1 can be used to identify the PET models that come into consideration. However, it is advised to evaluate the assumptions and limitations of the individual methods regarding their applicability in the given case. Please refer to the comments in Table A2 and the references in Table 1 for this purpose. We generally recommend applying all models that have been identified as suitable (PET model ensemble), but the purpose and specific implementation of such a multi-model approach will depend on the research context. Example 4 (Sect. 3.4) illustrated how PET model ensembles can be used to include model uncertainties in PET projections under warming and elevated atmospheric CO₂ concentration. Since the latter effect is frequently excluded in hydrological projections, Milly and Dunne (2016) and Yang et al. (2019) advocate the inclusion of the effect of elevated CO₂ on stomatal resistance when estimating PET under warming and elevated atmospheric CO₂ concentrations.

To evaluate the performance of the estimated PET against observed values or other PET methods, performance metrics are employed. In this paper, the model bias and coefficient of determination (R²) are used to provide readers and practitioners with a clear and concise assessment of the overall deviation in flux and the percentage of the variability explained by the model (Onyutha, 2024). Additionally, the Kling–Gupta Efficiency (KGE) (Gupta et al., 2009) is used, offering a comprehensive view by incorporating correlation, variability, and bias into a single statistic. KGE is particularly valued in hydrological modeling, as it addresses some limitations of the R² by providing a differentiated perspective on model performance. The choice of the most appropriate performance metric is important to reduce model uncertainty, as it affects the judgment of the performance of a method. For further guidance on the evaluation of various performance metrics, readers are referred to Krause et al. (2005) and Onyutha (2024).

To improve reliability and efficiency in estimating PET, it is crucial to use a reproducible workflow. Scripts provide an efficient way to document the modeling process and are an important step towards full reproducibility. As shown in the examples, Jupyter Notebooks (Kluyver et al., 2016) provide a solution for publishing code, results, and explanations in a single document. As such, the presented package and its application in this paper are in line with the steps suggested by Hutton et al. (2016) to improve reproducibility in hydrological studies. To speed up adaptation of the methods and allow a faster transfer between research teams, formal procedures such as benchmarking (e.g., Maxwell et al., 2014) can help to ensure confidence in key complex codes.

4.2 Building the PyEt community and outlook

As a community project, the success of PyEt depends on the uptake from and interaction with the community. This, in turn, depends on the ease of use and the trust in the project. Emphasis was put on designing user-friendly and well-documented software, including various user examples, and extensive benchmark testing using continuous integration. Since the initial launch of PyEt, the package has already seen a good community uptake. Apart from applications of the software in projects related to the authors, which include estimating PET under conditions of warming and elevated CO₂ concentrations, assessing potential crop evapotranspiration, and providing inputs for hydrological models (e.g., Vremec et al., 2022; Forstner et al., 2022; Collenteur et al., 2023; Jemeljanova et al., 2023), PyEt has also been independently used by other researchers. These studies have used PyEt for various purposes: integrating PET estimates with machine learning models for enhanced analytical capabilities (Vaz et al., 2022; Kajári et al., 2024) and combining them with software for computing groundwater recharge and the water balance (Hassanzadeh et al., 2024). The software has also played a role in evaluating hydroclimatic changes and generating regional and global PET estimates (e.g., Tercini and Mello Júnior, 2023; Aguayo et al., 2024; Ha et al., 2024). The quick uptake of the software by the community confirms the need for this software.

The primary channel for communication with the PyEt community is GitHub, which provides several options for discussions, tracking code issues, and code development. Users are encouraged to ask questions in GitHub discussions and to report potential issues, suggest improvements, and feature requests via the GitHub issue tracker. As a community project, we plan to continue to improve the existing code and develop new capabilities based on feedback and with help from the community. An example of developments that are currently underway is the adaptation of the current methods to also work for hourly data, allowing the estimation of hourly PET. Other future work will focus on improvements in usability and the inclusion of other alternative methods.

5 Conclusions

This paper introduced PyEt, a Python package for the estimation of daily potential evapotranspiration (PET). The package enables the inclusion of model uncertainty and climate change in the PET estimation in a consistent, tested, and reproducible environment. With PyEt, PET can be estimated using 20 different methods with just a few lines of Python code. Unlike existing tools for PET calculation, which are designed for station-based time series, PyEt can also be applied to gridded (3D) datasets. This is of great practical relevance, particularly in climate impact studies, where gridded datasets are often used. The examples described in this paper illustrate how PyEt can be used in geoscientific studies to (1) facilitate the characterization of model uncertainty using a multimodel approach (model ensembles), (2) calibrate PET models and apply them in data-scarce regions and time periods, and (3) include the effects of warming and elevated atmospheric CO₂ concentrations. The use of Python scripts and Jupyter Notebooks ensures reproducibility and provides a transparent report of the PET computation process. We believe that PyEt will help improve the handling of PET and allow a more sophisticated and comprehensive consideration of PET in environmental studies, particularly those related to climate change.

Appendix A: PET methods – climate suitability and applications

Table A1 List of variables and symbols used in the paper and appendix.

Table A2 Overview of PET methods: climate suitability, applications, and limitations. ET₀ refers to reference crop or surface ET, PET_C refers to potential crop or surface ET, PET_OW refers to potential ET for open water, PET_SL refers to potential ET for shallow lakes, PET_RR refers to potential ET for rainfall–runoff modeling. Based on Allen et al. (1998), McMahon et al. (2013), Jensen and Allen (2016), Yang et al. (2021), and Pimentel et al. (2023).

Appendix B: PET methods – equations

B1 Penman–Monteith (pm)

Through the introduction of the Penman–Monteith equation by Monteith (1965), a broad applicability of PET estimation to different surfaces and vegetation types was achieved (Jensen and Allen, 2016). This was done by implementing the plant aerodynamic resistance (r_a) and the surface resistance (r_s) in the PET formula:

$$\begin{array}{}\text{(B1)}& \text{PET}={\displaystyle \frac{\mathrm{1}}{\mathit{\lambda }{\mathit{\rho }}_{\mathrm{w}}}}\left[{\displaystyle \frac{\mathrm{\Delta }(R_{\mathrm{n}}-G)+{\mathit{\rho }}_{\mathrm{a}}{c}_p{K}_{\mathrm{u}}(e_{\mathrm{s}}-e_{\mathrm{a}})/r_{\mathrm{a}}}{\mathrm{\Delta }+\mathit{\gamma }\left(\mathrm{1}+\frac{r_{\mathrm{s}}}{r_{\mathrm{a}}}\right)}}\right].\end{array}$$

Users of PyEt can include leaf or canopy cover measurements (leaf area index, LAI) to calculate surface resistance (r_s), thereby accounting for the effects of crop management and phenology on PET. A modified stomatal resistance model also allows for the inclusion of the sensitivity of the stomatal resistance (r_l) to the atmospheric CO₂ concentration (e.g., Yang et al., 2019; Vremec et al., 2022):

$$\begin{array}{}\text{(B2)}& r_{\mathrm{s}}={\displaystyle \frac{r_{\mathrm{l}}\left({\mathrm{CO}}_{\mathrm{2}}\right)}{\mathrm{0.5}\mathrm{LAI}}}={\displaystyle \frac{r_{r_{\mathrm{l}}-\mathrm{300}}\left\{\mathrm{1}+S_{r_{\mathrm{l}-{\mathrm{CO}}_{\mathrm{2}}}}({\mathrm{CO}}_{\mathrm{2}}-\mathrm{300})\right\}}{\mathrm{0.5}\mathrm{LAI}}},\end{array}$$

where $S_{r_{\mathrm{l}}-\left[{\mathrm{CO}}_{\mathrm{2}}\right]}$ [ppm⁻¹] is the relative sensitivity of r_l to Δ [CO₂] and $r_{r_{\mathrm{l}}-\mathrm{300}}$ [s m⁻¹] is the reference stomatal resistance when atmospheric CO₂ concentration is 300 ppm. The relative sensitivity of r_l to Δ [CO₂] represents the change in r_l per ppm increase in CO₂ concentration.

If measurements of crop height exist, these data can be used to calculate the aerodynamic resistance to vapor and heat transfer (r_a) to represent the effects of crop phenology on PET:

$$\begin{array}{}\text{(B3)}& r_{\mathrm{a}}={\displaystyle \frac{\text{ln}\left[\frac{z_{\mathrm{m}}-d}{z_{\mathrm{om}}}\right]\text{ln}\left[\frac{z_{\mathrm{h}}-d}{z_{\mathrm{oh}}}\right]}{k^{\mathrm{2}}{u}_z}},\end{array}$$

where z_m is the reference level at which the wind speed is measured; z_h is the height of the temperature and humidity measurements; k is the von Karman constant (=0.41); u_z is the measured wind speed (Allen et al., 1998); d is the zero plane displacement height, taken as 0.67h_c [m]; z_om is the roughness parameter for momentum (=0.123h_c) [m]; and z_oh is the roughness parameter for heat and water vapor (=0.1z_om) [m] (Jensen and Allen, 2016).

Free parameters in the Penman–Monteith equation, available for calibration, include the bulk stomatal r_l=100 (ranging between 40–150) and surface resistance (ranging between 50–200), with values for specific surfaces or crops found in Jensen and Allen (2016). Additionally, one can adjust the surface albedo (α = 0.23, ranges between 0.04 for water surfaces and 0.9 for snow) (Jensen and Allen, 2016) to estimate net shortwave radiation (Eq. 38 in Allen et al., 1998), the empirical coefficients for net longwave radiation a=1.35 and $b=-\mathrm{0.35}$ (Eq. 39 in Allen et al., 1998), or the empirical coefficients for clear-sky radiation as1=0.25 and bs1=0.5 (Eq. 36 in Allen et al., 1998). Optional meteorological inputs include G, T_max, T_min, RH_max, RH_min, p, and N.

B2 ASCE-PM (pm_asce)

The ASCE Penman–Monteith equation for is computed after Walter et al. (2000) (Eq. 1):

$$\begin{array}{}\text{(B4)}& \text{PET}={\displaystyle \frac{\mathrm{0.408}\mathrm{\Delta }(R_{\mathrm{n}}-G)+\mathit{\gamma }\frac{C_{\mathrm{n}}}{T+\mathrm{273}}{u}_{\mathrm{2}}(e_{\mathrm{s}}-e_{\mathrm{a}})}{\mathrm{\Delta }+\mathit{\gamma }\left(C_{\mathrm{d}}{u}_{\mathrm{2}}\right)}},\end{array}$$

where C_n=900 and C_d=0.34 are empirical coefficients for short reference vegetation (grass), while for tall reference vegetation (alfalfa) C_n=1600 and C_d=0.38 can be specified. The free parameters a, b, as1, bs1, and α are consistent with those specified for the Penman–Monteith method, while the optional meteorological inputs remain the same.

B3 FAO-56 (pm_fao56)

The FAO-56 Penman–Monteith equation for reference crop evapotranspiration is computed after Allen et al. (1998) (Eq. 6):

$$\begin{array}{}\text{(B5)}& \text{PET}={\displaystyle \frac{\mathrm{0.408}\mathrm{\Delta }(R_{\mathrm{n}}-G)+\mathit{\gamma }\frac{\mathrm{900}}{T+\mathrm{273}}{u}_{\mathrm{2}}(e_{\mathrm{s}}-e_{\mathrm{a}})}{\mathrm{\Delta }+\mathit{\gamma }(\mathrm{1}+\mathrm{0.34}{u}_{\mathrm{2}})}}.\end{array}$$

The free parameters a, b, as1, bs1, and α are consistent with those specified for the Penman–Monteith method, while the optional meteorological inputs remain the same.

B4 Penman (penman)

The Penman's PET formulation is computed after Penman (1948):

$$\begin{array}{}\text{(B6)}& \text{PET}={\displaystyle \frac{\mathrm{1}}{\mathit{\lambda }{\mathit{\rho }}_{\mathrm{w}}}}\left[{\displaystyle \frac{\mathrm{\Delta }(R_{\mathrm{n}}-G)+\mathit{\gamma }(e_{\mathrm{s}}-e_{\mathrm{a}})(a_{\mathrm{w}}+b_{\mathrm{w}}{u}_{\mathrm{2}})}{\mathrm{\Delta }+\mathit{\gamma }}}\right],\end{array}$$

where free parameters for the Penman's wind function include a_w=1 and b_w=0.537 (Valiantzas, 2006), while Penman (1948) suggested values of a_w=2.626 and b_w=1.381. The free parameters a, b, as1, bs1, and α are consistent with those specified for the Penman–Monteith method, while the optional meteorological inputs remain the same.

B5 Priestley–Taylor (priestley_taylor)

Priestley–Taylor's PET formulation is computed after Priestley and Taylor (1972):

$$\begin{array}{}\text{(B7)}& \text{PET}={\mathit{\alpha }}_L{\displaystyle \frac{\mathrm{\Delta }(R_{\mathrm{n}}-G)}{\mathit{\lambda }{\mathit{\rho }}_{\mathrm{w}}(\mathrm{\Delta }+\mathit{\gamma })}},\end{array}$$

where α_L=1.26 is an empirical coefficient. The free parameters a, b, as1, bs1, and α are consistent with those specified for the Penman–Monteith method, while the optional meteorological inputs remain the same.

B6 Kimberly–Penman (kimberly_penman)

The Kimberly–Penman equation (Wright, 1982) is computed after Oudin et al. (2005):

$$\begin{array}{}\text{(B8)}& \text{PET}={\displaystyle \frac{\mathrm{\Delta }(R_{\mathrm{n}}-G)+\mathit{\gamma }(e_{\mathrm{s}}-e_{\mathrm{a}})w}{\mathit{\lambda }{\mathit{\rho }}_{\mathrm{w}}(\mathrm{\Delta }+\mathit{\gamma })}},\end{array}$$

where $w=u_{\mathrm{2}}\left[\mathrm{0.4}+\mathrm{0.14}\text{exp}\left(-{\left(\frac{(j-\mathrm{173})}{\mathrm{58}}\right)}^{\mathrm{2}}\right)\right]+\left[\mathrm{0.605}+\mathrm{0.345}\text{exp}\left(-{\left(\frac{j-\mathrm{243}}{\mathrm{80}}\right)}^{\mathrm{2}}\right)\right]$.

The free parameters a, b, as1, bs1, and α are consistent with those specified for the Penman–Monteith method, while the optional meteorological inputs remain the same.

B7 Thom–Oliver (thom_oliver)

Thom–Oliver's PET formulation is computed (Thom and Oliver, 1977) the same way as in Oudin et al. (2005)

$$\begin{array}{}\text{(B9)}& \text{PET}={\displaystyle \frac{\mathrm{\Delta }(R_{\mathrm{n}}-G)+\mathrm{2.5}\mathit{\gamma }(e_{\mathrm{s}}-e_{\mathrm{a}})\cdot a_{\mathrm{w}}(\mathrm{1}+b_{\mathrm{w}}{u}_{\mathrm{2}})}{\mathit{\lambda }{\mathit{\rho }}_{\mathrm{w}}(\mathrm{\Delta }+\mathit{\gamma }(\mathrm{1}+\frac{r_{\mathrm{s}}}{r_{\mathrm{a}}}\left)\right)}},\end{array}$$

where a_w=2.6 and b_w=0.536. The free parameters a, b, as1, bs1, α, r_l, r_s, and $S_{r_{\mathrm{l}}-\left[{\mathrm{CO}}_{\mathrm{2}}\right]}$ are consistent with those specified for the Penman–Monteith method, while the optional meteorological inputs remain the same.

B8 Blaney–Criddle (blaney_criddle)

Three different approaches can be taken to estimate the Blaney–Criddle PET, depending on the selected method.

$$\begin{array}{}\text{(B10)}& \mathrm{\text{PET}}=\left\{\begin{array}{ll}a+b{P}_y(\mathrm{0.46}T+\mathrm{8.13})& \text{if }\mathit{\text{Method}}=\mathrm{0}\\ & \text{ (after }\text{),}\\ k\cdot P_y(\mathrm{0.457}T+\mathrm{8.128})& \text{if }\mathit{\text{Method}}=\mathrm{1}\\ & \text{ (after }\text{, Eq. 6),}\\ k\mathrm{1}+b_{\mathrm{var}}\cdot P_y(\mathrm{0.46}T+\mathrm{8.13})& \text{if }\mathit{\text{Method}}=\mathrm{2}\\ & \text{ (after }\text{,}\\ & \text{Eqs. S9.7 and S9.8)},\end{array}\right.\end{array}$$

Here, k=0.65, $a=-\mathrm{1.55}$, and b=0.96 are empirical coefficients, while $k\mathrm{1}=\mathrm{0.0043}\cdot {\mathrm{RH}}_{\min }-\frac{n}{N}-\mathrm{1.41}$ and $b_{\mathrm{var}}=\mathrm{0.81917}-\mathrm{0.0040922}\cdot {\mathrm{RH}}_{\min }+\mathrm{1.0705}\cdot \frac{n}{N}+\mathrm{0.065649}\cdot u_{\mathrm{2}}-\mathrm{0.0059684}\cdot {\mathrm{RH}}_{\min }\cdot \frac{n}{N}-\mathrm{0.0005967}\cdot {\mathrm{RH}}_{\min }\cdot u_{\mathrm{2}}$.

B9 Hamon (hamon)

The PET formulation after Hamon (1963), as used in Oudin et al. (2005), is as follows:

$$\begin{array}{}\text{(B11)}& \text{PET}=k\cdot {\left({\displaystyle \frac{N}{\mathrm{12}}}\right)}^{\mathrm{2}}\cdot \exp \left({\displaystyle \frac{T}{\mathrm{16}}}\right),\end{array}$$

where k=1 is a calibration coefficient.

B10 Romanenko (romanenko)

Romanenko's PET formulation (Romanenko, 1961), as used in Oudin et al. (2005), is as follows:

$$\begin{array}{}\text{(B12)}& \text{PET}=k{\left(\mathrm{1}+{\displaystyle \frac{T}{\mathrm{25}}}\right)}^{\mathrm{2}}\cdot \left(\mathrm{1}-{\displaystyle \frac{e_{\mathrm{a}}}{e_{\mathrm{s}}}}\right),\end{array}$$

where k=4.5 is an empirical coefficient (Oudin et al., 2005). Optional meteorological inputs include T_max, T_min, RH_max, and RH_min.

B11 Linacre (linacre)

Linacre's PET formula (Linacre, 1977), as used in Oudin et al. (2005), is as follows:

$$\begin{array}{}\text{(B13)}& \text{PET}={\displaystyle \frac{\frac{\mathrm{500}\cdot T_{\mathrm{m}}}{\mathrm{100}-{\text{lat}}_{\text{deg}}}+\mathrm{15}\cdot (T-T_{\text{dew}})}{\mathrm{80}-T}},\end{array}$$

where $T_{\mathrm{m}}=T+\mathrm{0.006}\cdot \text{Elev}$.

B12 Haude (haude)

Haude's PET formulation (Haude, 1955), as used in Schiff (1975), is computed as follows:

$$\begin{array}{}\text{(B14)}& \text{PET}=k\cdot FK\cdot (e_{\mathrm{0}}-e_{\mathrm{a}})\cdot \mathrm{10},\end{array}$$

where k=1 is a calibration coefficient and FK represents Haude's monthly coefficients, as adapted by Schiff (1975).

B13 Turc (turc)

The PET formula, as derived from Turc (1961) and used in McMahon et al. (2013) (Eqs. S9.10 and S9.11), is computed as follows:

$$\begin{array}{}\text{(B15)}& \text{PET}=k\cdot c\cdot {\displaystyle \frac{T}{T+\mathrm{15}}}(\mathrm{23.88}{R}_{\mathrm{s}}+\mathrm{50}),\end{array}$$

where k=0.013 is an empirical coefficient and c, dependent on the relative humidity (RH), is defined as follows:

$$\begin{array}{}\text{(B16)}& c=\left\{\begin{array}{ll}\mathrm{1}+\frac{\mathrm{50}-\text{rh}}{\mathrm{70}}& \text{if }\text{RH}<\mathrm{50},\\ \mathrm{1}& \text{otherwise}.\end{array}\right.\end{array}$$

B14 Jensen–Haise (jensen_haise)

The PET according to the Jensen–Haise model (Jensen and Haise, 1963), varies depending on the chosen method.

$$\begin{array}{}\text{(B17)}& \mathrm{\text{PET}}=\left\{\begin{array}{ll}k\frac{R_{\mathrm{s}}}{\mathit{\lambda }{\mathit{\rho }}_{\mathrm{w}}}(T-T_x)& \text{if }\mathit{\text{Method}}=\mathrm{0}\\ & \text{ (after }\text{)},\\ k\frac{R_{\mathrm{a}}T}{\mathit{\lambda }{\mathit{\rho }}_{\mathrm{w}}}& \text{if }\mathit{\text{Method}}=\mathrm{1}\\ & \text{ (after }\text{)},\end{array}\right.\end{array}$$

Here, cr=0.025 is an empirical coefficient, while $T_x=-\mathrm{3}$, as used in Jensen and Allen (2016).

B15 McGuinness–Bordne (mcguinness_bordne)

McGuinness–Bordne's PET equation (McGuinness and Bordne, 1972), as used in Oudin et al. (2005), is as follows:

$$\begin{array}{}\text{(B18)}& \text{PET}=k{\displaystyle \frac{R_{\mathrm{a}}(T+\mathrm{5})}{\mathit{\lambda }{\mathit{\rho }}_{\mathrm{w}}}},\end{array}$$

where k=0.0147 is an empirical coefficient as suggested by Xu and Singh (2000).

B16 Hargreaves (hargreaves)

The Hargreaves PET equation (Hargreaves and Samani, 1982) is computed as follows:

$$\begin{array}{}\text{(B19)}& \text{PET}=k{c}_{\mathrm{HS}}{\displaystyle \frac{R_{\mathrm{a}}}{\mathit{\lambda }{\mathit{\rho }}_{\mathrm{w}}}}\sqrt{T_{\text{max}}-T_{\text{min}}}(T+\mathrm{17.8}).\end{array}$$

B17 FAO-24 (fao_24)

The FAO-24 PET equation from Doorenbos (1977); Jensen et al. (1990), as used in Xu and Singh (2000) (Eqs. 11 and 12), is as follows:

$$\begin{array}{}\text{(B20)}& \text{PET}=a+{\displaystyle \frac{\mathrm{\Delta }}{\mathrm{\Delta }+\mathit{\gamma }}}{\displaystyle \frac{R_{\mathrm{s}}}{\mathit{\lambda }{\mathit{\rho }}_{\mathrm{w}}}}(\mathrm{1}-\mathit{\alpha })b,\end{array}$$

where $a=-\mathrm{0.3}$ and $b=\mathrm{1.066}-\mathrm{0.13}\cdot \frac{\text{RH}}{\mathrm{100}}+\mathrm{0.045}{u}_{\mathrm{2}}-\mathrm{0.02}\cdot \frac{\text{RH}}{\mathrm{100}}{u}_{\mathrm{2}}-\mathrm{0.315}\cdot (\frac{\text{RH}}{\mathrm{100}}{)}^{\mathrm{2}}-\mathrm{0.0011}{u}_{\mathrm{2}}$. Free parameters include the surface albedo α, consistent with those used in the Penman–Monteith method.

B18 Abtew (abtew)

Abtew's PET equation (Abtew, 1996), as used in Xu and Singh (2000) (Eq. 14), is as follows:

$$\begin{array}{}\text{(B21)}& \text{PET}=k{\displaystyle \frac{R_{\mathrm{s}}}{\mathit{\lambda }{\mathit{\rho }}_{\mathrm{w}}}},\end{array}$$

where k=0.53 is an empirical coefficient as suggested by Xu and Singh (2000).

B19 Makkink (makink)

Makkink's PET equation (Makkink, 1957) is as follows:

$$\begin{array}{}\text{(B22)}& \text{PET}={\displaystyle \frac{\mathrm{\Delta }}{\mathrm{\Delta }+\mathit{\gamma }}}{\displaystyle \frac{R_{\mathrm{s}}}{\mathit{\lambda }{\mathit{\rho }}_{\mathrm{w}}}},\end{array}$$

where k=0.65 is the empirical coefficient recommended by Hiemstra and Sluiter (2011) and ranging between (0.61–0.77) (Jensen and Allen, 2016).

B20 Makkink–KNMI (makink_knmi)

The Royal Netherlands Meteorological Institute (KNMI) employs a slightly modified version of the Makkink equation, tailored specifically for conditions in the Netherlands, as described by Hiemstra and Sluiter (2011):

$$\begin{array}{}\text{(B23)}& \text{PET}=k{\displaystyle \frac{\mathrm{\Delta }}{\mathrm{\Delta }+\mathit{\gamma }}}{\displaystyle \frac{R_{\mathrm{s}}}{\mathit{\lambda }{\mathit{\rho }}_{\mathrm{w}}}},\end{array}$$

where k=0.65. The calculations for s, e_s, γ, and λ are specified as follows:

• $\mathrm{\Delta }=\frac{\mathrm{7.5}\cdot \mathrm{237.3}}{(\mathrm{237.3}+T{)}^{\mathrm{2}}}\ln \left(\mathrm{10}\right)\cdot e_{\mathrm{s}}$,

• $e_{\mathrm{s}}=\mathrm{0.6107}\cdot {\mathrm{10}}^{\frac{\mathrm{7.5}T}{\mathrm{237.3}+T}}$,

• $\mathit{\gamma }=\mathrm{0.0646}+\mathrm{0.00006}T$,

• $\mathit{\lambda }=(\mathrm{2501}-\mathrm{2.375}T)\mathrm{1000}$.

B21 Oudin (oudin)

According to Oudin et al. (2005) (Eq. 2), the potential evapotranspiration (PET) can be expressed as follows:

$$\begin{array}{}\text{(B24)}& \text{PET}=\left\{\begin{array}{ll}{R}_{\mathrm{a}}\frac{(T+k_{\mathrm{2}})}{\mathit{\lambda }{\mathit{\rho }}_{\mathrm{w}}{k}_{\mathrm{1}}}& \text{if }T+k_{\mathrm{2}}>\mathrm{0}\\ \mathrm{0}& \text{otherwise,}\end{array}\right.\end{array}$$

where k₂=5 and k₁=100 (ranging between 75–100) are empirical coefficients recommended by Oudin et al. (2005).

B22 Performance metrics

He we provide mathematical definitions of the performance metrics used to evaluate the PET models discussed in the main text. The model bias (mm d⁻¹) is calculated as the average difference between the PET values estimated using PyEt (${\widehat{y}}_i$) and the reference PET values, which include (i) literature values presented in Fig. 1 and (ii) PET values computed using the Penman–Monteith method in Example 3 (y_i) over n time steps:

$$\begin{array}{}\text{(B25)}& \text{Bias}={\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n({\widehat{y}}_i-y_i).\end{array}$$

The coefficient of determination (R²) was computed as follows:

$$\begin{array}{}\text{(B26)}& R^{\mathrm{2}}=\mathrm{1}-{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n(y_i-{\widehat{y}}_i{)}^{\mathrm{2}}}{{\sum }_{i=\mathrm{1}}^n(y_i-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{2}}}},\end{array}$$

where $\stackrel{\mathrm{‾}}{y}$ is the mean of the reference data.

The Kling–Gupta efficiency (KGE) is defined as follows:

$$\begin{array}{}\text{(B27)}& \text{KGE}=\mathrm{1}-\sqrt{(r-\mathrm{1}{)}^{\mathrm{2}}+(\mathit{\alpha }-\mathrm{1}{)}^{\mathrm{2}}+(\mathit{\beta }-\mathrm{1}{)}^{\mathrm{2}}},\end{array}$$

where r is the correlation coefficient between the reference and estimated data, α is the ratio of the standard deviation of estimated data to that of reference data, and β is the ratio of the mean of estimated data to that of reference data (bias ratio).