Annual Actual Evapotranspiration Estimation via GIS Models of Three Empirical Methods Employing Remotely Sensed Data for the Peloponnese, Greece, and Comparison with Annual MODIS ET and Pan Evaporation Measurements

: Actual evapotranspiration (ETa) has been insufﬁciently investigated in Greece. This study aimed to estimate annual ETa by empirical methods (Turc, modiﬁed Turc, and Coutagne) for the Peloponnese, Greece, a Mediterranean testbed, between 2016–2019, four of the warmest years since the preindustrial era, and compare them to MODIS ET. Furthermore, measurements of annual pan evaporation (Epan) were performed for two Class A pan stations in the Peloponnese with different reliefs and conditions. The empirical methods and statistical formulae (RMSD, MB, and NMB) were developed as models in ArcMap. The outcomes of the Turc method resembled MODIS ET ranges for all years, followed by those of Coutagne. The estimates by the modiﬁed Turc method were almost identical to MODIS ET. Therefore, the modiﬁed Turc method can be used as an alternative to MODIS ET (and vice versa) for the Peloponnese for 2016–2019. Moreover, the Epan at Patras University station (semiurban, low elevation) exhibited an upward trend resembling the trends of the empirical methods over the study years, whereas the Epan at Ladonas station (higher elevation, lakeside) required investigation on a monthly time scale. Additionally, the gradual decrease of pan-water icing at Ladonas in December (from 20 d in 2016 to 0 d in 2019) could imply an undergoing decrease in snowpack storage retention across the mountains of the Peloponnese.


