Reference Evapotranspiration (ETo) Methods Implemented as ArcMap Models with Remote-Sensed and Ground-Based Inputs, Examined along with MODIS ET, for Peloponnese, Greece

: The present study develops ArcMap models to implement the following three methods: Hansen, used a reference. three models implementing statistical indices (RMSD, MB, NMB) are also created. The purpose is threefold, as follows: to investigate the variability in the daily mean reference evapotranspiration (ETo) for the Decembers and Augusts during 2016–2019, over Peloponnese, Greece. Furthermore, to investigate the agreement between the methods’ ETo estimates, and examine the former along with MODIS ET (daily) averaged products. The study area is a complex Mediterranean area. Meteorological data from sixty-two stations under the National Observatory of Athens (NOA), and MODIS Terra LST products, have been employed. FAO PM is found sensitive to wind speed and depicts interactions among climate parameters (T, evaporative demand and water availability) in the frame of climate change. The years 2016–2019 are four of the warmest since the preindustrial era. Hargreaves–Samani’s estimations for the Decembers of 2016–2019 were almost identical to MODIS ET, despite their different physical meaning. However, for the Augusts there are considerable discrepancies between the methods’ and MODIS’s estimates, attributed to the higher evaporative demand in the summertime. The GIS models are accurate, reliable, time-saving, and adjustable to any study area.


Introduction
Evapotranspiration (ET) is a substantial parameter of the hydrological cycle, which determines, along with precipitation, water availability with implementations in irrigation and water management [1,2]. The reference ET (ETo) is a climate parameter that has been investigated as an indicator of climate change [3][4][5][6]. The difficulty of measuring ET led to the development of several ways of estimation, from simple empirical or physically based models to complex algorithms employing neuro-fuzzy or machine learning techniques [7][8][9][10][11]. The methods utilize data from meteorological stations as inputs, due to the scarcity or unavailability of the latter to remote sensed data. Because of the relationship between the reference evapotranspiration (ETo) and air temperature (T) (i.e., the latter as a proxy of the energy status of the system), T is for several formulae the only pre-required input. In addition, land surface temperature (LST) is a proxy for (near surface) air temperature. Therefore, LST is the most frequently satellite-derived product for ET estimation. MODIS daily overpassing frequency makes it the ultimate source of such data for ET estimation, with the shortcoming of low resolution (1 km). LST can be used as a proxy of air temperature since there is a strong relationship between LST or "skin temperature" and near-surface air temperature (T), although those two parameters have different physical most convenient). At last, although MODIS is a suitable remote sensed data source for ET studies due to its frequency of overpassing, there are few studies incorporating MODIS products for Greek areas. This is because those studies are not focused on ETo investigation as a primary objective [41][42][43][44]. Kitsara et al. (2018) [17] used MODIS LST in place of T in Hargreaves-Samani (HS), with Penman-Monteith as a reference for Greece. They found similar errors in ETo from the aforementioned methods to more complex algorithms in the literature. Paparrizos et al. (2017) [45] investigated several methods, among which was Hargreaves method to estimate the potential ET (ETp) (three river basins of Greece). They concluded that the temperature-based equations led to more accurate results, followed by the radiation-based equations. Mamassis and Panagoulia (2014) [46] verified that T (as a proxy of Rs) has a governing role on annual evaporation, while the relative humidity, wind speed and sunshine duration follow. During the summer months, T affects the most monthly evaporation. The influence of T and sunshine duration on evaporation is greater during the summer months, and the opposite. Efthimiou et al. (2013) [47] computed ETo with HS for NW Greece. Hargreaves-Samani (HS) exhibited generally good performance (low RMSE compared to FAO PM) in the winter season (October-February), but poor performance in the summer season (March-September). Diamantopoulou et al. (2011) [48] (North and central Greece in 2005) found RMSE between HS and FAO PM from 0.55 (northmost station) to 0.72 mm (southmost station), with HS underestimating ETo.
To our knowledge there are no reports of studies on ET specifying over Peloponnese peninsula, although it is an area that occupies almost 1/6th of Greece territory, with terrain combining wide plains with high mountains and populated cities. The sharp relief and the high spatial variability of the meteorological conditions over very short distances provide a suitable testbed for the comparison of the selected methods. The purpose of the study is threefold, as follows: to investigate the variability of the daily mean ETo for the Decembers and Augusts during 2016-2019, over Peloponnese, Greece. Furthermore, to examine the agreement between the methods, as well as the differences between the daily mean ETo and MODIS ET (daily) averaged products. The methodological novelty of the study is the way the methods are implemented and the statistical indices are computed, which is via GIS models with interactive character. Those models are developed as python scripts, which can run as they are for any study area by providing the corresponding inputs. The statistical models are user-friendly and applicable to any study area as well. They compute areal indices by processing large images with more than 20,000 valid values (e.g., MODIS products, ET distributions produced by the above-mentioned models, interpolated meteorological data), which provides accuracy and reliability. The study exhibits interdisciplinary interest regarding remote sensing, hydrology, climate change, water management, and sustainability. The years 2016 and 2019 are the warmest years on record since the pre-industrial era [49]. The Mediterranean regions are anticipated to be subjected foremost to climate change consequences in the near future, due to the undergoing increase in T, the decrease in rainfall in the summertime, and the shifting of precipitation patterns [50][51][52].
The study is organized into four sections after the introduction. In the Materials and Methods, the used equations and statistical measures are presented. Moreover, the followed steps of the methodology and data processing are elaborated. Section 3 provides the characteristics of the study area, which makes it a suitable testbed for the research. In Section 4, the statistical indices that have been implemented as GIS models are described. The results are presented in Section 5 with a pattern. In Section 5.1, the descriptive statistics of the areal daily ETo and MODIS ET estimates are presented. Section 5.2 includes the spatial distributions of ETo produced by the three methods for the 8 months of interest and the graphs of the statistical indices, which assess the agreement between the methods' estimates (in pairs). In Section 5.3 the spatial distributions of MODIS ET are displayed, followed by the graphs of statistical indices, which assess the differences between each method's ETo estimates against the MODIS ET estimates. Those indices are displayed as graphs of each statistical measure against the month of interest (i.e., Augusts or Decembers) of the period 2016-2019. The graphs facilitate the visual representation of the variability of estimations for the months of interest through the examined years. In Section 7, the results are discussed via the following two axes: (i) the ETo estimates of the methods separately for the Augusts and Decembers, since seasonality patterns occur; (ii) the similarities (or differences) between the ETo and MODIS ET estimates despite their different physical meaning. The limitations and future research follow, and at the end the main conclusions are deduced.

