A Comparative Study of Potential Evapotranspiration Estimation by Three Methods with FAO Penman–Monteith Method across Sri Lanka

: Among numerous methods that have been developed to estimate potential evapotranspiration (PET), the Food and Agricultural Organization Penman–Monteith model (FAO P–M) is often recognized as a standard method to estimate PET. This study was conducted to evaluate the applicability of three other PET estimation methods, i.e., Shuttleworth–Wallace (S–W) model, Thornthwaite (TW) and pan methods, to estimate PET across Sri Lanka with respect to the FAO P–M model. The meteorological data, i.e., temperature, relative humidity, wind speed, net solar radiation, and pan evaporation, recorded at 14 meteorologic stations, representing all climate and topographic zones of Sri Lanka, were obtained from 2009 to 2019. The models’ performances were assessed based on three statistical indicators: root mean squared error (RMSE), bias, and Pearson correlation coefﬁcient (R). In comparison with the FAO P–M model estimates, the seasonal and annual estimates of all three models show great differences. The results suggested that pan and S–W methods perform better in the dry zone of the country. Both S–W and pan methods underestimated PET over the entire county in all seasons. TW does not show consistent results over the country, thus being found as the least reliable alternative. Although S–W is highly correlated with the FAO P–M model, the application of the model in a data-scarce region is more constrained, as it requires more parameters than the FAO P–M model. Thus, the study suggests employing alternative methods based on the region of the country instead of one single method across the entire country.


Introduction
Evapotranspiration (ET) is the loss of surface and soil water to the atmosphere as water vapour by the combined actions of the two processes: evaporation from surface water bodies, bare soil and other surfaces that intercept rainwater, and transpiration from plants [1]. ET plays a significant role in the global water budget; hence, an accurate estimation of ET is essential for water resources management, ecosystem management and climatological studies. ET is primarily driven by soil condition, plant type, plant development stage and weather parameters such as solar radiation, wind speed, humidity and temperature [1]. Because of the large number of influencing factors and their heterogeneity over a watershed, accurate estimation of ET is challenging. Thus, estimation of potential evapotranspiration (PET), the potential amount of water that could evaporate and transpire from a landscape vegetated with unlimited water supply to the surface [2], or reference ET (RET), evapotranspiration of a hypothetical reference crop (defined by Allen et al. [3] as a crop having height of 0.12 m, surface resistance of 70 s.m −1 , and an albedo of 0.23) having unlimited water supply to the surface, have been widely used instead of actual Their study has shown that the better methods to substitute FAO P-M vary from one geomorphic unit to another; however, in general, radiation-based methods performed better than temperature-based methods. In contrast to the above studies, Lu et al. [4] conducted a comparative study in the southeastern United States evaluating six PET methods, including TW, with respect to actual evapotranspiration estimated by the water balance method instead of the FAO P-M method. Out of all, TW has been found by them as the worst performing model. Li et al. [21] compared remote sensing-based PET models as an alternative to data-rich models. Srivastava et al. [7] evaluated a satellite-based remote sensing method for a river basin which is rich with paddy lands, in eastern India.
Sri Lanka, an island nation in the Indian ocean, although blessed with heavy tropical rainfall, which unfortunately has highly uneven spatial distribution, was divided into three major climatic zones, i.e., the wet zone, intermediate zone and dry zone, based on the spatial distribution of rainfall. The dry zone of Sri Lanka can be characterized as a semi-arid region, where as per the literature, ET may have accounted for a loss in total rainfall ranging from 30% to more than 60% [22]. Hence, ET can be considered to be a major source of water loss in Sri Lanka's water budget, as the dry zone is composed of two-thirds of Sri Lanka (by area as in Figure 1). In addition, the dry zone is extensively irrigated. Extensive irrigation increases ET significantly, and therefore, the local water budget will be changed significantly [23]. The island nation's probable future vulnerability to climate changes and enhanced future agricultural needs to feed the rapidly growing population may result in far greater ET, by placing an enormous burden on the country's water budget. Thus, accurate quantification of ET over Sri Lanka is paramount for efficient and sustainable management of the country's water resources. The most accurate local ET can be obtained by using lysimeters, which are very costly and are not available in Sri Lanka. This necessitates identifying alternative PET methods which can estimate ET over Sri Lanka with greater certainty. No matter how important the ET estimation in Sri Lanka's water resources and irrigation management is, the attempts to identify better PET models are still scant in Sri Lanka.
Therefore, this study was formulated to compare the performances of two PET models, namely Shuttleworth and Wallace (S-W) and TW, across Sri Lanka with respect to the FAO P-M method. Following many previous studies, the FAO P-M model was considered to be the reference model for this comparative study. The pan method is one of the widely used methods when necessary meteorological data to estimate PET based on other equations are lacking [10]. Despite being the simplest low-cost technique to estimate PET [10], the pan method only provides evaporation from an open-water surface without accounting for transpiration [9]. However, as a general practice, the pan method has been commonly applied to account for evapotranspiration in hydrological analyses in Sri Lanka, as well as in some other countries [24][25][26]. By considering this fact, this study also evaluated the performance of the pan method, mainly to analyse the correlation of the pan method with estimated evapotranspiration values and then to investigate whether there is a necessity for computing coefficients that can bridge pan evaporation and evapotranspiration.
The S-W method, which takes both evaporation and transpiration into account by assuming two source-based (the crop and the substrate soil) schemes and balancing the energy exchanged between soil; canopy; air space between soil and canopy; and the atmosphere above the canopy, was also contrasted against the FAO P-M model. This attempt is unique to this study, as the data-rich S-W model has not been considered for previous comparative studies. Selection of the S-W model for the present study was followed by the study of Senathilake et al. [27], in which the S-W model was evaluated against pan method values and was found to be applicable to Sri Lanka. Since both FAO P-M and S-W models are complex with many meteorological parameters that are not easily obtainable, especially in data-scarce regions, and the pan method can be a tedious and time-consuming task [26], it was decided to examine the performances of an empirical equation, which has a minimum data requirement, with the standard method. Thus, the temperature-based TW method, which was developed by Thornthwaite [12], was evaluated in the present study. TW, which was developed for a humid climate, requires only the air temperature and duration of the site's daytime to estimate PET. With a lack of parameters that influence PET, TW estimates may be inaccurate in some regions, especially those with dissimilar climates to where TW was first established. However, if the TW model is proven to produce reliable estimates, it would be a better method for data-scarce regions, as temperature might be the only parameter being measured continuously in such regions. easily obtainable, especially in data-scarce regions, and the pan method can be a tediou and time-consuming task [26], it was decided to examine the performances of an empirica equation, which has a minimum data requirement, with the standard method. Thus, th temperature-based TW method, which was developed by Thornthwaite [12], wa evaluated in the present study. TW, which was developed for a humid climate, require only the air temperature and duration of the site's daytime to estimate PET. With a lack of parameters that influence PET, TW estimates may be inaccurate in some regions especially those with dissimilar climates to where TW was first established. However, i the TW model is proven to produce reliable estimates, it would be a better method fo data-scarce regions, as temperature might be the only parameter being measured continuously in such regions. Figure 1. Map of Sri Lanka indicating three major climatic zones, topographic zonation and th selected weather stations. Rainfall isohyets delineating climatic zones were drawn based on Climat Change Secretariat Sri Lanka [28].
Although several comparative studies have been conducted worldwide [4-7,10,19 20], there appear to be only a few previous studies comparing the performances o Although several comparative studies have been conducted worldwide [4][5][6][7]10,19,20], there appear to be only a few previous studies comparing the performances of different PET models that were carried out in the tropics, especially across South Asia. Therefore, not only could this study be employed by the water managers in Sri Lanka to identify the optimal method for the estimation of PET across Sri Lanka, but it would also contribute to addressing the aforementioned research gap because of its relevance to (1) many tropical Hydrology 2022, 9, 206 5 of 26 counties, as Sri Lanka is a tropical country featured by high temperatures, high humidity and unevenly distributed monsoon rainfall (both temporally and spatially); and (2) most of south Asia, as Sri Lanka, finds many similarities to the rest of South Asian countries in terms of tropical climate; environmental aspects; irrigation practices; data scarcity and socio-economic conditions.