Introduction
Evapotranspiration (ET), a parameter of major importance to the hydrological cycle, determines the availability of water and, subsequently, water management and irrigation [1][2][3][4]. Actual evapotranspiration (ETa) is a parameter that is affected by climate change and land use/land cover (LULC) change, and has been examined along with wildfires [5][6][7][8][9]. The difficulty in directly measuring ETa has led to the implementation of empirical methods, algorithms and models with different levels of complexity, incorporating data from meteorological stations and remotely sensed data [10][11][12][13]. It is known that air temperature determines ET processes, as it expresses the state of energy of the system. Based on the established relationship between land surface temperature (LST), or "skin temperature", and near-surface air temperature (Tair), LST can be used as a satisfactory proxy for Tair despite the different physical meaning and responses of the two parameters to atmospheric conditions. Moderate Resolution Imaging Spectroradiometer (MODIS) has a daily overpass frequency, which makes it suitable for ETa estimation, but it is limited by its coarse resolution (1 km). LST is a global remote sensing product, and as such, its accuracy relies on several limiting factors: cloudiness and occasionally heavy aerosol loadings, land cover and the stage of crop growth, solar radiation and seasonality, and topography and elevation [14][15][16][17][18]. According to Jin and Dickinson (2010) [19], LST and Tair differ by 3.5 to 5.5 • C at most for maximum values (LST > Tair). There is consensus in the literature that LST night products from Terra and Aqua exhibit agreement with Tmin values because the lack of solar radiation does not cause fluctuations between Tair and LST [14,16,20]. LST day and Tmax values exhibit less agreement. Precipitation is the limiting factor governing the range of ETa values. This is especially true for the southern Mediterranean region, where it is anticipated that rainfall will increase during summertime in the near future [21,22].
There is a wide range of methods for ETa estimation, from simple formulae to complex algorithms [23][24][25][26][27][28][29][30]. In the case of Greece, ETa has not been investigated adequately at the regional scale. Crop evapotranspiration (ETc) at the local scale (crop plots) and reference evapotranspiration (ETo) at the basin scale have mostly been addressed in productive regions, especially the Thessaly plain, due to concerns regarding irrigation [31][32][33][34][35][36][37][38][39][40]. Due to the interest in the local scale, such studies usually focus on one station with large time series data, or process data obtained from up to three meteorological stations. In studies where ET has been estimated for Greece, sparsely distributed stations or satellite-sensed data were employed for either a few convenient days or short intervals [41][42][43].
Mavromatis and Stathis (2011) [44] defined ETa based on potential evapotranspiration (ETp) using 17 stations for Greece (only three stations for the Peloponnese) and found that the ETa trend followed the trend of ETp rather than that of precipitation. Soulis et al. (2010) [7] estimated ETa from the ETo of a basin subjected to wildfires, analyzing the soil moisture and land cover characteristics. However, in those cases, the ETa estimation depended on the accuracy of the ETp or ETo methods, introducing more degrees of freedom to the methodology. Voudouris et al. (2013) and Gudulas et al. (2013) [45,46] estimated annual ETa for three basins in northern Greece with satisfactory accuracy, using the Coutagne, Turc, and modified Turc empirical methods with ground-based inputs. Demertzi et al. (2020) [47] compared the tested method to annual estimates by the Turc and Coutagne methods for Greece as a whole, highlighting the suitability of these two methods for nonirrigated areas. Thus, according to the latter, applying these methods to the Peloponnese is rational, since more than 90% of the area is not irrigated.
Remotely sensed data (MODIS and Landsat products) for ETa estimation of two days in 2007 were employed for Karla Lake (Thessaly) and combined with data from one meteorological station. For the same study area, ETa was computed via the METRIC algorithm as an intermediate step for the assessment of two irrigation strategies. For the Thessaly plain, ETa was estimated for 21 days (summer of 2001) via the Carlson-Buffum and Granger empirical methods, incorporating remotely sensed data and data from three meteorological stations. Moreover, MODIS NDVI products were derived to compute ETa over a coastal grassland in northern Greece for one hydrological year (2013-2014) via the eddy covariance and Priestley-Taylor methods [48][49][50][51].
Pan evaporation (Epan) is used to define ETo via coefficients regarding environmental conditions and can be a proxy for ETp, since, despite their differences, they both quantify the evaporative demand of the atmosphere [52][53][54]. The most common type of pan evaporimeter is the Class A pan (d = 120.1 cm), and specific guidelines to run and maintain it have been set by the Food and Agriculture Organization of the United Nations (FAO) [52]. In Greece, there have been few Class A pan studies, as pans are scarce and usually utilized to define the irrigation needs of specific crop plots. Kitsara et al. (2009Kitsara et al. ( , 2012 [55,56] investigated the trends of Epan for Greece as a whole for two 30-year periods, 1979-1999 and 1999-2004. They derived measurements from 16 and 13 stations, respectively, and found positive trends during the latter period. To our knowledge, there is a gap in the literature regarding ETa for the Peloponnese, which is representative of the Mediterranean area, occupying almost one-sixth of the Greek territory, as well as Epan studies for recent years in the Greek territory. Coexamination of the latter, with ETa at the annual scale, could help to draw some useful conclusions. Therefore, the purpose of this study was to estimate the ETa of the Peloponnese for 2016-2019 using established empirical methods that were satisfactorily applied previously in Greece (Turc, modified Turc, and Coutagne), with monthly LST MODIS Terra products and precipitation data from 62 meteorological stations under the National Observatory of Athens (NOA) in (or near) the region of interest. In addition, MODIS net ET products were acquired, and annual estimates computed via ArcMap 10.6 (https://www.esri.com (accessed on 15 January 2021)) for comparison with annual empirical estimates (see Appendix A, Table A1). Reliable statistical indices were used to assess the agreement of empirical methods with MODIS ET. Since ETa is a component of the hydrological cycle, in which there is great interdisciplinary interest, both methods and the formulae of the statistical indices were developed in the model builder with raster inputs (i.e., MODIS products and Thiessen polygons), aiming to create user-friendly tools applicable to any area of interest. In addition, spot comparisons of ETa and MODIS ET to Epan from two extra stations (Class A pans) were performed at the annual scale. The produced maps for annual ETa and MODIS derived ET for 2016-2019 can be utilized by researchers, water managers and policy makers.
According to the World Meteorological Organization (WMO) [57], the study years (2016-2019) were four out of the five warmest years since the preindustrial era, and the study area is a suitable southern Mediterranean testbed. It combines different relief and LULC over short distances, experiences precipitation scarcity in the context of climate change, and is anticipated to have considerable increases in minimum, maximum, and extreme Tair values in the near future [21,22]. Therefore, the results attract interdisciplinary interest in fields such as hydrology and water management, remote sensing, climate change, and sustainability. The study employed ETa, MODIS ET, and Epan and lays a foundation for refining the applied empirical methods, identifying coefficients for the two distinguished Class A pan sites, and investigating possible interconnections among ET types in the Peloponnese.
In this paper, the Materials and Methods section follows the Introduction, describing the three empirical equations selected with the ground-based (precipitation) and remotely sensed (LST) data employed, MODIS ET products, and pan evaporation measurements from two extra Class A pan stations. The third section presents the study area. The fourth section includes the equations of each method, and the corresponding GIS models that implement those equations are included in the Supplementary Materials. The results are presented in Section 5, starting with descriptive statistics of the areal ETa estimates, followed by the statistical measures that assess the agreement between the estimates of the three empirical methods and their spatial distributions. Then, the similarities of each empirical method to MODIS ET and MODIS ET maps are shown. In the sixth section, the results are discussed in three ways: a comparison of ETa values from the empirical methods, a coexamination of empirical estimates, and MODIS ET and pan evaporation. Lastly, the main conclusions are summarized.

Empirical Equations Employing Meteorological and MODIS Data Implemented via GIS Models
Three empirical methods for annual ETa estimation were used (Turc, modified Turc, and Coutagne). The corresponding models to implement the methods were developed in ModelBuilder. These methods produced satisfactory estimates for Greece according to the literature. One user-friendly model, along with the corresponding Python script for the Turc and Coutagne equations, was developed in ModelBuilder (ArcMap 10.6) (https://www.esri.com (accessed on 15 January 2021)). The more complex modified Turc method needed four models to implement. All models consist of a sequence of preprocessing steps (e.g., re-projection, masking) and intermediate equations computed via a raster calculator. In the case of Greece, the overpass time of Aqua is nearest to midday, when local Tair is maximum. However, overestimation of Aqua LST compared to Tair values was reported for Greece during the summer (up to 5 • C) [15]. In addition, Zhu et al. (2013) [64] reported considerably smaller RMSE for day and night Terra LST products than Aqua products, for an area with a similar latitude to the current study area. Therefore, Terra LST was preferred over Aqua to be used in place of Tmax and Tmin in this research. Terra LST day and night products for 2016-2019 (eight annual scale images) were derived from FORTH (http://rslab.gr/downloads_lst.html (accessed on 10 March 2021)). They were used as maximum and minimum annual air temperature values (Tmax, Tmin), based on Kitsara et al. (2018) [15], who proved that MODIS LST products can be satisfactorily used as Tair proxies for Greece.
In addition, the use of Thiessen polygons is a reliable method for integrating precipitation, which has been satisfactorily used for Greek areas (either on the basin scale or for the whole country) [43,65]. Thus, Thiessen polygons were applied in this study area to create raster distributions of precipitation. The latter were used as inputs for models producing annual ETa, or as conditionals (see Coutagne's equations). Daily precipitation as a weather parameter was available online from the meteosearch.gr database (https://meteosearchmeteo.gr/ (accessed on 16 January 2021)) for the 62 automatic stations under NOA. The measurements of these stations were digitized and used in ArcMap 10.6 to form Thiessen polygons over the study area. Only data from stations with no more than two missing days per month were employed.

MODIS ET Products
MODIS products of net ET (MOD16A2V6; 184 images) were derived from the EarthExplorer (https://earthexplorer.usgs.gov/ (accessed on 21 February 2021)) and EARTHDATA (https://urs.earthdata.nasa.gov/ (accessed on 21 February 2021)) platforms. The aim was to assess the agreement between MODIS ET estimates and ETa values by empirical equations utilizing meteorological (precipitation) and remotely sensed (LST) data. An updated algorithm based on the Penman-Monteith equation was used to produce MOD16A2V6 images. This algorithm incorporates remotely sensed vegetation inputs and reanalysis meteorological data (air temperature, radiation, humidity) from NASA's Global Modeling and Assimilation Office (GMAO, v.4.0.0) [66]. The 184 images make up an 8-day composite that includes, along with other layers, the ET 500 m. After being appropriately preprocessed and reprojected in ArcMap, annual ETa values for 2016-2019 were derived using cell statistics (sum) and ignoring no-data pixels.

Pan Evaporation Measurements
A Class A pan provided by the Institute for Environmental Research and Sustainable Development (IERSD) of NOA (https://www.iersd.noa.gr/en/ (accessed on 11 July 2021)) was placed and run according to FAO guidelines [52] at the meteorological station of the Laboratory of Atmospheric Physics on the campus of the University of Patras. This is a semiurban location at the northmost part of the Peloponnese, at an absolute height of 24.8 m. Hourly Epan measurements were taken over 4 years (2016-2019). An automatic potentiometer, which converted water depth values to voltage values after proper calibration and data logger units, was used. The voltage signs of the potentiometer were frequently compared against voltage signs of a portable voltmeter. The produced depths were also compared to the water depths measured by a ruler. A pipe for overflow and an inox bar with a spike edge were adjusted at the pan, at the same depth where the water was maintained (20.2 cm) [52]. The latter helped with visual control when the pan was manually refilled with water. The daily precipitation for the Patras University station, and the meteorological parameters for completing any missing daily Epan values, were provided by the Laboratory of Atmospheric Physics of the University of Patras (http://mymeasurements.eu/u/lapup/meteo.php?lang=el (accessed on 11 July 2021)).
A second Class A pan was located near the artificial lake of Ladonas (Ladonas station) in the central part of the Peloponnese, at an absolute height of 420 m. That pan has been kept running for decades by the technical staff of the hydroelectric plant. The daily ET measurements for 2016-2019 were acquired and digitized. The Epan measurements (from Patras University and Ladonas stations) are presented and plotted along with the MODIS ET values and ETa estimates by the empirical methods for the same (annual) time scale. The values of the pixels, including the spots where the pans are located, were used for coexamination of ETa, MODIS ET, and Epan. The values were derived from the raster images produced by the empirical equations and the (annually accumulated) MODIS ET products. There were no other evaporimeters available in Peloponnese.

Statistical Indices Computed via GIS Models
Descriptive statistics (max, min, mean, std value) for each raster were calculated in ArcMap 10.6 (via calculate statistics). In order to evaluate the agreement between values from different empirical methods, the statistical indices root mean square deviation (RMSD), mean bias (MB), and normalized mean bias (NMB) ( Table 1) were developed in ModelBuilder [67][68][69]. Those indices are considered the most reliable for the statistical analysis of hydrological parameters [70]. Additionally, unpaired Student's t-test (mean difference, standard error of difference, p-value) and Cohen's D coefficient (effect size) were computed for the empirical methods [43,[71][72][73].

Turc Method
Annual ET a was calculated via a model (see Turc method in Supplementary Materials) that used LST Terra day and night images and meteorological data (precipitation) as inputs to implement Turc's equation (Equations (1) and (2)) [74], which had been previously used for Greece [45,46]. The model incorporates submodels for the computation of Thiessen polygons and the L function (Equation (2)).

Coutagne's Method
ETa was derived via a model (see Coutagne method in Supplementary Materials) that implements the Coutagne method (Equations (3)-(6)) [75], based on conditional selection (via the Con function in ArcMap). The model employs LST Terra products and Thiessen polygons (precipitation data) [46].

Turc Modified Method
The difference between the Turc and modified Turc methods is the way Tmean is calculated (Equation (7)) [45,76]: where P i stands for monthly precipitation calculated via an iterative model based on Thiessen polygons, with ground-based precipitation data as inputs (see modified Turc methods a and b in Supplementary Materials), and T i stands for monthly mean air temperature derived via three submodels employing LST Terra day and night images as inputs (see modified Turc methods a, b, and c in Supplementary Materials). ETa is then derived via a final model (see modified Turc method d in Supplementary Materials), consisting of two submodels which calculate T mean,mod and the L function. Iteration is not an asset of ModelBuilder, as it demands several submodels to implement the sums of products. However, models in optical programming are user-friendly and do not require a background in programming to run, and thus have interdisciplinary interest.

Descriptive Statistics of Empirical Method Estimates and MODIS ET
As shown in Table 2, the Turc method produced higher estimates of annual ETa than the other methods, followed by Coutagne. The modified Turc method showed mean annual values closer to MODIS ET values.

Statistical Measures between Annual Estimates by Empirical Methods
The statistical measures computed for the annual ETa estimates suggest that the Coutagne and Turc methods produced closer annual ET estimations, with RMSD = 60.04 mm y

Spatial Distributions of Annual ETa by Empirical Methods
The minimum values of annual ETa were similar for 2016-2018 (Figure 4a-d). There were similarities in spatial distribution between the three methods for each year. The maximum values were seen over the central part of the Peloponnese (with consistently higher elevation) for the four years, because of the well-known relationship between elevation and precipitation. It is also rational that in a Mediterranean region, areas with higher precipitation will exhibit higher ETa values. The spatial similarities are attributed to the fact that all three methods used annual precipitation data that were interpolated via Thiessen polygons.

Spatial Distributions of Annual ETa by Empirical Methods
The minimum values of annual ETa were similar for 2016-2018 (Figure 4a-d). There were similarities in spatial distribution between the three methods for each year. The maximum values were seen over the central part of the Peloponnese (with consistently higher elevation) for the four years, because of the well-known relationship between elevation and precipitation. It is also rational that in a Mediterranean region, areas with higher precipitation will exhibit higher ETa values. The spatial similarities are attributed to the fact that all three methods used annual precipitation data that were interpolated via Thiessen polygons.

Statistical Measures between Annual Estimates by Empirical Methods and MODIS ET Products
Statistical indices between empirical methods for the annual ETa and MODIS ET estimates (the latter derived via the sum of MOD16A2V6 8-day composites, ignoring no-data pixels, in cell statistics) indicate that the estimation of modified Turc was closer to MODIS's, with a mean bias error (MBE = |MB|) between 3.10 and 33.35 mm d −1 and NMBE (|NMB|) between 0.01 and 0.06 ( Figure 5). This is probably because even though their minimum and maximum values deviated considerably, the majority of their cell-to-cell values were close. Aiming to further investigate this contradiction, unpaired Student's t-tests were performed between the estimates of each empirical method and MODIS for every year, based on the use of the test by Dalezios et al. (2002) and Papadavid and Hadjimitsis (2010) [43,73] to compare the annual ET outputs from different methods. The estimates by MODIS and modified Turc showed minimum mean differences for the four years, which were statistically significant at 95% (p < 0.05) ( Figure 6). The effect size of the differences was small (Cohen's D < 0.24), thus their differences were small. Effect size is a sensitive measure that is meaningful to estimate the order of magnitude of the difference in cases where the sample size (N) is large enough to lead to small p-values (in the present study, N > 20,000 pixels). The estimates by Coutagne

Spatial Distributions of Annual MODIS ET
Comparing the maps of annual MODIS ET (Figure 7) to those of the empirical model ETa (see Figure 4a-

Spatial Distributions of Annual MODIS ET
Comparing the maps of annual MODIS ET (Figure 7) to those of the empirical model ETa (see Figure 4a-d), it is obvious that the Turc method produced a wider range of values than the other empirical methods, closer to that of MODIS.

Annual Values
The datasets of daily measurements of the Class A Epan were acquired for 2016-2019 from both stations (Patras University and Ladonas) ( Table 3). The maximum Epan values

Annual Values
The datasets of daily measurements of the Class A Epan were acquired for 2016-2019 from both stations (Patras University and Ladonas) ( Table 3). The maximum Epan values appeared in 2017, followed by those in 2019, which were close to the corresponding values in 2016.

Graphs of Annual Values of Epan, Empirical ETa and MODIS ET
As shown in Figure 8, the annual Epan values of both stations were plotted with ETa by empirical methods and MODIS ET (in mm). Epan (raw) data, rather than data previously converted to ETo, were plotted to investigate the trend of the measured parameters over the studied years and identify any similarities to the ETa trend by empirical methods. The line of annual values of Epan for Patras University exhibits some differences compared to the corresponding lines of ETa from empirical methods and MODIS ET (Figure 8a

Graphs of Annual Values of Epan, Empirical ETa and MODIS ET
As shown in Figure 8, the annual Epan values of both stations were plotted with ETa by empirical methods and MODIS ET (in mm). Epan (raw) data, rather than data previously converted to ETo, were plotted to investigate the trend of the measured parameters over the studied years and identify any similarities to the ETa trend by empirical methods. The line of annual values of Epan for Patras University exhibits some differences compared to the corresponding lines of ETa from empirical methods and MODIS ET ( Figure  8a)

Annual ETa by Empirical Methods
For 2016 (see Figure 4a), Turc estimated an annual ETa of 1013.7 mm, and the estimates by modified Turc and Coutagne were close to each other (882 and 894 mm, respectively). The minimum values were very close for all methods (211.4-220.4 mm). The distributions of annual ETa show similar segmenting of the Peloponnese, since the three methods depend on annual precipitation. Thiessen polygons were the same for all methods and years, but for modified Turc, the final distribution was affected by the L function. The eastern part of the Peloponnese experiences very low precipitation (214-694 mm); however, the locally elevated mean annual Tair (around Nemea and Argos, ranging between 16 and 22.6 °C) did not alter the ETa distribution, which was mostly affected by the precipitation polygons. The Coutagne method relies more on the precipitation contribution, since ETa is defined based on the condition of precipitation (P < l/8, P > l/2, or in

Annual ETa by Empirical Methods
For 2016 (see Figure 4a), Turc estimated an annual ETa of 1013.7 mm, and the estimates by modified Turc and Coutagne were close to each other (882 and 894 mm, respectively). The minimum values were very close for all methods (211.4-220.4 mm). The distributions of annual ETa show similar segmenting of the Peloponnese, since the three methods depend on annual precipitation. Thiessen polygons were the same for all methods and years, but for modified Turc, the final distribution was affected by the L function. The eastern part of the Peloponnese experiences very low precipitation (214-694 mm); however, the locally elevated mean annual Tair (around Nemea and Argos, ranging between 16 and 22.6 • C) did not alter the ETa distribution, which was mostly affected by the precipitation polygons.
The Coutagne method relies more on the precipitation contribution, since ETa is defined based on the condition of precipitation (P < l/8, P > l/2, or in between) and the majority of the created Thiessen polygons comply with the second condition (l/8 ≤ P ≤ l/2). ETa for 2017 (see Figure 4b) had a narrower range compared to 2016. Although the minimum values were very close among the three methods (242.60-248.20 mm), the maximum values differ. Turc produced higher estimates, followed by Coutagne. The distributions of Turc and Coutagne seemed to be more affected by the precipitation levels. The central part (Stemnitsa, Andritsaina) had the highest ETa by Turc and Coutagne, whereas modified Turc assigned the highest values to the western and central parts. The eastern part (except Kavos Malea) and the southwestern Filiatra area had low values with the former methods, but medium to high values with modified Turc. It seems that modified Turc has a tendency to counterbalance the outcome when high precipitation is combined with lower Tair, and vice versa.
The annual ETa was higher in 2018 (see Figure 4c)

Estimates of Empirical Methods Compared to MODIS ET
MODIS ET is an estimate of ETa produced by a sophisticated algorithm. However, underestimation of ETa has occurred not only in our implementation of ET but also compared with sites of the AmeriFlux network in the literature (by 26% [77] and 25% [78]), which were used to calibrate the updated algorithm [66]. As a rule, global satellite products contain more noise than ground-based measurements [79]. For annual estimates, modified Turc values had the minimum mean difference from MODIS ET values which, based on Student's t-test, was statistically significant at 95% and showed a small effect size  Figure 5). It is known that Turc performs better for areas with a precipitation index greater than 700 mm y −1 [40]. Moreover, the need to modify the way Tmean was calculated, which resulted in the modified Turc method, was verified. Overall, Turc modified is recommended for annual areal ETa estimation in the Peloponnese for 2016-2019 as a preferable alternative to MODIS ET. In the recent literature, MODIS ET was integrated into empirical equations of water recharge [80,81].

Pan Evaporation
Patras University station is located in a semiurban, low-elevation area where wind speed values are often high. The annual trend of Epan showed some differences compared to the trend of ETa by empirical methods and MODIS ET. The main difference was that the Epan value was higher in 2017 than 2018, whereas the opposite was found for ETa and MODIS ET (see Figure 8a). This is completely rational, since actual ET is subject to water stress limitations in contrast with Epan. While the empirical method estimates formed a similar line, the differentiation between those estimates and MODIS ET estimates was that the value for 2019 was lower than that for 2018, probably attributed to uncertainty in the global products.
The Ladonas Class A pan is located on the bank of Ladonas artificial lake (420 m elevation) surrounded by high mountains with sparse sclerophyll vegetation. On only 10 days during December 2016 the water of the pan was in the liquid phase, and this was the reason why the corresponding Epan value was very low. The line for annual Epan values was almost the opposite to the line for MODIS ET values (except for 2016) (see Figure 8b). The empirical method estimates generally formed a similar line for Ladonas annually, with the modified Turc values closer to the MODIS values. A distinguished peak for the Epan value in 2017 was noted and should be compared to monthly values, aiming to identify whether the evaporative demand of the atmosphere exhibited seasonal peaks or was consistently elevated during 2017. This peak potentially occurred because 2016 was the warmest reported year since the preindustrial era [57].
A second observation from Ladonas station was the substantial change in ice cover in the pan for the four successive winter seasons, reported as "total number of days with ice cover" and as shifting of ice cover days from December to January. Specifically, the water in the pan was in the solid phase (ice) for 12 days in January 2016 and 21 days in December 2016 (33 days in total). The total number of days gradually decreased to 14 days in total in 2019 (i.e., 14 days in January 2019 and zero days in December 2019).
The sequence demonstrates higher minimum Tair values over the studied years in line with the WHO's (2020) report stating that 2016-2019 were four of the five warmest years since the preindustrial era. Furthermore, it suggests a decrease in snowpack storage retention from year to year across the central mountains of the Peloponnese. In that case, increased runoff is inevitable, with implications for aquifer recharging and potential flooding events in the lowland, with socioeconomic impacts [22]. These indications were more apparent considering the Epan values [82].
The main limitation of the study is related to missing Ladonas Epan values in cases where solidification of pan water into ice occurred (January and December). In addition, the lack of an overflow pipe on the pan (at Patras University station) led to a few missing daily values during the rainy months (November-January). Those days were ignored, since Epan in such conditions (days with heavy rainfall or ice cover, but average wind speed) cannot be considerable. In addition, due to constraints in iterative programming in the ModelBuilder environment, the modified Turc method needs four models to be implemented, but the inconvenience is not important as those models are quick and user-friendly.

Conclusions
The annual ETa for 2019 showed the largest discrepancies among all compared pairs of empirical methods. That was the second warmest year since the preindustrial era and the third sequential year between 2016 and 2018. The empirical methods produced the highest maximum estimates for 2016, the warmest recorded year, and the lowest maximum values for 2017, implying that the evaporative demand could not be satisfied by the water availability.
In a comparison of estimates at the annual scale, it was found that all empirical methods produced narrower ranges of values than MODIS, and the range of the Turc method was the closest to MODIS. However, statistical indices suggest that modified Turc could serve as an alternative to MODIS annual ET for the Peloponnese for 2016-2019.
The line of annual Epan values for Patras University station showed only a few discrepancies when compared to the lines of ETa values of the empirical methods. This is a semiurban station located at low altitude at the northmost edge of the Peloponnese that often records elevated wind speeds. For Ladonas station, the lines of annual ETa values for the empirical methods appeared to agree with each other, but showed substantial discrepancies with the line of Epan values.
The peak annual Epan value seen for 2017, which was the lowest annual ETa value for the empirical methods for the Peloponnese as a whole, and is probably due to the fact that 2016 was the warmest year on record. Therefore, the increased evaporative demand was apparent in the Epan measurements, which involves no water limitations. The annual value for 2019 was at the same level as 2017, indicating that it was the second warmest year on record and the last in a sequence of very warm years. Another observation from Ladonas station was the substantial change in ice cover in the pan between years, which may imply a decrease in snowpack storage retention across the mountains of the Peloponnese (and thus an increase in runoff, with implications in terms of aquifer charging and flooding).
Computing multiplying factors is proposed to convert Epan into ETo for both stations (which differ in terms of latitude, elevation, land cover, and microclimate) by several methods in the literature, and comparing the estimated results between the empirical methods and MODIS ET for finer time scales (daily, monthly).