Correcting Thornthwaite potential evapotranspiration using a global grid of local coefficients to support temperature-based estimations of reference evapotranspiration and aridity indices

Thornthwaite's formula is globally an optimum candidate for large-scale applications of potential evapotranspiration and aridity assessment at different climates and landscapes since it has lower data requirements compared to other methods and especially from the ASCE-standardized reference evapotranspiration (formerly FAO-56), which is the most data-demanding method and is commonly used as the benchmark method. The aim of the study is to develop a global database of local coefficients for correcting the formula of monthly Thornthwaite potential evapotranspiration (E_p) using as benchmark the ASCE-standardized reference evapotranspiration method (E_r). The validity of the database will be verified by testing the hypothesis that a local correction coefficient, which integrates the local mean effect of wind speed, humidity, and solar radiation, can improve the performance of the original Thornthwaite formula. The database of local correction coefficients was developed using global gridded temperature, rainfall, and E_r data of the period 1950–2000 at 30 arcsec resolution (∼ 1 km at Equator) from freely available climate geodatabases. The correction coefficients were produced as partial weighted averages of monthly $E_{\text{r}}/E_{\text{p}}$ ratios by setting the ratios' weight according to the monthly E_r magnitude and by excluding colder months with monthly values of E_r or E_p < 45 mm per month because their ratio becomes highly unstable for low temperatures. The validation of the correction coefficients was made using raw data from 525 stations of Europe; California, USA; and Australia including data up to 2020. The validation procedure showed that the corrected Thornthwaite formula E_ps using local coefficients led to a reduction of RMSE from 37.2 to 30.0 mm m⁻¹ for monthly step estimations and from 388.8 to 174.8 mm yr⁻¹ for annual step estimations compared to E_p using as a benchmark the values of the E_r method. The corrected E_ps and the original E_p Thornthwaite formulas were also evaluated by their use in Thornthwaite and UNEP (United Nations Environment Program) aridity indices using as a benchmark the respective indices estimated by E_r. The analysis was made using the validation data of the stations, and the results showed that the correction of the Thornthwaite formula using local coefficients increased the accuracy of detecting identical aridity classes with E_r from 63 % to 76 % for the case of Thornthwaite classification and from 76 % to 93 % for the case of UNEP classification. The performance of both aridity indices using the corrected formula was extremely improved in the case of non-humid classes. The global database of local correction factors can support applications of reference evapotranspiration and aridity index assessment with the minimum data requirements (i.e., temperature) for locations where climatic data are limited. The global grids of local correction coefficients for the Thornthwaite formula produced in this study are archived in the PANGAEA database and can be assessed using the following link: https://doi.org/10.1594/PANGAEA.932638 (Aschonitis et al., 2021).

1 Introduction

The assessment of potential or reference evapotranspiration is among the most important components for many hydro-climatic applications such as irrigation design and management, water balance assessment studies, and assessment of aridity classification and drought indices (Weiß and Menzel, 2008; Wang and Dickinson, 2012; McMahon et al., 2013; Aschonitis et al., 2017).

Such applications, and especially applications of aridity classification and drought indices (UNEP, 1997; Thornthwaite, 1948; Palmer, 1965; Holdridge, 1967; Beguería et al., 2014) that are usually employed at large scales, require estimations of potential or reference evapotranspiration of respective scale. The major problem in such applications is not only the limited availability of stations per se but also the limitation of many stations to provide data for a complete set of parameters (i.e., precipitation, temperature, solar radiation, wind speed, humidity). A complete set of climate parameters is prerequisite for accurate estimations of potential or reference evapotranspiration using integrated methods such as these of Penman (1948), Shuttleworth (1993), Allen et al. (1998, 2005), and others which are expressions of energy balance. Unfortunately, large-scale applications suffer from these limitations, and the common solution is to use temperature-based formulas (Thornthwaite, 1948; McCloud, 1955; Hamon, 1961, 1963; Baier and Robertson, 1965; Malmström, 1969; Hargreaves and Samani, 1982; Camargo et al., 1999; Droogers and Allen, 2002; Pereira and Pruitt, 2004; Oudin et al., 2005; Trajkovic, 2005, 2007; Trajkovic and Kolakovic, 2009a, b; Almorox et al., 2015; Aschonitis et al., 2017; Sanikhani et al., 2019; Quej et al., 2019; Trajkovic et al., 2020). However, extensive literature shows that temperature-based formulas are inherently of low performance because temperature cannot properly describe the evaporative flux, while various studies have shown differences among the Penman–Monteith-based and temperature-based potential evapotranspiration assessments such as the one of Thornthwaite (1948), which is the most popular in aridity and drought index applications (Sheffield et al., 2012; Dai, 2013; van der Schrier et al., 2013; Trenberth et al., 2014; Yuan and Quiring, 2014; Zhang et al., 2015; Asadi Zarch et al., 2015).