Study Area
Sri Lanka (please refer to Figure 1) is a tropical island in the Indian ocean, located between 5 • and 10 • N latitude and 79 • and 82 • E longitude, with a total geographical area of 65,610 km 2 , comprising 62,705 km 2 area of land and 2905 km 2 area of water bodies. These water bodies comprise 103 distinct natural river basins and an extensive network of tanks and reservoirs (about 13,000). Approximately two-thirds of the country's landmass is lowlands with elevations less than 100 m. For highlands, the elevations vary from 100 to 2500 m approximately, with the highest mountain peak of 2525 m located in the country's central part [29].
The rainfall distribution is governed by the two major monsoon seasons: southwest monsoon (SWM) from March to September and northeast monsoon (NEM) from December to February. The wet zone is separated by the 2000 m annual average rainfall isohyet and is on the western and southwestern slopes. In the wet zone, rainfall ranges from 2000 to over 5000 mm, with an annual average rainfall of about 2400 mm. On the other hand, in the dry zone, the annual average rainfall is about 1450 mm, and rainfall can even be below 1000 mm. In addition to the two major monsoons and the inter-monsoon rains, tropical depressions that originate in the Bay of Bengal frequently enter Sri Lanka, resulting in extreme rainfall events, which sometimes may exceed 500 mm d −1 . Nearly 35-45% of annual rainfall contributes to annual surface runoff. However, in most dry zone river basins, the runoff percentage is less than 35%, with the rest of the rainfall being lost as evaporation and groundwater recharge [29].
Mean annual temperature in the lowlands and highlands varies between 26.5 and 28.5 and 14.7-17.1 • C, respectively. In addition, the pan evaporation values show considerable temporal and spatial variations-varying between 795 and 1900 mm/year, with higher values recorded in the hotter dry zone. The climate of the country is characterised by high relative humidity, generally ranging between 75 and 95% [29].
The Department of Meteorology of Sri Lanka maintains only 22 main meteorological stations over the country. The measured meteorological data are expensive in Sri Lanka, and there is a scarcity of measured data. One of the major limitations was obtaining solar radiation data, as solar radiation is not measured at all 22 main meteorological stations. Moreover, some of the stations that were not selected had incomplete data due to various reasons, including the 30-year civil war (which ended in 2009), instrumental errors, and logistic errors. Therefore, only 14 stations were selected (refer to Figure 1) so that they mostly cover all climatic (wet, intermediate and dry) and topographic (hilly areas and lowlands) zones of the country. The FAO P-M equation for calculating RET can be expressed as follows (Equation (1)) [3].
where ET o is the total daily RET (mm), ∆ is the slope of saturation vapour pressure curve (kPa • C −1 ), R n is the net incoming radiation (kPa), G is the soil heat flux (MJ m −2 ), γ is the psychrometric constant (kPa • C −1 ), T is the average daily temperature ( • C), U 2 is the wind speed at 2 m height (ms −1 ), e s and e a are the saturation vapour pressure (kPa) and the ambient vapour pressure (kPa), respectively. The equations required to estimate various parameters are presented in the Appendix B. However, similarly in a few recent studies (e.g., [5,6].), this study considers RET estimated by the FAO P-M method to be the same as the PET.

