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 the lower data requirements compared to other methods and especially from the ASCE-standardized reference evapotranspiration (former FAO-56), which is the most data demanding 10 method and is commonly used as 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 (Ep) using as benchmark the ASCEstandardized reference evapotranspiration method (Er). 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 15 using global gridded temperature and Er data of the period 1950-2000 at 30 arc-sec resolution (~1 km at equator) from freely available climate geodatabases. The correction coefficients were produced as partial weighted averages of monthly Er/Ep ratios by setting the ratios’ weight according to the monthly Er magnitude and by excluding colder months with monthly values of Er or Ep <45 mm 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. 20 The validation procedure showed that the corrected Thornthwaite formula Eps using local coefficients led to a reduction of RMSE from 37.2 to 30.0 mm m for monthly and from 388.8 to 174.8 mm y for annual step estimations compared to Ep using as benchmark the values of Er method. The corrected Eps and the original Ep Thornthwaite formulas were also evaluated by their use in Thornthwaite and UNEP (United Nations Environment Program) aridity indices using as benchmark the respective indices estimated by Er. The analysis was made using the validation data of the stations and the results showed that 25 the correction of Thornthwaite formula using local coefficients increased the accuracy of detecting identical aridity classes with Er 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 indices assessment with the minimum data requirements (i.e. temperature) for locations where climatic data are limited. The global 30 https://doi.org/10.5194/essd-2021-115 O pe n A cc es s Earth System Science Data D icu ssio n s Preprint. Discussion started: 9 July 2021 c © Author(s) 2021. CC BY 4.0 License.