The formula of Thornthwaite (1948) was firstly proposed as the internal part of the respective Thornthwaite aridity–humidity index, and it was calibrated based on measured monthly evapotranspiration from some well-watered grass-covered lysimeters in the eastern and central USA (Willmott et al., 1985; Van Der Schrier et al., 2011). The specific formula overestimates the potential evapotranspiration in humid climates, and underestimates it in arid climates (Pereira and Pruitt, 2004; Castañeda and Rao, 2005; Trajkovic and Kolakovic, 2009a, b). Thus, a number of efforts have been made to amend the parameters or constants of the empirical formula to adapt it to various geographical zones (Jain and Sinai, 1985; Pereira and Pruitt, 2004; Castañeda and Rao, 2005; Zhang et al., 2008; Bakundukize et al., 2011; Yang et al., 2017). Indicative modifications were proposed by Willmott et al. (1985) using an additional parametrization presented for mean monthly temperature above 26.5 ^∘C and an adjustment for variable daylight and month lengths. Camargo et al. (1999) substituted the mean monthly temperature by another factor called effective temperature considering the amplitude between maximum and minimum temperature. Jain and Sinai (1985) modified the constant in the general formula based on the min–max range of the annual mean air temperature to calculate the evapotranspiration for semiarid conditions. Pereira and Pruitt (2004) proposed an adaptation of the Thornthwaite scheme to estimate the daily reference evapotranspiration on two contrasting environments in the USA and Brazil. Castañeda and Rao (2005) recalibrated the coefficient of the general formula based on estimations of potential evapotranspiration using the FAO Penman–Monteith method in southern California while Bautista et al. (2009) performed a similar procedure for stations located in a coastal semiarid climate and inland tropical subhumid climate in regions of Mexico. Zhang et al. (2008) used a modified formula to estimate the actual evapotranspiration in cropland, shrubland, and forest located in the subalpine region of southwestern China. Bakundukize et al. (2011) used two modifications of the original Thornthwaite method for groundwater recharge estimations in the inter-lacustrine zone of east Africa. Yang et al. (2017) presented a method to quantitatively identify the differences in the spatiotemporal variabilities of global drylands between the Thornthwaite and Penman–Monteith parameterizations. Trajkovic et al. (2019, 2020) provided successful corrections of the original formula based on the FAO Penman–Monteith method for stations located in Hungary, Serbia, Romania, Croatia, and Slovakia.

In recent years, advanced interpolation techniques, climatic models, and other methods have achieved gridded datasets of various climatic parameters (Hijmans et al., 2005; Sheffield et al., 2006; Osborn and Jones, 2014; Harris et al., 2014; Brinckmann et al., 2016; Liu et al., 2020), facilitating attempts to develop global maps of potential/reference evapotranspiration and to investigate the accuracy of formulas of reduced parameters versus benchmark methods at global scale (Droogers and Allen, 2002; Weiß and Menzel, 2008; Zomer et al., 2008; Aschonitis et al., 2017). A similar attempt is performed in this study, aiming to develop a global database of local correction coefficients for the original Thornthwaite formula. This attempt aims to support all hydro-climatic applications and specifically to support large-scale applications of aridity indices, which are highly affected by the use of different potential evapotranspiration methods (Proutsos et al., 2021). The hypothesis that is tested in this work is that a global grid of local correction coefficients that integrates the local mean effects of wind speed, humidity, and solar radiation can improve the performance of the original potential evapotranspiration formula of Thornthwaite by converting it into a formula of reference evapotranspiration for short reference crop based on the FAO Penman–Monteith concept.

2 Data and methods

2.1 Data

The methodological steps of the next sections are used to develop a global map of local coefficients for correcting the original potential evapotranspiration formula of Thornthwaite following a calibration and a validation procedure.

Figure 1 (a) Mean annual temperature for the period (Hijmans et al., 2005), (b) mean annual precipitation for the period (Hijmans et al., 2005), and (c) mean annual reference evapotranspiration of ASCE-standardized method for short reference crop for the period (Aschonitis et al., 2017) of 1950–2000.

The derivation–calibration procedure was performed at a global scale using global gridded data from two databases. The first database of Hijmans et al. (2005) provides gridded data of mean monthly precipitation P and mean monthly temperature T of the period 1950–2000 (WorldClim version 1.2) at 30 arcsec spatial resolution (∼ 1 × 1 km at the Equator) (Fig. 1a, b). The second database is of Aschonitis et al. (2017) and provides gridded data of mean monthly reference evapotranspiration E_r of the period 1950–2000 at five different resolutions (30 arcsec, 2.5 arcmin, 5 arcmin, 10 arcmin, and 0.5^∘) (Fig. 1c). The method used for estimating E_r is the ASCE-standardized method (formerly FAO-56 Penman–Monteith), which estimates reference evapotranspiration for short, clipped grass (Allen et al., 2005). The database of E_r (Aschonitis et al., 2017) was built using the temperature from the first database of Hijmans et al. (2005) at 30 arcsec resolution, and for this reason the two gridded databases are compatible.