Methods
Reference evapotranspiration (ETo) represents "the evaporation of a surface of green grass of uniform height (0.12 m), growing with adequate water" [21]. The estimations of ETo have been produced by three methods, namely, the physically based method FAO-56 Penman-Monteith (FAO PM) [21], which has been used as the reference method, and the temperature-based methods Hansen [53], and Hargreaves-Samani (HS) [54]. Those methods are well established in literature and have been previously used for other areas in Greece with satisfying performance [17,32,45,55,56]. Another reason for the selection of these methods is that they demand as inputs parameters that are available (from groundbased or remote sensed datasets) or that could be estimated via the available parameters according to FAO guidelines [21].

FAO 56 Penman-Monteith
The method used as reference is FAO PM (Equation (1)) [21]. FAO PM employs air temperature (T), net radiation (Rn), vapor pressure and wind speed at 2 m from soil surface (u 2 ). Despite Rs, Rn at a given location depends on latitude, Julian day, albedo (a), cloudiness and elevation (z). The parameters ∆, γ, e s -e a are functions of T and/or pressure (P), with pressure being a function of z (Table 1).

Hansen Equation
Hansen method is a simple method with Rs and T inputs (Equation (2)) [53], which has been previously applied in Greece [34].

Hargreaves-Samani Equation
According to the literature, Hargreaves-Samani (HS) [54] is a simple method that can be applied with only measured temperature data and indirectly incorporates humidity. The estimates for periods of more than 5 days compare favorably with the estimates by FAO PM [54,57]. The HS method has been satisfactorily previously used for areas in Greece and neighboring countries [17,32,55]. The implemented Equation (5) is produced from Equations (4) and (11) ( Table 1) for KT = 0.17 [58].