The Shuttleworth-Wallace Model
In the S-W model, total PET is computed as the summation of two major evapotranspiration components: soil evapotranspiration and transpiration from the dry canopy. The S-W model for estimating PET is presented in Equations (2)-(9) [17]: R c = (∆ + γ)r ac + γr c (8) where λ is the latent heat of vaporisation (MJ kg −1 ), ET o is the total daily PET (mm), ρ is the air density (kg m −3 ), C p is the specific air heat at constant pressure (= 1.013 × 10 −3 MJ kg −1 • C −1 ), D is the water vapour deficit at the reference height (kPa), r as is the aerodynamic resistance between the soil surface and canopy air space (s m −1 ), R ns is the net radiation at the substrate surface (MJ m −2 ), r a is the aerodynamic resistance between canopy source and reference level (s m −1 ), r s is the soil resistance (s m −1 ), and r ac is the bulk boundary-layer resistance of the canopy (s m −1 ). The necessary relationships to estimate the parameters of the S-W model are given in the Appendix B.

The Thornthwaite Model
The Thornthwaite [12] model is based on the length of daytime and monthly mean air temperature and is presented in Equation (10).
where ET o is the total monthly PET (cm), L d is the daytime length, T m is the monthly mean air temperature ( • C), I is the annual heat index computed using Equation (11), whereas a is computed using Equation (12).
where T j is the mean temperature ( • C) for the month j (j = 1, 2, . . . ,12). a = 6.75 × 10 −7 I 3 − 7.71 × 10 −5 I 2 + 0.01791I + 0.49239 (12) Hydrology 2022, 9, 206 7 of 26 Daytime length is the time from sunrise to sunset in multiples of 12 h. Being a tropical country throughout the year, the sunrise to sunset duration can be taken as an average of 12 h, which gives L d = 1 for all months in this study.