The validation procedure was performed using raw data of stations from three different databases. The first database is the CIMIS database (California Irrigation Management System – CIMIS, https://cimis.water.ca.gov/stations.aspx, last access: 1 October 2020), which includes stations from California, USA, and it was selected because it provides a dense and descriptive network of stations for a specific region that combines semiarid/temperate coastal, plain, and mountain environments. In total 60 stations (Fig. 2a) were used from the CIMIS database that have at least 15 years of observations, with a significant part of their observations after 2000. The second database is the AGBM database (Australian Government – Bureau of Meteorology, http://www.bom.gov.au, last access: 1 June 2020). This database includes many stations from Australia and was selected because the station's network covers a large territory with a large variety of climate classes from desert to tropical climate. The selection of stations was performed in order to cover all the possible existing Köppen–Geiger climatic types (Peel et al., 2007) and altitude ranges that exist in the Australian territory. In total 80 stations were used (Fig. 2b) that have at least 15 years of observations, with a significant part of their observations after 2000. The third database is the ECAD database (European Climate Assessment & Database, https://www.ecad.eu, last access: 1 February 2020). This database is a network that contains more than 20 000 stations throughout Europe and provides daily observations of climatological parameters. In this study, a final number of 385 stations (Fig. 2c) was selected from this database because they contained complete data of precipitation, temperature, solar radiation, relative humidity, and wind speed for a period of at least 20 years with a significant part of their observations after 2000. Some additional stations from the three databases (CIMIS, AGBM, ECAD), which do not have at least 15 years of observations, were selected due to their special Köppen–Geiger climate class or the high altitude of their location. The total final number of stations used in the study from the three databases is 525, and their full description is given in Table S1 of the Supplement.

Figure 2 (a) The 60 stations in California from the CIMIS database, (b) 80 stations in Australia from the AGBM database, and (c) 385 stations in Europe from the ECAD database.

2.2 Derivation and validation of Thornthwaite correction coefficients for short reference crop based on ASCE-standardized method

The monthly potential evapotranspiration E_p using the Thornthwaite (1948) method after its adjustment for variable daylight and month lengths (Willmott et al., 1985) is estimated as follows.

$$\begin{array}{}\text{(1)}& E_{\text{p}}=\mathrm{16}\cdot {\left({\displaystyle \frac{\mathrm{10}\cdot T_{\text{mean}}}{J}}\right)}^a\cdot {\displaystyle \frac{N\cdot n}{\mathrm{365}}}\end{array}$$

$$\begin{array}{}\text{(2a)}& {\displaystyle }J={\sum }_{i=\mathrm{1}}^{\mathrm{12}}{j}_i\text{(2b)}& {\displaystyle }{j}_i={\left({\displaystyle \frac{T_{\text{mean},i}}{\mathrm{5}}}\right)}^{\mathrm{1.514}}\text{(2c)}& {\displaystyle }\begin{array}{rl}\mathit{\alpha }& =\left(\mathrm{6.75}\cdot {\mathrm{10}}^{-\mathrm{7}}\right)\cdot J^{\mathrm{3}}-\left(\mathrm{7.71}\cdot {\mathrm{10}}^{-\mathrm{5}}\right)\cdot J^{\mathrm{2}}\\ & +\left(\mathrm{1.79}\cdot {\mathrm{10}}^{-\mathrm{2}}\right)\cdot J+\mathrm{0.492}\end{array}\end{array}$$

$$\begin{array}{}\text{(3a)}& {\displaystyle }N={\displaystyle \frac{\mathrm{24}}{\mathit{\pi }}}\cdot {\mathit{\omega }}_{\text{s}}\text{(3b)}& {\displaystyle }{\mathit{\omega }}_{\text{s}}={\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}}-\arctan \left[{\displaystyle \frac{-\tan \left(\mathit{\phi }\right)\cdot \tan \left(\mathit{\delta }\right)}{X^{\mathrm{0.5}}}}\right]\end{array}$$

where $X=\mathrm{1}-\left[\tan \right(\mathit{\phi }){]}^{\mathrm{2}}\cdot [\tan \left(\mathit{\delta }\right){]}^{\mathrm{2}}$, if X≤0 then X=0.00001

$$\begin{array}{}\text{(3c)}& \mathit{\delta }=\mathrm{0.409}\cdot \sin \left(\mathrm{2}\cdot \mathit{\pi }\cdot {\displaystyle \frac{d_j}{\mathrm{365}}}-\mathrm{1.39}\right)\end{array}$$

Here E_p is the mean monthly potential evapotranspiration or potential evapotranspiration of month i (millimeters per month), T_mean,i is the mean monthly temperature (^∘C), n is the number of days in the month, N is the mean length of daylight of the days of the month (hours), J is the annual heat index, j_i is the monthly heat index, α is the function of the annual heat index, and d_j is the Julian day.

The benchmark method that was used for developing correction coefficients for the temperature-based method of Thornthwaite E_p is the ASCE-standardized method (formerly FAO-56 Penman–Monteith), which estimates reference evapotranspiration from short, clipped grass as follows (Allen et al., 2005):

$$\begin{array}{}\text{(4)}& E_{\text{r}}={\displaystyle \frac{\mathrm{0.408}\cdot \mathrm{\Delta }\cdot \left(R_{\text{n}}-G\right)+\frac{\mathit{\gamma }\cdot u_{\mathrm{2}}\cdot \left(e_{\text{s}}-e_{\text{a}}\right)\cdot C_n}{\left(T_{\text{mean}}+\mathrm{273.16}\right)}}{\mathrm{\Delta }+\mathit{\gamma }\cdot \left(\mathrm{1}+C_d\cdot u_{\mathrm{2}}\right)}},\end{array}$$

where E_r is the reference evapotranspiration (mm d⁻¹), Δ is the slope of the saturation vapor pressure–temperature curve (kPa ^∘C⁻¹), R_n is the net radiation at the crop surface (MJ m⁻² d⁻¹), G is the soil heat flux density at the soil surface (MJ m⁻² d⁻¹), γ is the psychrometric constant (kPa ^∘C⁻¹), u₂ is the wind speed at 2 m above the soil surface (m s⁻¹), e_s is the saturation vapor pressure (kPa), e_a is the actual vapor pressure (kPa),T_mean is the mean daily air temperature (^∘C), and C_n and C_d are constants, which vary according to the time step and the reference crop type and describe the bulk surface resistance and aerodynamic roughness. Equation (3) can be applied for two types of reference crop (i.e., short and tall). The short reference crop (ASCE-short) corresponds to clipped grass of 12 cm height and surface resistance of 70 s m⁻¹, where the constants C_n and C_d have the values 900 and 0.34, respectively (Allen et al., 2005). The use of Eq. (3) in daily or monthly steps for short reference crop is equivalent to the FAO-56 method (Allen et al., 1998), and this is how it is used in this study.