Data and Models of the Three Methods
The preparation of the inputs has been made in two sections. The first section is the calculation of the equations presented in Table 1 using the ground-based measurements of wind speed (u) measured at different heights from soil surface, Tmin and Tmax, as inputs. These equations (Table 1) are provided by FAO guidelines to compute missing or intermediate parameters [21]. Ground-based data were available by the National Observatory of Athens (NOA) (https://meteosearch.meteo.gr/ (accessed on 15 January 2021)) for sixty-two meteorological stations for the study area (Table A1) in daily scale. For the computation of the missing data equations shown in Table 1 air temperature values (Tmin, Tmax, Tmean) and wind speed (uh) derived from 62 meteorological stations under NOA were used (Table A1). Then all the missing parameters (Table 1) were calculated for every station in daily scale, for every December and August of the 4 study years. Then the daily values were averaged for every station (for every month). Those average values (Rs, Rn, es-ea and u2) were then interpolated over the study area using ordinary kriging, after establishing normality of the data (via log transformation) where needed, in ArcMap 10.6 (https://www.esri.com (accessed on 28 March 2021)). For example, Rs was missing data, computed via Equation (10) ( Table 1), where daily T range from meteorological stations was used to account for cloudiness. Solar radiation (Rs) has been computed the same way [21,26,58] for all the three methods used in the study. For each one of the 8 examined months, the sixty-two (daily averaged) Rs values of the stations were interpolated (ordinary kriging) for the study area and used as input to the corresponding models (see Supplementary Materials). Rn was computed from Rs via Equations (11)-(14) ( Table 1).
Each equation (FAO PM, HS and Hansen) produces the distribution over the study area of the mean daily ETo (mm) of one month (e.g., mean daily  [17], who proved that those parameters can be satisfactorily used that way for Greece. The rest inputs are some of the interpolated parameters (Rs, Rn, e s -e a and u 2 ) from the first section, based on each equation. Models include, amongst others, a sequence of intermediate stages where all needed preparation is performed (e.g., the user can define reprojections, cell size, mask etc.; see Supplementary Materials) via interactive ArcToolbox commands. All intermediate equations (e.g., ∆, γ functions needed in FAO PM formula) are executed via raster calculator. At the final stage of the algorithm raster calculator computes the ETo formula of each method, producing the mean daily ETo distribution of the examined month.
LST day and night products for 2016-2019 (16 monthly scale images) were derived from FORTH (http://rslab.gr/downloads_lst.html (accessed on 9 March 2021)). August and December have been selected for years 2016-2019 as representative months for summer and winter based on the meteorological conditions of the study area (e.g., in December there is less snowfall than in January) and data availability for the majority of the stations (almost no missing data for August compared to June and July). Moreover, according to FAO, during those months the soil is already warm or cold, respectively, thus soil heat flux beneath the grass reference surface is relatively small, so that it can be ignored G day ∼ = 0 [21]. Ground-based data are available online by meteosearch.gr database (https://meteosearch.meteo.gr (accessed on 15 January 2021)) for all automatic stations under NOA.

Statistical Models
Descriptive statistics (i.e., areal values of max, min, mean, std) of each ETo distribution have been calculated in ArcMap 10.6 (via calculate statistics). In order to investigate the agreement between different methods, statistic indices, namely, mean bias (MB), normalized mean bias (NMB) and root-mean-square deviation (RMSD) [59,60] were computed via developing the corresponding models in model builder. These models execute the equations of each statistical index presented in Table 2. Despite that, they include steps of data preprocessing preparation (e.g., reprojections, masks) aiming at exactly the same cells of the participating images to have valid values in order for accurate cell-to-cell computations to take place. Raster calculator is used for all the intermediate equations as well.
For every month (e.g., August 2019) we have run the 3 statistical models for all possible pairs of methods (e.g., RMSD values between the following: HS and FAO PM, Hansen and FAO PM, Hansen and HS). Then, we created the graph of every statistical index separately for Augusts and Decembers during 2016-2019 (e.g., graph of RMSD values for Augusts 2016-2019) to observe their variability with years.

MODIS ET Products
MODIS products of net ET (MOD16A2V6) have been acquired (32 images) from EarthExplorer and EARTHDATA platforms (https://earthexplorer.usgs.gov/ (accessed on 11 February 2021); https://urs.earthdata.nasa.gov/ (accessed on 11 February 2021)). Evapotranspiration MOD16A2V6 products are not raw data, but are produced based on a recently updated complex algorithm [61] based on Penman-Monteith equation, which employs reanalysis meteorological data (humidity, air temperature, radiation) from NASA's Global Modeling and Assimilation Office (GMAO, v.4.0.0) and remote sensed vegetation data. Each one of the 32 images constitutes an 8-day composite of ET and includes, amongst other layers, net ET (cell size of 500 m). The images were then appropriately prepared in ArcMap 10.6 (reprojection, masked with the outline of the study area, cell size of 1000 to be compatible with the distributions of ETo distributions (see Section 2.2)). Then the images of each month were averaged (using cell statistics (mean) command, set to ignore no data pixels). Thus, the mean daily MODIS ET (mm) distributions of the study area for Augusts and Decembers of years 2016-2019 were produced.
MODIS ET is an estimate of actual ET, whereas the three methods produce ETo. Those parameters have different physical meanings and different responses to weather conditions (e.g., precipitation, wind speed). This study is the first part of a thorough research project investigating the agreement between several methods estimating ETo or ETa and the differences among the aforementioned estimates to MODIS ET estimates, and to pan evaporation measurements. This methodology presents novelty in the way methods are executed, via ArcMap models with an interactive character of letting users adjust properly the settings of their study area. They would be available for usage by other researchers and they are prone to creative expansion by researchers (both in python and in model builder environment), which makes them flexible and user-friendly. Another novelty is the creation of statistical models of RMSD, MB and NMB providing the areal value. Those models are developed in a way that they can compute indices by processing large images (with more than 20,000 values) cell-to-cell, instead of the time-consuming process of extracting points from the distributions in ArcMap and using external statistical software to compute the areal index in a series of steps, which demands expertise.

Study Area
Peloponnese is a peninsula to the south of Greece, occupying 21,439 km 2 , which is about 1/6 of the Greek territory, separated from the mainland, with a population of 1,086,935 (census 2011; https://www.statistics.gr/el/statistics/-/publication/SAM03/20 11 (accessed on 20 March 2021)). A large part is covered by high hills and mountains running NW to SE-the highest is Taygetus (2407 m, S). The wider plain (Ilia plain) lies over the west coastal part (Lappa to Kyparissia area). The biggest city is Patra at the north edge. The hydrographic network is well developed [62], though with few large rivers (e.g., Alfeios, Evrotas, Pineios). The main land use/land cover (LULC) types are agricultural patterns, forest and transitional vegetation (Figure 1a), a large part of which (Western Peloponnese) had been subjected to the most severe wildfires of Greece in 2007 [63].
Lithology, climate, and tectonic activity define the relief formation of the area (Figure 1b) [64]. The bedrock of Peloponnese consists of the formations of geotectonic zones (west to east), as follows: the Ionian zone, Gavrovo-Tripolis zone, Olonos-Pindos zone (External Hellenides) and Pelagonian zone (Internal Hellenides) (Figure 1b). The climate of Peloponnese is a Mediterranean warm temperate climate with a dry summer (Csa) [65,66]. The annual average temperature ranges from below 8 • C to 20 • C, precipitation ranges between 400 to over 2000 mm, and sunshine between 1900 and 3100 h (http://climatlas.hnms.gr/sdi/?lang=EN (accessed on 20 March 2021)). However, during the years 2016-2019, which are four out of the five of the warmest years since the pre-industrial era [31], there have been substantial departures from the average T (https://meteosearch.meteo.gr/ (accessed on 15 January 2021)).

Statistical Measures
There are specific statistical indices used in the literature for the assessment of the agreement between the different methods of ET estimation [22,47,59,60]. According to the literature, the coefficient of determination (R 2 ) has been proven as oversensitive to extreme values and insensitive to additive and proportional differences between the compared values, thus it is characterized as a rather misleading goodness of fit for hydrological parameters [69]. For these reasons, the root-mean-square difference (RMSD), mean bias (MB) and normalized mean bias (NMB) were used, which are satisfactorily reliable [70]. The rootmean-square difference expresses the average difference between the compared values of ET (in mm). The mean bias (MB) demonstrates the bias or the mean bias error (MBE) (in mm) ( Table 2). The normalized mean bias (NMB) or normalized mean bias error (NMBE) ( Table 2) are both dimensionless measures, so that direct comparisons of ET values between different time scales are possible. Aiming to calculate RMSD, MB and NMB values between pairs of methods, models have been developed in ArcMap 10.6 (model builder). Those models are easy and quick to run, while they process raster images (213*218) with more than 20,000 (valid) pixels. The plots display MB and NMB rather than MBE and NMBE, respectively, since they have the same absolute value, but the former are signed. A plus sign (+) expresses the overestimation of the method compared to the reference, whereas a minus sign (−) expresses underestimation. The model of RMSD results in a raster of squared differences (RSMD 2 ) calculated cell-to-cell and their mean value, the square root of which is RSMD. The MB model results in a raster of differences (cell-to-cell), the mean value of which represents MB. The model producing NMB masks the denominator raster, using the nominator raster as a mask in order to acquire the same size (213*218) and exactly the same no data pixels as the nominator raster.

Statistical Measures
There are specific statistical indices used in the literature for the assessment of the agreement between the different methods of ET estimation [22,47,59,60]. According to the literature, the coefficient of determination (R 2 ) has been proven as oversensitive to extreme values and insensitive to additive and proportional differences between the compared values, thus it is characterized as a rather misleading goodness of fit for hydrological parameters [69]. For these reasons, the root-mean-square difference (RMSD), mean bias (MB) and normalized mean bias (NMB) were used, which are satisfactorily reliable [70]. The root-mean-square difference expresses the average difference between the compared values of ET (in mm). The mean bias (MB) demonstrates the bias or the mean bias error (MBE) (in mm) ( Table 2). The normalized mean bias (NMB) or normalized mean bias error (NMBE) ( Table 2) are both dimensionless measures, so that direct comparisons of ET values between different time scales are possible. Aiming to calculate RMSD, MB and NMB values between pairs of methods, models have been developed in ArcMap 10.6 (model builder). Those models are easy and quick to run, while they process raster images (213*218) with more than 20,000 (valid) pixels. The plots display MB and NMB rather than MBE and NMBE, respectively, since they have the same absolute value, but the former are signed. A plus sign (+) expresses the overestimation of the method compared to the reference, whereas a minus sign (−) expresses underestimation. The model of RMSD results in a raster of squared differences (RSMD 2 ) calculated cell-to-cell and their mean value, the square root of which is RSMD. The MB model results in a raster of differences (cell-to-cell), the mean value of which represents MB. The model producing NMB masks the denominator raster, using the nominator raster as a mask in order to acquire the same size (213*218) and exactly the same no data pixels as the nominator raster.

Descriptive Statistics (Mean, SD) of Areal Daily ETo and MODIS ET
As presented in Table 3, the daily mean ETo estimates by FAO PM are greater, followed for December by the Hansen's and HS's estimates. The MODIS ET mean values are closer to the Hansen mean values. For August, FAO PM produces higher estimates followed by HS and Hansen. The MODIS ET values are substantially lower than empirical methods' estimates, and similar to the FAO PM estimates for December (Table 3).

Statistical Measures between Estimates by Empirical Methods
It is obvious by the statistical indices (RMSD, MB and NMB) that Hansen and HS display RMSD ≤ 0.11 mm in the daily mean ETo estimation of December for the 4 years ( Figure 4). Concerning the reference method (FAO PM), Hansen is closer with the latter underestimating ETo (RMSD between 0.41 and 0.67 mm d −1 ), whereas the corresponding RMSD for HS lies within 0.50-0.74 mm d −1 , with the maximum daily mean ETo appearing for December 2017. For August, the statistical indices indicate that the HS and FAO PM estimates are consistently close, with the former underestimating ETo with an RMSD between 0.72 and 0.76 mm d −1 (Figure 5). It is obvious that for August 2017 the estimation by Hansen shows significant departure, thus explaining the differentiation between Hansen and HS, which generally display close estimations (RMSD ≤ 0.66 mm d −1 ) for the Augusts (slight underestimation by Hansen) of the other 3 years.

Spatial Distributions of Daily Mean MODIS ET for Decembers and Augusts
As shown in Figure 6a

Statistical Measures between Estimates by Empirical Methods
It is obvious by the statistical indices (RMSD, MB and NMB) that Hansen and HS display RMSD ≤ 0.11 mm in the daily mean ETo estimation of December for the 4 years ( Figure 4). Concerning the reference method (FAO PM), Hansen is closer with the latter underestimating ETo (RMSD between 0.41 and 0.67 mm d −1 ), whereas the corresponding RMSD for HS lies within 0.50-0.74 mm d −1 , with the maximum daily mean ETo appearing for December 2017. For August, the statistical indices indicate that the HS and FAO PM estimates are consistently close, with the former underestimating ETo with an RMSD between 0.72 and 0.76 mm d −1 (Figure 5). It is obvious that for August 2017 the estimation by Hansen shows significant departure, thus explaining the differentiation between Hansen and HS, which generally display close estimations (RMSD ≤ 0.66 mm d −1 ) for the Augusts (slight underestimation by Hansen) of the other 3 years.

Spatial Distributions of Daily Mean MODIS ET for Decembers and Augusts
As shown in Figure 6a

Spatial Distributions of Daily Mean MODIS ET for Decembers and Augusts
As shown in Figure 6a

Statistical Measures for Investigating the Difference between Empirical ETo and MODIS ET Estimates for Decembers and Augusts
As shown in Figure

Parameters Differentiating the Distributions of Reference Evapotranspiration (ETo) and MODIS ET
For December 2016 (see Figure 2a), although all of the methods assign the lowest values to the northcentral part the physically based FAO PM differentiated from HS and Hansen by assigning the maximum values in the following two patterns: where the maximum u2 (>12.50 ms −1 ) and elevated Rn occurs (eastern part), or when the maximum Rn (>2.53 MJ m −2 ) and elevated u2 (>7.30 ms −1 ) is encountered (SE part). For December 2017 (see Figure 2b), Hansen and HS demonstrate almost identical distributions. FAO PM differentiates the northcentral edge due to a very high mean monthly u2 (>9.50 ms −1 ). The sensitivity of FAO PM towards u2 and Rn is responsible for a maximum of 1 mm greater than those of the rest methods assigned to the area around Argos, where VPD is the maximum (0.48 kPa) (u2: 7.99-13.46 ms −1 , Rn: 2.56-2.63 MJ m −2 and VPD: 0.30-0.35 kPa). For December 2018 (see Figure 2c), Hansen and HS exhibit similar contributions and a similar range of values. The main difference is the green area, which FAO PM assigns at the northwestern part, which is due to the very low monthly mean u2 reported at the station (<0.09 ms −1 ) since the other parameters are rather elevated (VPD = 0.38 kPa and Rn = 2.17 MJ m −2 ), expressing the sensitivity of FAO PM to wind speed (u2) not only for the maximum, but also for the minimum ETo values assignments. For December 2019 (see Figure 2d), FAO PM displays the maximum due to increased Rn (>2.61 MJ m −2 ) and u2 (>7.78 ms −1 ). The high u2 (>9.65 ms −1 ) at the northcentral edge is responsible for the yellow color in the FAO PM distribution in contradiction with the green color assigned by the other methods that do not incorporate u2.

Parameters Differentiating the Distributions of Reference Evapotranspiration (ETo) and MODIS ET
For December 2016 (see Figure 2a), although all of the methods assign the lowest values to the northcentral part the physically based FAO PM differentiated from HS and Hansen by assigning the maximum values in the following two patterns: where the maximum u2 (>12.50 ms −1 ) and elevated Rn occurs (eastern part), or when the maximum Rn (>2.53 MJ m −2 ) and elevated u2 (>7.30 ms −1 ) is encountered (SE part). For December 2017 (see Figure 2b), Hansen and HS demonstrate almost identical distributions. FAO PM differentiates the northcentral edge due to a very high mean monthly u2 (>9.50 ms −1 ). The sensitivity of FAO PM towards u2 and Rn is responsible for a maximum of 1 mm greater than those of the rest methods assigned to the area around Argos, where VPD is the maximum (0.48 kPa) (u2: 7.99-13.46 ms −1 , Rn: 2.56-2.63 MJ m −2 and VPD: 0.30-0.35 kPa). For December 2018 (see Figure 2c), Hansen and HS exhibit similar contributions and a similar range of values. The main difference is the green area, which FAO PM assigns at the northwestern part, which is due to the very low monthly mean u2 reported at the station (<0.09 ms −1 ) since the other parameters are rather elevated (VPD = 0.38 kPa and Rn = 2.17 MJ m −2 ), expressing the sensitivity of FAO PM to wind speed (u2) not only for the maximum, but also for the minimum ETo values assignments. For December 2019 (see Figure 2d), FAO PM displays the maximum due to increased Rn (>2.61 MJ m −2 ) and u2 (>7.78 ms −1 ). The high u2 (>9.65 ms −1 ) at the northcentral edge is responsible for the yellow color in the FAO PM distribution in contradiction with the green color assigned by the other methods that do not incorporate u2. Throughout the years, August 2016 exhibits higher maximum values at the southmost part (around Geraki). Regarding the distribution by FAO PM, it is obvious that the green central lane from the northmost to southmost edges in 2016 is eventually diminished to just the central mountainous part in 2019. This spreading of the redish areas with the years can also be noticed from the distributions of the other two methods with a lower consistency over the years. Our findings are in line with the seasonality, which makes empirical methods deviate more during the summertime, as several researchers have noted not only for the South Mediterranean areas in general, but also specifically for Greece [47]. Both Hansen and HS exhibit good performance (low RMSD) in the winter season, but poor performance in the summer season. Their performance has been reported to be the opposite for the winter and summer of neighboring areas with humid climates (e.g., Serbia [55])

Differences and Similarities between Estimates of Reference Evapotranspiration (ETo) and MODIS ET
MODIS ET is a product of a sophisticated algorithm; however, it produces considerable underestimation in the literature (such as by 26%; [72]; and by 25% uncertainty [73]) compared with sites of the AmeriFlux network, which has been used for the calibration of the updated algorithm [61] attributed to the need of further refinement. Besides, global satellite products generally contain more noise than ground-based measurements lacking in accuracy [33,74]. Westerhoff (2015) [74] has conservatively set the uncertainty of satellite to ground-based estimations equal to 1.5. For December, FAO PM showed the greatest difference to MOD16A2V6 for all of the statistical indices (e.g., RMSD: 0.52-0.84 mm d −1 , MB: 0.47-0.78 mm d −1 ). HS appeared to have a minimal difference to MODIS ET (e.g., RMSD: 0.15-0.17 mm d −1 , MB: -0.02-0.09 mm d −1 ) and Hansen followed (RMSD: 0.16-0.22 mm d −1 , MB: 0.03-0.18 mm d −1 ). Moreover, all of the empirical methods compared to MODIS ET, except from Hansen, exhibited the greatest departures (maximums) in statistical indices (demonstrating lowest agreement) for 2017. Relative studies in Greece report accuracy of empirical methods' performances during the winter. Aiming at investigating whether good performance is a seasonal characteristic also for 2016-2019, the research should be expanded to the rest of the months of this period. The estimates of ETo are close to MODIS ET in the Decembers of 2016-2019, because although the former reflects the evaporative demand of the atmosphere, the latter is an estimate of actual ET subjected to water stress. However, during December the evaporative demand of the atmosphere is low, the VPD is low (thus relative humidity is elevated), and December is the month with the third highest precipitation (in mm) in Greece, so the evaporative demand is satisfied.
The opposite happens during the Augusts when the evaporative demand is elevated (high Rn, T, VPD and low u 2 ), but due to water stress during the summertime it cannot be adequately satisfied. Thus, the MODIS ET values are considerably lower than the ETo estimates, with the former appearing values almost at the levels of December (i.e., FAO PM values for December). FAO PM differs the most from MODIS (as also in December) (RMSD: 3.87-4.06 mm d −1 , MB: 3.26-4.06 mm d −1 ) due to its incorporation of more climate parameters, and the Hansen estimates are the closest to the MODIS ones (RMSD: 1.51-2.91 mm d −1 , MB: 0.96-1.71 mm d −1 ) followed by HS; however, with considerable discrepancies in both cases. The Hansen ETo values exhibit the greatest difference from the MODIS actual ET estimates for the two warmest years (2016 and 2019). Considering the physical meaning of ETo and MODIS ET, this might imply Hansen being more affected by Rs (and subsequently by T) than the rest. Another observation is that Hansen bears similar departures in the estimates to MODIS ET for August 2017, which should be further investigated in terms of sensitivity.
The main limitations are the following: Rs is not measured during the study years but computed under the FAO guidelines for missing data, with a KT (0.17) coefficient used for both coastal and mainland areas. Moreover, the constant value (0.23) used for albedo, the replacement of which is with MODIS albedo products, would refine the estimates. Model builder in ArcMap cannot process raster images and point features together, thus the RMSD model bears the limitation of producing the squared areal RMSD value (RMSD 2 ). However, the square root of a number can be trivially computed.
The primary contribution of the study is that the application of ETo models and statistical models addresses international and interdisciplinary interest. Moreover, the produced maps of ETo and MODIS ET of Peloponnese can be used for comparisons against the local ETo and ETa estimates produced by other methods or from data acquired from alternative sources. In particular, FAO PM distributions could be employed as a reference by researchers. Moreover, considering their physical meaning, ETo and MODIS ET maps of the most warm recent years provide water managers of Peloponnese with a general picture, which can be integrated to define the irrigation needs during the Augusts [34,51], and the flooding risks during the Decembers [50].

Conclusions
Our findings are in line with the seasonality, which makes the empirical methods more accurate in the winter (December) and deviate more during the summertime (August), both for the South Mediterranean areas and specifically for Greece. For the Decembers of 2016-2019, the estimates by HS and Hansen are almost identical (NMB up to 0.11 mm), but FAO PM produces higher estimates attributed to the wind speed (u 2 ) sensitivity of the latter. For August 2016-2019, HS and FAO PM present almost identical estimates, probably because the wind speed values can be neglected for the majority of the cells. The Hansen method produces departures (lower values) for August 2017 compared to the rest. These departures are also observed for MODIS ET, thus the relations would be investigated in terms of physical or computational factors, which affect alike the sensitivity of the two different approaches.
The areal daily mean MODIS ET, as an estimate of the actual ET, was examined along with the methods' ETo estimates. Since they have different physical meanings, the ETo values for the Augusts 2016-2019 were much higher than MODIS anticipated. However, for the Decembers 2016-2019, HS produced almost identical estimates (NMB between −0.02 and 0.11) with the MODIS ET averaged products. This is attributed to the low evaporative demand during December, which can be satisfied by rainfall. HS could be used satisfactorily regarding the accuracy for the estimation of the daily mean actual ET of Peloponnese for the Decembers 2016-2019. The latter could be further investigated for a larger period of years for the rest of the winter months.
Partitioning Peloponnese into rough segments with similar climate characteristics and LULC, in order to calibrate the empirical methods for the local and regional spatial scales for the monthly time scale, is proposed for future research, since it could assess the performance of empirical methods during different seasons and LULC types. Moreover, examination of ETo along with ETa, MODIS ET and Epan would contribute to the perception of their complex relations despite their substantial differences.