Data Acquisition
The study requires both meteorological and spatial data ( Table 1). The meteorological data, i.e., time series data of maximum and minimum daily temperatures, relative humidity (day and night), solar net radiation, wind speed, and pan evaporation measured at the selected meteorological stations from 2009 to 2019 were purchased from the Department of Meteorology of Sri Lanka.  (Table A1) Zhou et al. [26] Spatial data were obtained from various sources as described below. Elevation data are required for the S-W model to calculate the atmospheric pressure (please see Equation (A4) in the Appendix B), which is then used to compute ρ (Equation (A3) in the Appendix B). Elevation values were extracted from a digital elevation model (DEM) with a spatial resolution of 6 × 6 km, which was obtained from the Department of Survey, Sri Lanka. The land-use and soil maps of Sri Lanka, which were originally in vector format, were obtained from the Department of Survey, Sri Lanka, and then were transformed to raster data layers with each having a resolution of 6 × 6 km. In addition, the NDVI values (with a resolution of 6 × 6 km) were computed using remotely sensed Landsat images from 2009 to 2019, acquired from the United States Geological Survey Earth Explorer website. Furthermore, the monthly root zone soil moisture values were retrieved from the NASA Earth Data website (https://earthdata.nasa.gov/, accessed on 16 October 2020 from 2009 to 2019, and similar for the other data, these were converted to a grid of 6 × 6 km. The land cover threshold parameters required by the S-W model were gathered from Zhou et al. [30] and are given in the Appendix B (please refer to Table A1).

Model Building
Firstly, as required, the values of spatial data layers were extracted at the 14 selected stations. Then, all three models were developed in the Microsoft Excel office package, so that the models estimated monthly (for all 10 years) point PET at 14 selected stations.
Secondly, an attempt was made to develop spatial distributions of PET over the country, as estimated by each method. For this purpose, the entire country was represented by a 6 × 6 km grid. Since the selected 14 meteorological stations are located one dozen kilometres apart, having a grid with a finer resolution would not necessarily have improved the model's performance. Furthermore, a coarser resolution than the selected could not cover the entire country, with a substantial amount of square kilometres missing in the coastal areas. However, the selected grid can still be seen as a finer grid compared to the spatial distribution of meteorological stations, which must be noted as one of the major limitations of this attempt. Owing to this limitation, more reliable results can be expected near the stations, and further away from a station, any model may have to be applied with lesser confidence.
Since the meteorological data were point data, as the first step in spatial modelling, thematic layers (having grid size of 6 × 6 km) of each required meteorological data were prepared in ArcGIS 10.5. The thematic layers were prepared based on the average monthly data (averaged from 2009-2019) at 14 selected stations, as there are 12 (for each month) thematic layers of each meteorologic parameter. The thematic layers of temperature, relative humidity and evaporation were generated by inverse distance-weighted (IDW) interpolation [27]. The Universal Kriging interpolation method was used to obtain the spatial distribution of solar net radiation and wind speed [27]. Thematic layers of meteorological data for June (a mid-year month was selected for visualization) are shown in Figure 2.
Each and every grid cell value of all the required parameters were extracted into spreadsheets to create a spatial database, and all three models were again developed considering the spatial base. This process generated monthly PET (only for an average year) for each grid cell. Then, the outputs of all three methods were again exported to ArcGIS 10.5 in order to generate final PET distribution maps. Thus, the generated PET distribution layers also contain grids of 6 × 6 km spatial resolution covering the entire country. Then, these estimations and pan evaporations were compared with each other.

Performance Evaluation
To evaluate the performances of other PET methods in comparison to FAO P-M, three statistical indicators, namely, root mean squared error (RMSE), bias, and Pearson correlation coefficient (R), were calculated. Only the mean monthly point PET values at the 14 stations were considered under statistical performance evaluation. The governing equations for RMSE and bias are given in Equations (13) and (14), respectively.
where n is the number of PET estimates, PET o is the PET estimated by the other three methods (mm), PET PM is the PET estimated by FAO P-M. Pearson correlation coefficients were computed in R statistical programming language. A sensitivity analysis was carried out to assess the influence of each parameter on FAO P-M model estimates by increasing each parameter by 5%, 10%, 15%, 20% and 25%, one at a time, while others were kept at original values and computing per cent error. Due to the complexity of the S-W model, no sensitivity analysis was performed on the model. All the models were developed in MS Excel, and all performance evaluations except for Pearson correlation were also carried out in MS Excel. Hydrology 2022, 9, x FOR PEER REVIEW 9 of 28

Performance Evaluation
To evaluate the performances of other PET methods in comparison to FAO P-M, three statistical indicators, namely, root mean squared error (RMSE), bias, and Pearson correlation coefficient (R), were calculated. Only the mean monthly point PET values at the 14 stations were considered under statistical performance evaluation. The governing equations for RMSE and bias are given in Equations (13) and (14), respectively.
where is the number of PET estimates, is the PET estimated by the other three methods (mm), is the PET estimated by FAO P-M. Pearson correlation coefficients were computed in R statistical programming language.
A sensitivity analysis was carried out to assess the influence of each parameter on FAO P-M model estimates by increasing each parameter by 5%, 10%, 15%, 20% and 25%, one at a time, while others were kept at original values and computing per cent error. Due to the complexity of the S-W model, no sensitivity analysis was performed on the model. All the models were developed in MS Excel, and all performance evaluations except for Pearson correlation were also carried out in MS Excel.

Results and Discussion
Due to the unavailability of long-term data, especially solar radiation, calibration of the FAO P-M model was not carried out. As per the literature, the FAO P-M method is considered to be the standard RET (PET as considered by this study) method, since the FAO P-M can be used globally without any need for extra adjustments to model parameters [4]. The same reason negates the need for calibration of the FAO P-M model for different study areas. Therefore, as a general practice, many previous comparative studies [5,6,10,18,19] compared the performances of various PET models with respect to the FAO P-M model, taking the FAO P-M model as the standard method without calibrating. Furthermore, the range of PET values estimated by the FAO P-M model, 4.0-6.5 mm d −1 , was within the range of RET (2.9 to 10.8 mm d −1 ) obtained in Tamil Nadu, India, by Subburayan [31]. Although Tamil Nadu is characterized by similar climatology as in Sri Lanka, (as per Subburayan [31]) the range of (max to min) temperature and humidity are respectively greater and lower than those of Sri Lanka. Furthermore, the wind speed, which is 7.7 ms −1 [31], is significantly larger than in Sri Lanka. Most importantly, average number of sunshine hours is 7.8 h [31], whereas in Sri Lanka, it is around 12 h. Hence, the net solar radiation in Tamil Nadu is probably lower than in Sri Lanka. All above mentioned variations in Tamil Nadu could be the reason for a wider range of RET estimated by Subburayan [31]. However, another study, Bapuji Rao et al. [32], which also utilized FAO P-M to estimate PET across India, have found that the annual PET ranges between 1700 and 1950 mm in southern India. The annual PET by FAO P-M across Sri Lanka varies between 1460 and 2380 mm, which is wider than in southern India as estimated by Bapuji Rao et al. [32]. However, as the boundaries of the Sri Lankan PET (as estimated by FAO P-M) range are associated with extremes such as in Nuwara Eliya (which is discussed later), the average PET range in Sri Lanka is well within this South Indian PET range. Furthermore, Owusu-Sekyere et al. [20] found that in Cape Coast, Ghana (located between 5 • and 6 • N latitude), RET ranges between 3.5 and 4.42 mm d −1 . This range is similar to the current study's FAO P-M model estimates near the Sri Lankan southwest coast, which is located in between the same latitudes as Cape Cost. As given by FAO [33], annual RET across Sri Lanka varies approximately from 1000 to 2000 mm. As this study estimated, PET at only three locations falls outside of the annual RET range of Sri Lanka, which is given by FAO as a range from 1000 to 2000 mm. Therefore, with minor deviations at the two ends of the range, PET estimated by the FAO P-M model in Sri Lanka agrees with previous studies conducted in similar study areas, confirming the correctness of the FAO P-M model results of this study.
The results of the sensitivity analysis are shown in Figure 3. Table 2 summarises the Pearson correlation coefficients calculated by two correlation tests, which were performed to compare how each meteorological parameter affects the FAO P-M and S-W estimates. The FAO P-M shows the least sensitivity to the wind speed, which shows the least correlation to PET estimated by the FAO P-M. This suggests that wind speed has a comparatively small impact on FAO P-M estimates. This observation is in agreement with previous studies, such as Ngongondo et al. [19] and Nandagiri and Kovoor [34]. Similar to the FAO P-M model, the least correlation with the S-W model estimates are shown by the wind speed.
The p values prove that despite the models' least sensitivities to wind speed, the correlation of the parameter with PET estimates is significant in both FAO P-M and S-W models. All the climatic parameters in the FAO P-M model were found to be significantly correlated with the model's estimates in all climatic zones in India as well [34].   The key role played by the net solar radiation in PET estimation by the FAO P-M model is reflected by the model's highest sensitivity ( Figure 3) and correlation (Table 2) to the net solar radiation. Thus, in other words, net solar radiation is the major driving factor in the FAO P-M model in the Sri Lankan context. Solar radiation is hardly measured at weather stations, even in developed countries such as the USA [35]; thus, obtaining directly measured solar radiation data in data-scarce regions is exceptionally challenging. A similar fact, i.e., the FAO P-M model is highly sensitive to the most challenging parameter to measure, has been mentioned by Owusu-Sekyere et al. [20], confirming the current study's findings. Although not as high as in Sri Lanka, FAO P-M has shown the highest sensitivity to solar radiation in the Yanhe River Basin, China [6]. This high sensitivity basically limits the application of the FAO P-M model with high confidence, especially in data-scarce regions.
Nevertheless, there are several empirical methods suggested by previous studies (e.g., [30,34]) to estimate solar radiation rather than relying on direct measurements; some of the other parameters required by these methods are also difficult to obtain. Contrary to the FAO P-M model results, the S-W model shows, although positive, a minimum correlation with the net solar radiation ( Table 2), indicating that net solar radiation is not the primary mechanism driving the S-W model estimates. The FAO P-M model shows second and third higher sensitivities to relative humidity and temperature (Figure 3).  The key role played by the net solar radiation in PET estimation by the FAO P-M model is reflected by the model's highest sensitivity ( Figure 3) and correlation (Table 2) to the net solar radiation. Thus, in other words, net solar radiation is the major driving factor in the FAO P-M model in the Sri Lankan context. Solar radiation is hardly measured at weather stations, even in developed countries such as the USA [35]; thus, obtaining directly measured solar radiation data in data-scarce regions is exceptionally challenging. A similar fact, i.e., the FAO P-M model is highly sensitive to the most challenging parameter to measure, has been mentioned by Owusu-Sekyere et al. [20], confirming the current study's findings. Although not as high as in Sri Lanka, FAO P-M has shown the highest sensitivity to solar radiation in the Yanhe River Basin, China [6]. This high sensitivity basically limits the application of the FAO P-M model with high confidence, especially in data-scarce regions.
Nevertheless, there are several empirical methods suggested by previous studies (e.g., [30,34]) to estimate solar radiation rather than relying on direct measurements; some of the other parameters required by these methods are also difficult to obtain. Contrary to the FAO P-M model results, the S-W model shows, although positive, a minimum correlation with the net solar radiation ( Table 2), indicating that net solar radiation is not the primary mechanism driving the S-W model estimates. The FAO P-M model shows second and third higher sensitivities to relative humidity and temperature (Figure 3). Contrary to the present study's findings, Nandagiri and Kovoor [34] discovered that the maximum temperature significantly impacts FAO P-M model output; however, solar radiation was not incorporated in the factor correlation study of Nandagiri and Kovoor [34].
The FAO P-M model estimates are positively correlated with net radiation and wind speed ( Table 2), which is in agreement with Luo et al. [6], indicating a proportionate relationship between said parameters and the model estimates. Generally, an inverse relationship between relative humidity and PET is expected, and has been found by previous researchers (e.g., [6]). The negative correlation coefficients demonstrated by this study confirmed that the said expected relationship is valid for both the FAO P-M and S-W models. The average temperature was found to have a negative correlation with the FAO P-M model estimates, which has also been observed in a few locations across India by Nandagiri and Kovoor [34].

Comparison of PET Estimates at the Selected Weather Stations
Averaged (from 2009-2019) monthly PET estimates (by all four methods) were summed up annually and seasonally (two monsoon seasons) at fourteen meteorological stations and are shown in Table 3. Spatial distributions of average annual PET estimates (by all four methods) are displayed in Figure 4. Table 4 tabulates all three performance indicators: RMSE, bias and R-computed based on monthly PET estimates. A cross-comparison of the four methods was performed by finding the correlation coefficient of the monthly PET estimated by each method, (1) over the entire country, (2) for the wet zone, (3) for the dry zone and (4) for the intermediate zone. Hence, as in Table 4, the results of cross-comparison reveal the degree to which each model agrees with each other in different regions. All these statistical analyses were carried out only at meteorological stations, where input data are more accurate, hence lessening the error resulting from interpolation, as explained in Section 3.2.1. Performance statistics (Table 4) and spatial distributions (Figure 4) seem to be contradicting each other at some regions, as the first ones were derived from monthly series of PET, whereas the latter was derived from annual average series. Hence, Table 4 should be referred to obtain more accurate results, while Figure 4 is used to mainly qualitatively observe the spatial distribution of average annual PET.
A relatively higher influence of solar radiation on FAO P-M model estimates is reflected by the estimates showing a spatial variation similar to net solar radiation across Sri Lanka. Despite a few exceptions, i.e., Jaffna, Nuwara Eliya, and Colombo, dry zone estimates are the highest and wet zone ones are the lowest, while intermediate zone estimates range from mid to high (Table 3). Nuwara Eliya, which is characterised by comparatively (with the rest of the country) higher relative humidity throughout the year (Figure 2), shows the maximum annual (2349.1 mm) PET predicted by the FAO P-M model (Table 3) controverting to the expected low PET. These unexpected results can be explained by observing the net solar radiation and temperature. Nuwara Eliya records the highest net solar radiation in Sri Lanka throughout the year (Figure 2). Thus, solar radiation is the most influential parameter in the model with positive correlation results in larger PET estimates in Nuwara Eliya. When altitudes of weather stations are considered, Nuwera Eliya station is the highest (Figure 2); hence, the lowest mean annual temperature is recorded at Nuwera Eliya station. As the temperature inversely affects the model estimates (Table 2), the lowest temperature also has contributed to the significant deviation of the FAO P-M at Nuwera Eliya from the expected-lower than the rest of the country (Figure 4). Furthermore, as wind speed's effect on ET increases as the temperature decreases [1], higher PET estimates of the FAO P-M model can be supported by comparatively high wind speed recorded in Nuwera Eliya (Figure 2). The viability of this reason was supported by Owusu-Sekyere [20]. This study identified relatively low wind speed to maybe be the reason for low PET range estimated for a hot and humid study area in Ghana versus PET estimated by Subburayan [31] for a similar climate in India with higher wind speed. Galle has the lowest annual PET value predicted by the FAO P-M model (15,020.1 mm, Table 3), which can be attributed to the lowest mean annual net radiation and highest mean annual relative humidity (Figure 2).   Compared to the FAO P-M model estimates, the seasonal and annual estimates of the S-W model show great differences ( Table 3). As explained before, no sensitivity analysis was carried out for the S-W model; furthermore, the correlation between model results and land use/vegetation parameters was not analysed. Thus, the answer to the question "which parameter/s the S-W model is sensitive to the most?" was deducted based on a former study and qualitative observations in Figure 4b. Senathilake et al. [27] observed that the annual PET values estimated by the S-W model follow a pattern similar to the annual crop growth pattern, both temporarily and spatially, especially in extensively irrigated areas, such as Polonnaruwa, and no apparent similarities with the temporal or spatial distributions of other parameters were identified. Figure 4 also provides proof of S-W model estimates' closer relationship with the spatial distribution of cultivation in the country, as the highest mean annual PET values are distributed across a belt covering extensively cultivated areas (e.g., Anuradhapura and Polonnaruwa). Based on Senathilake et al.'s [27] study and Figure 4b, it was concluded that the S-W model estimates must have a higher correlation with land use/vegetation parameters, which can be employed to explain the vast difference between the estimates of S-W and FAO P-M models. As per Table 2, the S-W model results are also highly correlated with the average temperature. However, by observing Table 3 and temperature records at the weather stations, a somewhat proportionate pattern between temperatures and the estimates can be detected. For example, the highest annual and NEM PET values are simulated in Hambantota and Puttalam (Table 3) (Table 4) in Puttalam, Jaffna, Kurunegala, Hambantota, Rathnapura, and Polonnaruwa are less than 0.2. Pearson correlations (R) in all the areas mentioned above, except in Hambantota, are greater than 0.9. These statistics imply that the S-W model better fits the FAO P-M model in the areas above. When the spatial distribution of those is considered, it can be concluded that PET yield by the S-W model is closer to FAO P-M in the northwest of the country. This conclusion is further confirmed by the correlation coefficients in Table 5, according to which the S-W model and the FAO P-M model are highly correlated in the dry zone, as the north-western part of the country falls into the dry zone. The S-W model estimates being closer to FAO P-M estimates in the northwest is also valid for SWM and NWM, except for one significant discrepancy (R = 0.13) in Jaffna during NWM, for which no proper explanation can be found. Except for Colombo (which in Table 3 displays the overall highest annual bias and the highest annual RMSE), the annual highest bias and highest annual RMSE values are in the central hilly area of the country, suggesting that the S -W model would be the least reliable in this region when the FAO P-M model is considered as the standard model. The R value of the S-W model (0.49) computed for the wet zone in the region-wise analysis (Table 5), which is the second least correlation (second only to the intermediate zone, where R = 0.48) between the S-W model and the standard method, further supported the above judgement as nearly two-thirds of the upcountry belongs in the wet zone, and nearly one-third of the wet zone comprises hilly areas. As discussed earlier, the S-W method is mainly affected by vegetation cover but not solar radiation as in the FAO P-M method. Therefore, as the wet zone is rich with vegetation said weaker correlation between S-W and FAO P-M methods could be expected.
The TW model estimates do not have consistency in calculated bias values (Table 4); thus, it was concluded that the TW method could not capture the spatial distribution of FAO P-M model estimates. When comparing Table 4 with Table 5, it is apparent that the TW model overestimates annual PET in the regions where the annual average temperature is greater than 25 • C, with two exceptions in Hambantota and Polonnaruwa, which have recorded average annual temperature values greater than 27 • C, but with positive biases. Annual PET values are underestimated in Nuwara Eliya, Badulla, Bandarawela and Katugastota and hence in the central hilly areas, where the average annual temperature is lower than 25 • C. Generally, in the lowlands of Sri Lanka, mean annual temperatures vary between 25 and 32 • C, and in hilly regions mean annual temperatures drop below 25 • C. Thus, a conclusion was drawn that the TW model overestimates annual PET in lowlands and underestimates PET values in high lands. This inconsistency may be affected by the altitude. Many previous comparative studies, such as Chen et al. [9], Trajkovic and Kolakovic, [36] and Trajkovic et al. [37], have also observed inconsistent spatial bias in the TW model with respect to FAO P-M, but most researchers have identified no specific reason. Calculated statistics suggest that the PET estimate by TW for SWM follows a spatial distribution similar to the model's annual PET, with a minor anomaly in Hambantota. The TW model's PET predictions for NWM are highly irregular and do not seem to follow similar spatial patterns as the annual or SWM PET values. During the NWM, there is very little correlation between the TW and FAO P-M. The reason for this observation may be because of the possible superior influence of the SWM on the annual patterns than that of NWM, which could be attributed to several facts, such as the SWM season is longer than the NWM, brings in the most of annual rainfall and is the major cultivation season. This observation on the other hand implies that the TW estimates do not follow the temporal distribution of FAO P-M model estimates, which is not limited to the study area, but which has also been detected in other regions by previous studies [9,36,37]. Further, it can be seen that PET is overestimated in the SWM in the lowlands and is overestimated in the hills. Valle Júnior [8] claimed that the TW model was highly overestimated during the wet season-to which the present study agrees in the lowlands and during the SWM, the main wet season in Sri Lanka. The reason for the higher overestimation in wet seasons may be because the TW model does not consider humidity [8]. Chen et al. [9] also mentioned that inconsistent results (both spatial and temporal) derived from TW may be from the TW method's inability to account for any other PET driving factors, such as wind speed and humidity, which also reports vast variations spatially and temporally. Among all three methods, the pan method displays the most significant deviations from the FAO P-M estimates (Table 4), giving the lowest PET values ( Table 3). The pan evaporation values have a consistent positive bias (Table 4) over all the locations and all seasons, indicating underestimations of PET. This finding does not agree with many literature studies (e.g., [38]), which claim that actual (or potential) evapotranspiration is less than pan evapotranspiration. However, similar to the current study's findings, Tabari et al. [10], Nandagiri and Kovoor [34] and Xing et al. [39] also observed that the pan method underestimated the FAO P-M model consistently over their respective study areas in Iran, India, and Canada. Eagleman [40], explains the possibility of both underand overestimations observed by different researchers in different regions. As upper plant leaves absorb more energy, they evaporate more water, and this amount should be approximately the same, as the evaporation from shaded leaves is limited by lesser energy. However, with the presence of advected heat, lower leaves are supplied with more energy; thus, higher evapotranspiration (than evaporation from water surface) is possible in extensively vegetated areas. As the majority of the land mass of Sri Lanka is covered with vegetation, and advection heat transfer is also possible, as mentioned by Eagleman [40], an underestimation can be expected in Sri Lanka. Generally, said underestimations range from 10 to 20%, whereas in this study, the underestimations varied from 11 to 40% with two extremes of around 60% and 55% in Nuwera Eliya and Badulla, respectively. However, these larger underestimations were unsurprising when compared with the findings of a comparative study conducted in several humid and arid areas in southern India by Nandagiri and Kovoor [34] who found that the pan method produced extremely low PET compared to FAO P-M. In three southern states of India (namely Kerala, Andhra Pradesh and Karnataka), mean daily PET values by FAO P-M and pan methods respectively varied between 4.7-5.1 and 2.8-3.9 md −1 . In the present study, those values varied between 4.0-6.5 and 2.4-4.2 md −1 . These overlapping ranges of PET values by the two methods in the two countries clearly indicate that underestimations by the pan method are also within the same range of percentages. Looking at the possible unacceptance of the higher percentages of underestimation by the pan method, Nandagiri and Kovoor [34] suggested reassessing the pan coefficients for India. Based on extremely high underestimation by the pan method, this study also suggests re-evaluation of pan coefficients in the Sri Lankan context.
When the pan method is considered, the highest annual average PET values (pan method, Table 3) are in Hambantota, Puttalam, Polonnaruwa, Jaffna and Colombo (in the decreasing order)-stations in the hottest and driest areas of the country-while the lowest PET values are observed in Nuwera Eliya, Badulla, Rathnapura, Galle and Bandarawela (in the increasing order)-stations representing colder and humid climates of Sri Lanka. Therefore, it can be seen that the pan method undoubtedly reflects the temperature and humidity distribution of the country. Chen et al. [9] observed the pan method to reflect on different humidity conditions even across China [9]. The highest bias values ranging from 63% to 45%, are observed in the stations located in the mountainous region (Table 4): Nuwera Eliya (62.52%), Badulla (54.62%), Bandarawela (43.46%), and Rathnapura (44.60%). At low elevations, pan evaporation was underestimated compared to the FAO P-M model by lesser percentages ranging from 9% to 40%. When the stations are ranked from high to low based on the calculated bias (pan method, Table 3 Table 4 further confirm that pan evaporation values and the FAO P-M model outputs have a fairly good correlation in low elevations and a poor correlation in higher elevations. This is reinforced by the cross-comparison of the four methods (Table 5), which produced an R value of 0.78 in the dry zone, demonstrating a high correlation between the pan and FAO P-M methods. To further elaborate on PET values and biases, we can confidently state that the pan method works satisfactorily in hotter and drier environments. Temperature and humidity have more distinct variations in high and lowlands compared to other parameters in the FAO P-M model. Thus, those two could be identified as the parameters that determine the distinctive greater and lower correlations in the high and lowlands, respectively. However, within the scope of this study, no attempt was made to identify the dominant parameter. Weerasinghe [41] conducted a study to find how pan method PET values correlate with climatic parameters that could possibly influence PET, and a weaker correlation between temperature and pan evaporations was found, as well as the strongest correlation between humidity and pan evaporations. As the study was conducted in one of the hottest and driest areas of the country, it can be assumed that humidity is the leading parameter.
Although the PET observed via the pan method is a representation of the combined actions of solar radiation, wind, temperature, and humidity on evaporation, no clear relationship can be identified between the other parameters and pan evaporations, mainly because those parameters do not display a systematic distinct spatial distribution over the country, as temperature and humidity do.
Disregarding the Pearson correlation coefficient in Galle, the pan method has the largest annual R values, which are greater than 0.7, all over the country. The high correlations coupled with systematic deviations and greater correlation coefficients make the pan method a better substitute for the FAO P-M model. However, it is recommended to reassess the pan coefficients to establish appropriate zonal and seasonal corrections to the pan method errors to enhance performance. It seems that the requirement to adjust pan coefficients at regular intervals, probably more frequently due to rapid climatic changes, has been identified by many former studies worldwide (e.g., [9,35]).

Cross Comparison of Four PET Methods
This study classifies the performances of each PET method with respect to the FAO P-M model in different regions based on the R values as very good if R ≥ 0.90, good if 0.70 ≤ R < 0.90, moderate if 0.50 ≤ R < 0.70, low 0.30 ≤ R < 0.50 and poor if R < 0.30 (as per classification given by Mukaka [42]). The classifications are shown in Table 5. In the wet zone, the highest correlation of 0.49 with the FAO P-M has been exhibited by the S-W method; however, the highest correlation itself is weak. TW has a minus correlation in the wet zone, whereas the pan method shows a minimal positive correlation. Therefore, none of the methods can be considered a suitable replacement for FAO P-M in the wet zone. Even if all possible pairs of PET models were considered without limiting the comparison to FAO P-M, none of the methods show any strong relationship with each other in the wet zone. In the dry zone, the pan method has the strongest correlation (R = 0.78) with FAO P-M, which was labelled as good performance with respect to FAO P-M. As the S-W method is also highly correlated with the standard model, with an R value of 0.76, S-W was also identified as a good estimator of PET (Table 6). Having a moderate correlation (R = 0.67), the TW method can be categorised as a moderate performer, which can be employed instead of the standard method. These better correlations were showcased by considering only six stations spread over the dry zone, which account for nearly two-thirds of the country. Although the better correlations in the dry zone compared to other regions seem questionable with only six stations, given the data scarcity, the study has to rely on the obtained results. In the intermediate zone, the highest correlation with FAO P-M is shown by the pan method. Although not a stronger correlation, the correlation is still stronger than any correlations in the wet zone. Thus, it can be deduced that the pan method performs moderately in the intermediate zone (Table 6). Although it is not clear why some deviations exist among different models, it is plausible since these PET models have been developed for various climatic conditions and not for the Sri Lankan context. Thus, it is evident that the calibration of PET models is required [8] but may be hindered by data scarcity [30]. When considering the entire country, the TW model has the weakest correlation (R = 0.22) with the standard model. As region-wise correlation analyses revealed, the TW model performs the worst with respect to FAO P-M, except in the dry zone. These analyses indicate that the TW model would not be a good substitute for the FAO P-M model in Sri Lanka. Previous studies conducted in different parts of the world [4,9,37] have come to a similar conclusion-that the TW model performed worst. The poorest performances of TW may be because TW uses a local thermal index established for the east-central USA, where the model was developed [36], and that the index may be applicable to the climatic conditions of Sri Lanka. Thus, the TW model should not apply as is but must be calibrated to fit the local climate. Considering the entire country, the highest correlation (0.61) with the FAO P-M model is shown by the S-W model, which is the second strongest correlation the S-W model has with other models (Table 5) and seconds only to its correlation with the pan method (R = 0.67). Thus, performances of S-W against both the pan and FAO P-M methods were identified as moderate under this study. Therefore, the S-W was identified as a possible alternative in the Sri Lankan context. However, with S-W being the most complex method, the reliability of results depends on several variables, which are not easily obtainable; hence again, it cannot be recommended as a good alternative. As per the classification adhered to in this study, the pan method shows a moderate correlation of 0.55 with the FAO P-M method over the entire country. However, Lu et al. [4]. Considered that an R value of 0.57 indicates high correlation. Thus, considering both the classification mentioned above and the results tabulated in Tables 3-5, it would be reasonable to recommend the pan method as a better estimator-but not the best method. Furthermore, its performance can be enhanced by adjusting the pan coefficients [9,34].
Since the general practice in Sri Lanka is to employ pan measurements in hydrological analyses, the applicability of the other methods instead of the pan method was also evaluated in this study (Table 5). Over the entire country, among all the possible pairs of PET models, the TW and pan methods have the highest correlation (R = 0.69). The reason for this correlation may be that the pan method tends to follow temperature patterns (as explained before), which is the only parameter in the TW method. Despite the higher correlation between the TW and pan methods, given the underestimations by the pan method with respect to the standard method, the study suggests employing the FAO P-M method as the general practice in Sri Lanka whenever possible. This suggestion is because the three methods exhibit great differences compared to the FAO P-M method, making it difficult to agree on the best substitute for the FAO P-M. If the radiation data are uncertain, following Chiew et al. [43], the study recommends establishing appropriate zonal and seasonal corrections to the errors in the pan method at the meteorological stations so that the computed pan coefficient could be applicable to the zone in which the station belongs to. However, as this study pointed out, selecting a model that performs well in the region of focus is always good practice.

Conclusions
FAO P-M estimates showed the highest correlation with net solar radiation with a higher Pearson correlation coefficient value (R = 0.88), and the model is highly sensitive to net solar radiation. The TW model generated inconsistent results with minimal correlations; thus, it was identified as the least reliable method to use in place of the FAO P-M model. Both S-W and pan methods have shown consistent positive bias, indicating underestimation of PET over the entire county in all seasons. In low elevations, especially in the dry zone, all three models performed well compared to the models' performances in the central hilly area. No model was found to be strongly correlated with the FAO P-M model or each other in the wet zone. The pan measures were found to have the strongest correlation with the FAO P-M in the dry zone. In different parts of the country, different methods seem to produce a better correlation with the FAO P-M model; thus, the study recommends selecting different alternative methods depending on the focused area of the country instead of recommending a single alternative for the FAO P-M. In the case of island-wide hydrological analysis requirements, the S-W model can be recommended, as the S-W model showed the highest correlation with the FAO P-M over the entire country. Nevertheless, the S-W model, being a complex one, which requires a large number of meteorological data, some of which are not readily available in data-scarce regions, limits the application of the S-W model with certainty. Further, the study suggests replacing the pan method, which is the widely used method to estimate PET in Sri Lanka, with the standard FAO P-M model. FPAR min = 0.001, FPAR max = 0.95, SR min is the SR estimated for NDVI at 5% vegetation population (NDVI at 5% = 0.039 globally), SR max is the SR estimated for NDVI at 95% vegetation population (NDVI at 95%; see Table A1). LAI max values for different vegetation types are given in Table A1.
Estimation of the water vapour deficit at the reference height Zhou et al. [ θ c is equal to 0.75 times the saturated soil moisture Estimation of r s -the soil resistance Shuttleworth and Wallace [17] r s = 0 sm −1 at saturation point   [30]. Tall vegetation, 1 to 3; rest short vegetation.