The derivation of a correction coefficient for Eq. (1) using as a benchmark the values of Eq. (3) is performed based on the same procedure proposed by Aschonitis et al. (2017) that has been used before for developing partially weighted annual correction coefficients for Priestley–Taylor and Hargreaves–Samani evapotranspiration methods. The procedure starts with the derivation of the monthly coefficient c_th,i for each month i based on Eq. (5). Applying this procedure, 12 values of monthly c_th,i are produced. The 12 monthly c_th coefficients are then used to build mean annual coefficients. As was mentioned in Aschonitis et al. (2017), the efficiency of mean annual correction coefficients is mainly associated with their ability to better describe the larger values of the dependent variable (i.e., the values of E_r during summer/hot months) and not the smaller values during cold periods when the absolute errors ($e_i=E_{\text{r},i}-E_{\text{p},i}$) are smaller. For this reason, weighted annual averages based on the monthly c_th,i coefficients are estimated considering the participation weight of each month in the annual E_r. Moreover, under cold conditions, the monthly coefficients c_th,i may present unrealistic values that significantly affect the weighted averages. To solve this problem, threshold values for the monthly E_p,i and E_r,i were used before the inclusion of their c_th,i in the weighted average estimations. Preliminary analysis showed that when the mean monthly E_p,i and/or E_r,i values are below ∼ 45 mm per month (∼ 1.5 mm d⁻¹), then unrealistic mean monthly c_th,i values occur (as unrealistic values are considered those that are at least 1 order of magnitude larger or smaller than 1). Taking into account the above, the following procedure was performed in order to obtain a partially weighted average based on monthly c_th,i values after excluding those months with E_r and/or E_p≤45 mm per month as follows.

$$\begin{array}{}\text{(5)}& {\displaystyle }{c}_{\text{th},i}=E_{\text{r},i}/E_{\text{p},i}\text{(6)}& {\displaystyle }\begin{array}{rl}& \text{If}{E}_{\text{r},i}>\mathrm{45}\text{mm per month, then }{F}_{\text{r},i}=\mathrm{1}\\ & \text{ otherwise}=\mathrm{0}\end{array}\text{(7)}& {\displaystyle }\begin{array}{rl}& \text{If}{E}_{\text{p},i}>\mathrm{45}\text{mm per month, then }{F}_{\text{m},i}=\mathrm{1}\\ & \text{ otherwise}=\mathrm{0}\end{array}\text{(8)}& {\displaystyle }{E}_{\text{r},i}^{\text{adj}}=E_{\text{r},i}\cdot F_{\text{r},i}\cdot F_{\text{m},i}\text{(9)}& {\displaystyle }{\mathrm{AE}}_{\mathrm{r}}^{\mathrm{adj}}={\sum }_{i=\mathrm{1}}^{\mathrm{12}}\left(E_{\text{r},i}^{\text{adj}}\right)\text{(10)}& {\displaystyle }{C}_{\text{th}}={\sum }_{i=\mathrm{1}}^{\mathrm{12}}\left({\displaystyle \frac{E_{\text{r},i}^{\text{adj}}}{{\mathrm{AE}}_{\mathrm{r}}^{\mathrm{adj}}}}\cdot c_{\text{th},i}\right)\end{array}$$

Here c_th,i is the monthly correction coefficient, F_r,i is the filter function for the reference method (ASCE) with values of 0 or 1, F_m,i is the filter function for the understudy model (Thornthwaite formula) with values of 0 or 1, $E_{\text{r},i}^{\text{adj}}$ is the adjusted monthly value of E_r,i from the ASCE-short method that becomes 0 when F_r,i or F_m,i is 0, ${\mathrm{AE}}_{\text{r}}^{\text{adj}}$ is the annual sum of the monthly $E_{\text{r},i}^{\text{adj}}$ adjusted values, C_th is the annual partially weighted average (p.w.a.) of the monthly c_th,i coefficients for short reference crop, and i is the index of each month. Considering the above, the final corrected Thornthwaite formula for monthly calculations is given by the following equation:

$$\begin{array}{}\text{(11)}& E_{\text{ps},i}=C_{\text{th}}\cdot E_{\text{p},i},\end{array}$$

where E_ps,i is the corrected temperature-based short reference crop evapotranspiration (millimeters per month) of month i.

The above procedure was followed to calibrate the annual partially weighted average C_th (Eq. 10) for every location on the globe based on mean monthly E_r and E_p of 1950–2000 using

• the gridded mean monthly temperature data of Hijmans et al. (2005) that were further used to estimate the original mean monthly gridded Thornthwaite E_p (Eq. 1) for the period 1950–2000 (in the form of 12 raster datasets of E_p for each month) and

• the respective mean monthly grids of E_r based on ASCE-standardized for short reference crop (Eq. 1) from Aschonitis et al. (2017) (in the form of 12 raster datasets of E_r for each month).

The validation procedure with the data of the 525 stations was performed by comparing the mean monthly and the mean annual benchmark values of E_r (Eq. 4) versus the original E_p (Eq. 1) and versus the corrected E_ps Thornthwaite formula (Eq. 11) considering the annual partially weighted average coefficients C_th at the location of each station. The validation was performed separately for each database of stations (ECAD, AGBM, CIMIS) but also all together using the following five statistical criteria.