grids of local correction coefficients for Thornthwaite formula produced in this study are archived in PANGAEA database and can be assessed using the following link: https://doi.pangaea.de/10.1594/PANGAEA.932638 (Aschonitis et al., 2021).

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 35 classification and drought indices (Weiß and Menzel, 2008;Wang and Dickinson, 2012;McHanon, 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 40 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), Allen et al. (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 45 (Thornthwaite, 1948;McCloud, 1955;Hamon 1961Hamon , 1963Baier and Robertson, 1965;Malmström, 1969;Hargreaves et al., 1982;Camargo et al., 1999;Droogers and Allen 2002;Pereira and Pruitt 2004;Oudin et al., 2004;Trajkovic 2005Trajkovic , 2007Trajkovic 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 describe properly the evaporative flux, while various studies have shown differences among the 50 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 indices 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 internal part of the respective Thornthwaite aridity/humidity index and it was calibrated based on measured monthly evapotranspiration from some well-watered grass-55 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;Castaneda and Rao, 2005;Trajkovic and Kolakovic, 2009), and 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;Castaneda and Rao, 2005;Zhang et al., 2008;Bakundukize et al., 2011;Yang et al., 2017). Indicative 60 modifications were proposed by Willmott et al. (1985) using an additional parametrization presented for mean monthly temperature above 26.5 o C and an adjustment for variable daylight and month lengths. Camargo et al. (1999)  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 65 adaptation of the Thornthwaite scheme to estimate the daily reference evapotranspiration on two contrasting environments in USA and Brazil. Castaneda and Rao (2005)  The last years, advanced interpolation techniques, climatic models and other methods have achieved to generate gridded datasets of various climatic parameters (Hijmans et al., 2005;Sheffield et al., 2006;Osborn and Jones, 2014;Harris et al., 75 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 that will better support all hydro-climatic applications and especially to support large scale applications of aridity 80 indices, which are highly prone to data limitations. The hypothesis that is tested in this work is that a local correction coefficient that integrates the local mean effect of wind speed, humidity and solar radiation can improve the performance of the original Thornthwaite formula and to convert it at the same time to a formula of reference evapotranspiration for short reference crop.

Data 85
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.
The derivation/calibration procedure was performed at 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 arc-sec spatial resolution (~1×1 km at the equator) (Fig.1a,b). The 90 second database is of Aschonitis et al. (2017), which provides gridded data of mean monthly reference evapotranspiration Er of the period 1950-2000 at five different resolutions (30 arc-sec, 2.5 arc-min, 5 arc-min, 10 arc-min and 0.5 deg) (Fig.1c). The method used for estimating Er is the ASCE-standardized method (former FAO-56), which estimates reference evapotranspiration for short clipped grass (Allen et al., 2005). The database of Er (Aschonitis et al., 2017)  temperature from the first database of Hijmans et al. (2005) at 30 arc-sec resolution and for this reason the two gridded 95 databases are compatible.

[FIGURE 1]
The validation procedure was performed using raw data of stations from three different databases. The first database is 100 the CIMIS database (California Irrigation Management System -CIMIS, http://www.cimis.water.ca.gov), 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 semi-arid/temperate coastal, plain, mountain environments. In total 60 stations ( Fig.2a) were used from 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). This database 105 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). This 110 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 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 climate Köppen-Geiger class or the high 115 altitude of their location). The total number of stations used in the study from the three databases is 525 and their full description is given in Table S1 of the Supplementary material.

Derivation and validation of Thornthwaite correction coefficients for short reference crop based on ASCEstandardized method 120
The monthly potential evapotranspiration Ep using the Thornthwaite (1948) method after its adjustment for variable daylight and month lengths (Willmott et al., 1985) is estimated as follows: (1) The benchmark method that was used for developing correction coefficients for the temperature-based method of Thornthwaite Ep is ASCE-standardized method (former FAO-56), which estimates reference evapotranspiration from short clipped grass, as follows (Allen, et al., 2005): where Er: the reference evapotranspiration (mm d -1 ), Δ: the slope of the saturation vapour pressure-temperature curve (kPa o C -1 ), Rn: the net radiation at the crop surface (MJ m -2 d -1 ), G: the soil heat flux density at the soil surface (MJ m -2 d -1 ), γ: the 130 psychrometric constant (kPa o C -1 ), u2: the wind speed at 2 m height above the soil surface (m s -1 ), es: the saturation vapour pressure (kPa), ea: the actual vapour pressure (kPa), Tmean: the mean daily air temperature ( o C), Cn and Cd: constants, which vary according to the time step and the reference crop type and describe the bulk surface resistance and aerodynamic roughness.
Eq.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 -1 where the constants Cn and Cd have the values 900 and 0.34, 135 respectively. (Allen et al., 2005). The use of Eq.3 in daily or monthly step for short reference crop is equivalent to 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 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 partial weighted annual correction coefficients for Priestley-Taylor and Hargreaves-Samani evapotranspiration methods. The procedure starts with the 140 derivation of the monthly coefficient cth,i for each month i based on Eq.5. Applying this procedure, twelve values of monthly cth,i are produced. The 12 monthly cth coefficients are then used to build mean annual coefficients. As it was mentioned in Aschonitis et al. (2017), the efficiency of mean annual correction coefficients is mainly associated to their ability to better = ( 1.514 = (6.75 • 10 −7 ) • 3 − (7.71 • 10 −5 ) • 2 + (1.79 • 10 −2 ) • + 0.492 monthly cth,i coefficients are estimated considering the participation weight of each month in the annual Er. Moreover, under cold conditions, the monthly coefficients cth,i may present unrealistic values that significantly affect the weighted averages. To solve this problem, threshold values for the monthly Ep,i and Er,i were used before the inclusion of their cth,i in the weighted average estimations. Preliminary analysis showed that when the mean monthly Ep,i and/or Er,i are below ~45 mm month -1 (~1.5 mm d -1 ), then unrealistic mean monthly cth,i values occur (as unrealistic values are considered those, which are at least one 150 order of magnitude larger or smaller from 1). Taking into account the above, the following procedure was performed in order to obtain a partial weighted average based on monthly cth,i values after excluding those months with Er and/or Ep ≤ 45 mm month -1 as follows: If , > 45 mm month −1 then , = 1 else = 0 If , > 45 mm month −1 then , = 1 else = 0 where Eps,i: the corrected temperature-based short reference crop evapotranspiration (mm month -1 ) of month i. 160 The above procedure was followed in order to calibrate the annual partial weighted average Cth (Eq.10) for every location on the globe based on mean monthly Er and Ep of 1950-2000 using: • the gridded mean monthly temperature data of Hijmans et al. (2005) that were further used to estimate the mean monthly The validation procedure with the data of the 525 stations was performed by comparing the mean monthly and the mean annual benchmark values of Er (Eq.4) versus the original Ep (Eq.1) and versus the corrected Eps Thornthwaite formula (Eq.11) taking into account the annual partial weighted average coefficients Cth at the location of each station. The validation was made separately for each database of stations (ECAD, AGBM, CIMIS) but also all together using the following five statistical 170 criteria:

Evaluating the use of correction coefficients in aridity indices based on stations data
The role of the new corrected formula of Thornthwaite (Eq.11) as internal parameter of aridity indices was also evaluated against the original method (Eq.1). For this purpose, the AIUNEP (UNEP, 1997) and AITH (Thornthwaite, 1948) aridity indices were used. The difference between the two indices is that AIUNEP does not consider seasonality. The two indices estimated 180 based on Er (Eq.4) were used as benchmark in order to compare the respective indices calculated with the original Thornthwaite Ep (Eq.1) and the corrected Eps (Eq.11) using the 525 stations data. The evaluation was performed: • by comparing the estimated aridity classes of 525 stations produced by the benchmark AIUNEP and AITH values using Er versus the classes of the two indices using Ep and Eps, respectively.
• by comparing the respective values of the indices using 1:1 plots and the statistical metrics of Eqs. 12-16. 185 The AIUNEP aridity index is the simpler method for hydroclimatic analysis and it is given by the following equation: The AITH aridity index is calculated as follows: where

Derivation and validation of the Cth correction coefficients
The global map of the Cth correction coefficient was developed following the procedure described in Section 2.2 and it is given in Fig.3. The validation of the derived Cth 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- The aridity classes of 525 stations given by the benchmark AIUNEP using Er were identical at 76% with the classes of the AIUNEP using Ep and 93% identical with the classes of the AIUNEP using Eps. Similarly, the aridity classes of 525 stations given by the benchmark AITH using Er were identical at 52% with the classes of the AITH using Ep and 58% identical with the classes 230 of the AITH using Eps. Eps showed better performance compared to Ep at correctly identifying the aridity classes in both indices. The lower percentages of success in the case of AITH for both Ep and Eps are due to the double number of classes of AITH in comparison to AIUNEP. Merging the B and A classes of AITH to one Humid class, as in the case of AIUNEP, the successful identical codes are raised to 63% for Ep and 76% for Eps. The 1:1 log-log plots of AIUNEP using Er versus the AIUNEP using Ep and Eps are given in Fig.5a,b, respectively, while the 235 same comparisons using AITH are given in Fig.6a,b. The visual inspection of Figs.5,6 clearly shows that Eps outperforms the Ep in the range of non-humid classes of both AIUNEP and AITH. In order to highlight this result, the statistical metrics (Eqs.12-16) were estimated after splitting the stations in two groups (non-humid and humid) based on the respective thresholds of humid classes of each index calculated using Er (Table 2). Table 2  On the other hand, the statistics showed that Ep showed better performance in both AIUNEP and AITH aridity indices for their respective humid classes. This result is of less importance since Eps showed better performance compared to Ep 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 AIUNEP, there is only one Humid class (AIUNEP>0.65) and thus there is no point to compare the performance of Ep and Eps from a statistical point of view since their values will always lead 245 to the same classification code/characterization (i.e. Humid). In the case of AITH>20, the same justification of AIUNEP could be used since the detailed division of five humid classes (B1, B2, B3, B4, A) provided by AITH was proposed for the alternative use of the index as "humidity index" (Thornthwaite, 1948).

Validity of the derived Cth for periods beyond the calibration period
The derivation of local Cth coefficients at 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 255 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 Cth coefficients are the following: • They are in the form of high-resolution grids (30 arc-sec, ~1 km at equator), which have been developed using interpolation techniques that include the effects of latitude, longitude and elevation. These grids allow to derive more representative Cth that the Cth values of the specific territories range between 0.9-1.1 for 1950-2000 (Fig.3), it is not only a verification of the Cth derivation methodology but also an additional indication of a generalized temporal stability of Cth.
In the case of Fig.S3b, there is a distinctly deviated Cth pair of values from the 1:1 line (point indicated by a red arow), which is associated to a specific station belonging to the Centro de Investigación Atmosférica de Izaña. This station is an 275 exceptional case since it is at the top of a mountain at 2371 m a.s.l. in Tenerife island. The derived Cth of this station from the grid of the period 1950-2000 is almost half (Cth value equal to 1.37) from the one estimated using stations' raw data (Cth 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 Cth value of the grid corresponds to an area of ~1 km while the specific position of the station is at a very unique position, which can be described as the most extreme position within this pixel. There are also 3 stations in Tenerife island at 280 lowland areas where the derived Cth values of 1950-2000 are in agreement with those estimated by the stations' raw data.

Scale and other effects on the accuracy of the derived Cth
The case of Izana station in 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 % 285 deviation of Eps from Er, it was observed a relative systematic deviation in some stations of CIMIS-California database, which are concentrated in 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. In the specific stations, the average of Cth values of the period 1950-2000 from the position of these stations was 1.85, while the average of Cth values using their raw data was estimated at 1.46. Apart from the large topographic variation, another reason for the Cth 290 differences in these stations could be the bias that has been removed by clearing extreme flagged wind values in the data of CIMIS database, which are probably associated to hurricane or other extreme events in this region. The coastline region of California is strongly affected by hurricanes and the higher wind speeds during such events may have not been removed by the Sheffield et al. (2006) wind grids that were used by Aschonitis et al. (2017)  An additional analysis based only on the stations of California was made to show that a wider regional mean value of Cth coefficient could also be an additional option, especially when the whole territory is described by local Cth coefficients that The mean monthly and mean annual Eps values of these stations were computed using Cth=1.66 for all of them and compared with the respective Er values estimated with stations' raw data (Fig.S4a,b, supplementary material). The results of Fig.S4a,b showed that even a regional average of Cth values for California can lead to better results of Eps compared to Ep as it was given for monthly and annual estimations in Figs.S1a and S2a, respectively.

Justifications about the methodology for deriving annual Cth correction coefficients based on partial weighted averages
The initial trials to derive annual correction coefficients Cth of this study were made using the average value of the twelve monthly cth,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  310 (supplementary material), which corresponds to a position close to Garda Lake in Italy (10.124 o E, 45.45 o N). According to Fig.S5a, the annual average of monthly cth,i values for this location is equal to 2.4 due to the extremely high values during winter and especially January. Using the 2.4 value as annual correction coefficient, the Eps,i value of July becomes equal to 338 mm, which is 203 mm larger from the respective Er,i value of July (Fig.S5a). The specific procedure for deriving annual Cth coefficients was rejected due to this problem. A second approach was to use the 12 pairs of monthly Er and Ep for each 315 position on the grid in order to perform regression analysis based on the form y = a·x without intercept based on the form of Er = Cth·Ep. An example of the specific procedure is given in Fig.S5b using the data of Fig.S5a, where the annual Cth value was found equal to 0.98. The specific procedure provides annual Cth values, which are always closer to the monthly coefficients of the warmer months since optimization algorithms try to minimize the total error, which is mainly originated by the months that show larger evapotranspiration values. Despite the fact that the specific procedure pays less attention to the monthly cth,i 320 values of colder months, it was considered acceptable since most of the hydroclimatic applications require higher accuracy to the larger evapotranspiration values rather to 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-325 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 total size of ~70GB. In order to solve this problem, the method of partial weighted 330 average (Eqs.5-10) developed by Aschonitis et al. (2017) was used, which provides similar results to the regression analysis of y = a·x but allows to perform calculations step by step with a conventional computer unit in GIS environment using large gridded databases. For the data of Fig.S5a, the partial weighted average method provided a Cth value equal to 0.99, which is almost equal to 0.98 of Fig.S5b. The method of partial weighted average 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 335 one by one using the threshold of 45 mm month -1 , 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 month -1 was derived empirically after analysing many datasets using monthly and mean monthly data. In the case of monthly data, a representative example is given in Figs show that the monthly cth,i values of months with Er,i < 45 mm month -1 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 The produced global database of local Cth coefficients of this study has been archived in PANGAEA and can be assessed using the following link: https://doi.pangaea.de/10.1594/PANGAEA.932638 (Aschonitis et al., 2021). The database is provided at 5 different resolutions (30 arc-sec, 2.5 arc-min, 5 arc-min, 10 arc-min, 0.5 deg). The coarser resolutions are provided in order to 355 cover the observed resolution range in the initial climatic data used for developing the published Er gridded data by Aschonitis et al. (2017) (e.g., the temperature data of Hijmans et al. (2005) were provided at 30 arc-sec resolution, while the solar radiation, humidity and wind speed data of Sheffield et al. (2006) were provided at 0.5 deg resolution and rescaled to 30 arc-sec using bilinear interpolation). The data of different resolutions can be used as a tool to assess uncertainties associated to temperature variation effects within different resolution pixels or to estimate average values of the coefficients for larger territories, which 360 have problems at coarse resolutions (e.g., coastlines or islands that do not exist in 0.5 degree resolution) taking into account the concept and concerns of Daly et al. (2006).

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 365 for developing the correction coefficients was based on partial weighted averages of their respective mean monthly values estimated as the monthly ratios between the benchmark ASCE-standardized Er method (former FAO-56) versus the original Thornthwaite Ep. The correction coefficients were produced as partial weighted averages of monthly Er/Ep ratios by setting the ratios' weight according to the monthly Er magnitude and by excluding colder months because the Er/Ep ratio becomes highly