$$\begin{array}{}\text{(12)}& {\displaystyle }\text{MAE}={\displaystyle \frac{\mathrm{1}}{N}}{\sum }_{i=\mathrm{1}}^N\left|S_i-O_i\right|\text{(13)}& {\displaystyle }\text{ME}={\displaystyle \frac{\mathrm{1}}{N}}{\sum }_{i=\mathrm{1}}^N\left(S_i-O_i\right)\text{(14)}& {\displaystyle }\text{RMSE}=\sqrt{{\displaystyle \frac{\mathrm{1}}{N}}{\sum }_{i=\mathrm{1}}^N{\left(S_i-O_i\right)}^{\mathrm{2}}}\text{(15)}& {\displaystyle }{R}_{\text{Sqr}}={\left[{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^N\left(O_i-{\stackrel{\mathrm{‾}}{O}}_i\right)\left(S_i-{\stackrel{\mathrm{‾}}{S}}_i\right)}{\sqrt{{\sum }_{i=\mathrm{1}}^N{\left(O_i-{\stackrel{\mathrm{‾}}{O}}_i\right)}^{\mathrm{2}}{\sum }_{i=\mathrm{1}}^N{\left(S_i-{\stackrel{\mathrm{‾}}{S}}_i\right)}^{\mathrm{2}}}}}\right]}^{\mathrm{2}}\text{(16)}& {\displaystyle }d=\mathrm{1}-{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^N{\left(S_i-O_i\right)}^{\mathrm{2}}}{{\sum }_{i=\mathrm{1}}^N{\left(\left|S_i-{\stackrel{\mathrm{‾}}{O}}_i\right|+\left|O_i-{\stackrel{\mathrm{‾}}{O}}_i\right|\right)}^{\mathrm{2}}}}\end{array}$$

Here MAE is the mean absolute error, ME the mean error, RMSE the root-mean-square error, R_Sqr the coefficient of determination, d the index of agreement, O the observed or benchmark value (i.e., E_r), S the value simulated by the model (i.e., E_p or E_ps), N the number of observations, and i the subscript referring to each observation. The value of perfect fit is 0 for the criteria MAE, ME, and RMSE while 1 is a perfect fit for the criteria R_Sqr and d. The values of the MAE, ME, and RMSE criteria have the same units as the observed and simulated data while R_Sqr and d are unitless.

2.3 Evaluating the use of correction coefficients in aridity indices based on station data

The role of the new corrected formula of Thornthwaite (Eq. 11) as an internal parameter of aridity indices was also evaluated against the original method (Eq. 1). For this purpose, the AI_UNEP (UNEP, 1997) and AI_TH (Thornthwaite, 1948) aridity indices were used. The difference between the two indices is that AI_UNEP does not consider seasonality. The two indices estimated based on E_r (Eq. 4) were used as the benchmark in order to compare the respective indices calculated with the original Thornthwaite E_p (Eq. 1) and the corrected E_ps (Eq. 11) using the 525 stations' data. The evaluation was performed

• by comparing the estimated aridity classes of 525 stations produced by the benchmark AI_UNEP and AI_TH values using E_r versus the classes of the two indices using E_p and E_ps, respectively, and

• by comparing the respective values of the indices using 1 : 1 plots and the statistical metrics of Eqs. (12)–(16).

The AI_UNEP aridity index is the simpler method for hydroclimatic analysis, and it is given by the following equation:

$$\begin{array}{}\text{(17)}& {\text{AI}}_{\text{UNEP}}={\displaystyle \frac{P_y}{E_y}},\end{array}$$

where P_y is mean annual precipitation (mm yr⁻¹) and E_y is mean annual potential evapotranspiration (mm yr⁻¹). The values of Eq. (16) are classified according to the following (UNEP, 1997; Cherlet et al., 2018):

• AI_UNEP < 0.05 → hyper-arid

• 0.03 ≤ AI_UNEP < 0.2 → arid

• 0.2 ≤ AI_UNEP < 0.5 → semiarid

• 0.5 ≤ AI_UNEP < 0.65 → dry subhumid

• 0.65 < AI_UNEP → humid

The classes for AI_UNEP > 0.65 are usually given as one humid class. The UNEP index does not consider the effect of seasonal variation in precipitation and potential evapotranspiration.

The AI_TH aridity index is calculated as follows:

$$\begin{array}{}\text{(18a)}& S=\sum _{i=\mathrm{1}}^{\mathrm{12}}\left(P_i-E_i\right)\end{array}$$

and

$$\begin{array}{}\text{(18b)}& D=\sum _{i=\mathrm{1}}^{\mathrm{12}}\left(E_i-P_i\right),\end{array}$$

$$\begin{array}{}\text{(19)}& {\text{AI}}_{\text{TH}}=\mathrm{100}{\displaystyle \frac{S-\mathrm{0.6}D}{E_y}},\end{array}$$

where P_i and E_i are the monthly precipitation and potential evapotranspiration of month i, respectively. S (mm yr⁻¹) considers only the positive values of (P_i−E_i) > 0, while (P_i−E_i) < 0 values are set to 0. In the case of D (mm yr⁻¹), only the positive values of (E_i−P_i) > 0 are considered while for (E_i−P_i) < 0 they are set to 0. The various climatic types according to AI_TH values are the following.

• −60 > AI_TH → hyper-arid (HE)

• −60 ≤ AI_TH < −40 → arid (E)

• −40 ≤ AI_TH < −20 → semiarid (D)

• −20 ≤ AI_TH < 0 → dry sub-humid (C1)

• 0 ≤ AI_TH < 20 → moist sub-humid (C2)

• 20 ≤ AI_TH < 40 → low humid (B1)

• 40 ≤ AI_TH < 60 → moderate humid (B2)

• 60 ≤ AI_TH < 80 → highly humid (B3)

• 80 ≤ AI_TH < 100 → very humid (B4)

• 100 ≤ AI_TH → hyper-humid (A)

3 Results

3.1 Derivation and validation of the C_th correction coefficients

The global map of the C_th correction coefficient was developed following the procedure described in Sect. 2.2, and it is given in Fig. 3. The validation of the derived C_th coefficients was performed for each one of the three datasets of stations (California, CIMIS; Australia, AGBM; Europe, ECAD), separately, by comparing the performance of mean monthly values (Fig. S1a–f, Supplement) and the performance of mean annual values (Fig. S2a–f, Supplement) of E_p (Eq. 1) and E_ps (Eq. 11) versus the benchmark values of E_r (Eq. 4). The statistical criteria (Eqs. 12–16) for both monthly and annual comparisons for each one of the three datasets of stations are given in Table 1. The respective monthly and annual comparisons after merging all the stations from the three datasets are also presented in Fig. 4a–d. From the results shown in Figs. S1, S2, and 4 and Table 1, a much better performance of E_ps compared to the original Thornthwaite formula E_p is observed in all cases, providing not only better monthly but also better annual reference evapotranspiration estimations that approximate the values of ASCE for short reference grass.

Figure 3 Global map of the annual partially weighted average C_th coefficients.

Table 1 Statistical metrics (Eqs. 12–16) for the comparisons between E_p vs. E_r and E_ps vs. E_r for CIMIS (California), AGBM (Australia), and ECAD (Europe) stations (the unit for MAE, ME, and RMSE is millimeters per month).

Figure 4 (a) The 1 : 1 plots of mean monthly E_p versus mean monthly E_r, (b) mean monthly E_ps versus mean monthly E_r, (c) mean annual E_p versus annual monthly E_r, and (d) mean annual E_ps versus annual monthly E_r, using the data of all 525 stations from the three databases CIMIS, AGBM, and ECAD.

3.2 Evaluating the use of C_th coefficient in AI_UNEP and AI_TH aridity indices

The use of C_th coefficients in AI_UNEP and AI_TH aridity indices was also evaluated based on the raw data of all 525 stations (California, CIMIS; Australia, AGBM; Europe, ECAD).

The aridity classes of 525 stations given by the benchmark AI_UNEP using E_r were 76 % identical with the classes of the AI_UNEP using E_p and 93 % identical with the classes of the AI_UNEP using E_ps. Similarly, the aridity classes of 525 stations given by the benchmark AI_TH using E_r were 52 % identical with the classes of the AI_TH using E_p and 58 % identical with the classes of the AI_TH using E_ps. E_ps showed better performance compared to E_p at correctly identifying the aridity classes in both indices. The lower percentages of success in the case of AI_TH for both E_p and E_ps are due to the double number of classes of AI_TH in comparison to AI_UNEP. Merging the B and A classes of AI_TH to one humid class, as in the case of AI_UNEP, the successful identical codes with E_r are raised to 63 % for E_p and 76 % for E_ps.

Figure 5 (a) The 1 : 1 log-log plots of AI_UNEP using mean monthly E_p versus AI_UNEP using mean monthly E_r and (b) AI_UNEP using mean monthly E_ps versus mean monthly E_r using the data of all 525 stations from the three databases CIMIS, AGBM, and ECAD.

Figure 6 (a) The 1 : 1 log-log plots of AI_TH using mean monthly E_p versus AI_TH using mean monthly E_r and (b) AI_TH using mean monthly E_ps versus mean monthly AI_TH using the data of all 525 stations from the three databases CIMIS, AGBM, and ECAD.

Table 2 Statistical metrics (Eqs. 12–16) for the comparisons between E_p vs. E_r and E_ps vs. E_r when they are applied in the (a) AI_UNEP and (b) AI_TH aridity indices by dividing the 525 stations into two groups based on non-humid or humid classes of each index (MAE, ME, and RMSE are unitless as the indices).

^∗ Estimated by the benchmark E_r.

The 1 : 1 log-log plots of AI_UNEP using E_r versus the AI_UNEP using E_p and E_ps are given in Fig. 5a, b, respectively, while the same comparisons using AI_TH are given in Fig. 6a, b. The visual inspection of Figs. 5 and 6 clearly shows that E_ps outperforms E_p in the range of non-humid classes of both AI_UNEP and AI_TH. To highlight this result, the statistical metrics (Eqs. 12–16) were estimated after splitting the stations into two groups (non-humid and humid) based on the respective thresholds of humid classes of each index calculated using E_r (Table 2). Table 2 verifies the better performance of E_ps compared to E_p in both AI_UNEP and AI_TH aridity indices for the non-humid classes.

On the other hand, the statistics showed that E_p showed better performance in both AI_UNEP and AI_TH aridity indices for their respective humid classes. This result is of less importance since E_ps showed better performance compared to E_p at correctly identifying the aridity classes in both indices based on all stations despite the fact that the stations belonging to humid classes were more in both indices (Table 2). Moreover, in the case of AI_UNEP, there is only one humid class (AI_UNEP > 0.65), and thus there is no point in comparing the performance of E_p and E_ps from a statistical point of view since their values will always lead to the same classification code/characterization (i.e., humid). In the case of AI_TH > 20, the same justification of AI_UNEP could be used since the detailed division of five humid classes (B1, B2, B3, B4, A) provided by AI_TH was proposed for the alternative use of the index as a “humidity index” (Thornthwaite, 1948).

4 Discussion

4.1 Validity of the derived C_th for periods beyond the calibration period

The derivation of local C_th coefficients at a global scale was performed using the mean monthly grid datasets of 1950–2000 assuming stationary climate conditions, while the validation was performed using stations' raw data from California and Australia that are expanded up to 2016 and stations' raw data from Europe that are expanded up to 2020 (Table S1). The reasons for choosing the specific grid datasets for the derivation of C_th coefficients are the following.

• They are in the form of high-resolution grids (30 arcsec, ∼ 1 km at Equator), which have been developed using interpolation techniques that include the effects of latitude, longitude, and elevation. These grids allow us to derive more representative C_th values for every position even when weather stations do not locally exist.

• They cover a large period of time (i.e., 1950–2000), so they can provide more representative mean annual p.w.a. C_th values. The upper threshold of the year 2000 of these grids also allows the validation dataset of stations to be more valid since the larger part of their data is after 2000, and this reduces the possibility of having been used in grids' development.

On the other hand, several works have shown climate differences after 2000 (Hansen et al., 2010; McVicar et al., 2012a, b; Wild et al., 2013; Willet et al., 2014; Sun et al., 2017). Such changes could possibly affect the validity of C_th coefficients and the final estimated values of E_r for periods beyond 2000. For this reason, the C_th values and the mean monthly E_r values of the grids of Aschonitis et al. (2017) of the period 1950–2000 were extracted from the positions of all 525 stations and compared with the respective values of computed E_r and C_th using stations' raw data, which go beyond 2000. The results of this comparison are given in Fig. S3a, b (Supplement) and clearly show that the gridded E_r data and C_th of 1950–2000 do not show serious deviations from their respective values for periods beyond 2000, allowing their safe use. Moreover, the fact that the original Thornthwaite (1948) formula was built before 1950 using data from the eastern and central USA and that the C_th values of the specific territories range between 0.9–1.1 for 1950–2000 (Fig. 3), it is not only a verification of the C_th derivation methodology but also an additional indication of a generalized temporal stability of C_th.

In the case of Fig. S3b, there is a distinctly deviated C_th pair of values from the 1 : 1 line (point indicated by a red arrow), which is associated with a specific station belonging to the Centro de Investigación Atmosférica de Izaña. This station is an exceptional case since it is at the top of a mountain at 2371 m a.s.l. on Tenerife (Canary Islands). The derived C_th of this station from the grid of the period 1950–2000 is almost half (C_th value equal to 1.37) of the one estimated using stations' raw data (C_th value equal to 2.44). This large difference is not the result of climate difference before and after 2000, but it is fully justified by the fact that the C_th value of the grid corresponds to an area of ∼ 1 km² while the specific position of the station is unique, which can be described as the most extreme position within this pixel. There are also three stations on Tenerife in lowland areas where the derived C_th values of 1950–2000 are in agreement with those estimated by the stations' raw data.

4.2 Scale and other effects on the accuracy of the derived C_th

The case of Izana station on Tenerife was the perfect example for triggering further investigation for the possible effects of scale in similar environments with extremely variable topography. Investigating the individual stations with the larger percentage of deviation of E_ps from E_r, a relative systematic deviation was observed at some stations of the CIMIS (California) database, which are concentrated at the coastline between Los Angeles and San Diego. The specific region is a narrow (∼ 20–30 km), highly urbanized coastal zone of ∼ 200 km, which is enclosed between the coastline and a hilly/mountainous zone. At the specific stations, the average of C_th values of the period 1950–2000 from the position of these stations was 1.85, while the average of C_th values using their raw data was estimated at 1.46. Apart from the large topographic variation, another reason for the C_th differences at these stations could be the bias that has been removed by clearing extreme flagged wind values in the data of the CIMIS database, which are probably associated with frequent extreme events in this region (extreme winds, droughts including wildfires, and heavy precipitation). This could justify the fact that the gridded C_th values of 1950–2000 at the positions of the stations are greater than the C_th values estimated by their raw data from CIMIS after removing flagged extreme values.

An additional analysis based only on the stations of California was made to show that a wider regional mean value of C_th coefficient could also be an additional option, especially when the whole territory is described by local C_th coefficients that are only > 1 or only < 1 (in California all local C_th coefficients are > 1). For this analysis, the average value of C_th=1.66 was estimated based on the values of local C_th coefficients of 1950–2000 from the locations of all stations of CIMIS (California). The mean monthly and mean annual E_ps values of these stations were computed using C_th=1.66 for all of them and compared with the respective E_r values estimated with stations' raw data (Fig. S4a, b, Supplement). The results of Fig. S4a, b showed that even a regional average of C_th values for California can lead to better results of E_ps compared to E_p as it was given for monthly and annual estimations in Figs. S1a and S2a, respectively.

4.3 Justifications about the methodology for deriving annual C_th correction coefficients based on partially weighted averages

The initial trials to derive annual correction coefficients C_th of this study were made using the average value of the 12 monthly c_th,i values of each i month. This procedure led to unreasonably high values due to the extreme high values during winter. An example of this problem based on the gridded data used in the calibration–derivation procedure is given in Fig. S5a (Supplement), which corresponds to a position close to Lake Garda in Italy (45.45^∘ N, 10.124^∘ E). According to Fig. S5a, the annual average of monthly c_th,i values for this location is equal to 2.4 due to the extremely high values during winter and especially during January. Using the 2.4 value as the annual correction coefficient, the E_ps,i value of July becomes equal to 338 mm, which is 203 mm larger than the respective E_r,i value of July (Fig. S5a). The specific procedure for deriving annual C_th coefficients was rejected due to this problem. A second approach was to use the 12 pairs of monthly E_r and E_p for each position on the grid in order to perform regression analysis based on the form $y=a\cdot x$ without intercept based on the form of $E_{\text{r}}=C_{\text{th}}\cdot E_{\text{p}}$. An example of the specific procedure is given in Fig. S5b using the data of Fig. S5a, where the annual C_th value was found equal to 0.98. The specific procedure provides annual C_th values, which are always closer to the monthly coefficients of the warmer months since optimization algorithms try to minimize the total error, which mainly originates in the months that show larger evapotranspiration values. Despite the fact that the specific procedure pays less attention to the monthly c_th,i values of colder months, it was considered acceptable since most of the hydroclimatic applications require higher accuracy for the larger evapotranspiration values rather than the lower ones.

A similar approach with the one of Fig. S5b was performed by Cristea et al. (2013) for deriving annual correction coefficients for the Priestley–Taylor method for 106 stations across the contiguous USA. The correction coefficients were estimated for each station by minimizing the sum of the squared residuals between Priestley–Taylor and the benchmark FAO-56, considering data only for the period April–September (warmer semester). The obtained optimized values of the correction coefficients for each station were then interpolated to produce a map of the Priestley–Taylor correction coefficients. For our study, the specific procedure was found to be extremely demanding in computing requirements since it was impossible to be performed pixel by pixel (777.6 million pixels) with a conventional computer unit for the whole globe using as input 24 rasters of extremely high resolution (∼ 1 km) with a total size of ∼ 70 GB. In order to solve this problem, the method using partially weighted averages (Eqs. 5–10) developed by Aschonitis et al. (2017) was used, which provides similar results to the regression analysis of $y=a\cdot x$ but allows us to perform calculations step by step with a conventional computer unit in the GIS environment using large gridded databases. For the data of Fig. S5a, the partially weighted average method provided a C_th value equal to 0.99, which is almost equal to 0.98 of Fig. S5b. The method of partially weighted averages is also extremely efficient since it is not restricted only to the warmer semester or to any other predefined period like the case of Cristea et al. (2013) since it controls all months one by one using the threshold of 45 mm per month, which is more appropriate for global applications and especially for applications of high-resolution data, giving the appropriate weight to the months with significant values of evapotranspiration.

The threshold of 45 mm per month was derived empirically after analyzing many datasets using monthly and mean monthly data. In the case of monthly data, a representative example is given in Fig. S6a, b (Supplement) using the monthly data of Embrun station in France (44.57^∘ N, 6.50^∘ E) 1980–2020. Figure S6a shows the box–whisker plots of monthly E_r,i values of the station, while Fig. S6b shows the respective box–whisker plots of monthly c_th,i values. The maximum c_th,i values of December, January, and February are outside the plot of Fig. S6b, with values of 30.1, 129.4, and 210.1, respectively. Figure S6a, b show that the monthly c_th,i values of months with E_r,i < 45 mm per month are extremely unstable, and their mean monthly value, even if it seems normal, cannot guarantee its safe use. In the case of mean monthly data, a representative example is given in Fig. 7, where the 6300 mean monthly c_th,i values derived by the raw data of the 525 stations were plotted against their respective mean monthly E_r values using a 2D density scatter plot. Figure 7 shows that the mean monthly c_th,i values of the stations start to exhibit extremely high dispersion below the threshold of 45 mm per month, with values reaching 1 order of magnitude larger than unity. In the case where there is a location where all months show E_r or E_p values below 45 mm per month, it is suggested to either use the non-zero C_th value of the closer location in the map of Fig. 3 or directly use the original Thornthwaite formula without correction.

Figure 7 The 2D scatter density plot between the 6300 mean monthly c_th,i values versus the respective mean monthly E_r values derived by the raw data of the 525 stations ($c_{\text{th},i}=\mathrm{0}$ or non-defined due to E_r and/or E_p=0 not being included in the graph).

5 Data availability

The produced global database of local C_th coefficients of this study has been archived in PANGAEA and can be assessed using the following link: https://doi.org/10.1594/PANGAEA.932638 (Aschonitis et al., 2021). The database is provided at five different resolutions (30 arcsec, 2.5 arcmin, 5 arcmin, 10 arcmin, 0.5^∘). The coarser resolutions are provided in order to cover the observed resolution range in the initial climatic data used for developing the published E_r gridded data by Aschonitis et al. (2017) (e.g., the temperature data of Hijmans et al., 2005, were provided at 30 arcsec resolution, while the solar radiation, humidity, and wind speed data of Sheffield et al., 2006, were provided at 0.5^∘ resolution and rescaled to 30 arcsec using bilinear interpolation). The data of different resolutions can be used as a tool to assess uncertainties associated with temperature variation effects within different resolution pixels or to estimate average values of the coefficients for larger territories, which have problems at coarse resolutions (e.g., coastlines or islands that do not exist in 0.5^∘ resolution), taking into account the concept and concerns of Daly (2006).

6 Conclusions

A global database of local correction coefficients for improving the performance of the monthly temperature-based Thornthwaite potential evapotranspiration method was built using gridded data covering the period 1950–2000. The method for developing the correction coefficients was based on partially weighted averages of their respective mean monthly values estimated as the monthly ratios between the benchmark ASCE-standardized E_r method (formerly FAO-56) versus the original Thornthwaite E_p. The correction coefficients were produced as partially weighted averages of monthly $E_{\text{r}}/E_{\text{p}}$ ratios by setting the ratios' weight according to the monthly E_r magnitude and by excluding colder months because the $E_{\text{r}}/E_{\text{p}}$ ratio becomes highly unstable for low temperatures. The correction coefficients were validated using raw data from 525 stations of California, Australia, and Europe that include independent data beyond 2000 up to 2020. The results showed that the correction coefficients significantly improved the monthly and annual results of the original Thornthwaite method E_p. The use of E_p with or without correction coefficients was also evaluated through their use in the aridity indices of Thornthwaite and UNEP versus the respective indices estimated based on the benchmark ASCE-standardized E_r. The results showed again that the correction coefficients significantly improved the performance of the indices compared to the original Thornthwaite method, especially in non-humid environments. The global database of local correction coefficients supports applications of reference evapotranspiration and aridity index assessment with minimum data requirements (i.e., mean temperature) for locations where climate data are limited. Uncertainties in the values of correction coefficients were observed in regions of high topographic variability, and a recommendation for such cases is the use of a regional average of correction coefficients or the use of local C_th values based on the available coarser resolutions provided in the database. The methods and results presented in this study and the observed uncertainties can be used as a base for future works focusing on (a) the validation of the correction coefficients for other places in the world, (b) comparison with other models of low data requirements, and (c) use of the p.w.a. method for recalibrating correction coefficients using station or climate models' data of recent periods.