An Intercomparison of Satellite-Based Daily Evapotranspiration Estimates under Different Eco-Climatic Regions in South Africa

Knowledge of evapotranspiration (ET) is essential for enhancing our understanding of the hydrological cycle, as well as for managing water resources, particularly in semi-arid regions. Remote sensing offers a comprehensive means of monitoring this phenomenon at different spatial and temporal intervals. Currently, several satellite methods exist and are used to assess ET at various spatial and temporal resolutions with various degrees of accuracy and precision. This research investigated the performance of three satellite-based ET algorithms and two global products, namely land surface temperature/vegetation index (TsVI), Penman–Monteith (PM), and the Meteosat Second Generation ET (MET) and the Global Land-surface Evaporation: the Amsterdam Methodology (GLEAM) global products, in two eco-regions of South Africa. Daily ET derived from the eddy covariance system from Skukuza, a sub-tropical, savanna biome, and large aperture boundary layer scintillometer system in Elandsberg, a Mediterranean, fynbos biome, during the dry and wet seasons, were used to evaluate the models. Low coefficients of determination (R2) of between 0 and 0.45 were recorded on both sites, during both seasons. Although PM performed best during periods of high ET at both sites, results show it was outperformed by other models during low ET times. TsVI and MET were similarly accurate in the dry season in Skukuza, as GLEAM was the most accurate in Elandsberg during the wet season. The conclusion is that none of the models performed well, as shown by low R2 and high errors in all the models. In essence, our results conclude that further investigation of the PM model is possible to improve its estimation of low ET measurements.


Introduction
As an essential climate variable (ECV), evapotranspiration (ET) plays a critical role as a link of the energy, carbon and water cycles, with latent heat of vaporization serving as the largest sink of heat in the atmosphere.It is, therefore, important for disciplines such as agriculture, hydrology, meteorology and climatology.Because of its high spatio-temporal variability, it is a challenge to directly measure this biophysical variable.Remote sensing remains the only feasible means of spatially estimating ET over varying spatial and temporal extents.Several authors have reviewed the remote sensing approaches used to estimate ET [1][2][3][4] based on their structural complexities, theories and underlying assumptions, parameterizations, and uncertainties and limitations, and classified them as: (i) empirical methods involving the use of statistically-derived relationships between ET and vegetation indices such as the normalized difference vegetation index (NDVI) or the enhanced vegetation index (EVI) [5][6][7]; (ii) residual surface energy balance modeling (single-and dual-source models), which include the Surface Energy Balance Algorithm over Land (SEBAL)/Mapping EvapoTranspiration at high Resolution with Internalized Calibration (METRIC) [8][9][10][11], Surface Energy Balance System (SEBS) [12,13]; (iii) physically-based methods involving the Penman-Monteith (PM) [14][15][16] and Priestley-Taylor (PT) equations [17][18][19]; and (iv) data assimilation methods applied to the heat diffusion equation and radiometric surface temperature sequences.Global satellite-based ET products have been produced using these algorithms, the MOD16 ET product that is estimated using the Penman-Monteith equation by Mu et al. [20] and Mu et al. [21]; the Meteosat Second Generation ET product (MET) derived using the physically-based Soil-Vegetation-Atmosphere-Transfer (SVAT) model Tiled ECMWF (European Centre for Medium-Range Weather Forecasts) Surface Scheme of Exchange processes at the Land surface (TESSEL) that uses a combination of atmospheric model outputs and Meteosat Second Generation's Spinning Enhanced Visible and Infrared Imager MSG-SEVIRI remote sensing data [22][23][24]; and the Global Land-surface Evaporation: the Amsterdam Methodology (GLEAM) based on the Priestley-Taylor equation and the Gash analytical model of forest rainfall interception [25,26].
Validation and model comparison studies have been done at different locations, time steps and periods.Yang et al. [27] examined the performances and model physics of three dual-source ET models, the Hybrid dual-source scheme and Trapezoid framework-based ET Model (HTEM), the Two-Source Energy Balance (TSEB) model, and the MOD16 ET (PM) algorithm in the Heihe River Basin in Northwest China.They reported that HTEM outperformed the other models, whereas PM performed worst; and the reason for poor performance of PM was the model does not effectively capture soil moisture restriction on ET.Singh and Senay [28] compared four models, i.e., METRIC, SEBAL, SEBS, and the Operational Simplified Surface Energy Balance (SSEBop) models, at three AmeriFlux cropland sites in the Midwestern United States.Their results showed all models recorded R 2 of 0.81, with METRIC and SSEBop having low RMSE (<0.93 mm/day) and a high Nash-Sutcliffe coefficient of efficiency (>0.80), whereas SEBAL and SEBS gave relatively higher bias.Bhattarai et al. [29] evaluated five models (SEBAL, METRIC, SSEBop, SEBS and Simplified Surface Energy Balance Index (S-SEBI)) at four sites (marsh, grass, and citrus surfaces), to identify the most appropriate for use in the humid southeastern United States.SEBS generally outperformed the other with RMSE of 0.74 mm/day, whereas SSEBop was consistently the worst performing model (RMSE = 1.67 mm/day).They observed that for short grass conditions, SEBAL, METRIC, and S-SEBI worked much better than SEBS.Ershadi et al. [30] assessed the performance of four models, i.e., SEBS, PM, Priestley-Taylor Jet Propulsion Laboratory (PT-JPL), and Advection-Aridity (AA), in twenty FLUXNET towers covering different biomes that included grassland, cropland, deciduous broadleaf forest and evergreen needleleaf forest, mainly in Europe and North America, at half-hourly, hourly and monthly time steps.Their results showed that PT-JPL outperformed the other models, followed by SEBS, PM, and lastly AA.Their overall findings, however, were that there was no model that consistently performed well across all the biomes.A study by Ha et al. [31] in semi-arid pine forests with variable disturbance history, over a period of four years and at monthly time steps, also showed that the PT model gave the best results, with the PM model and MOD16 ET product under predicting ET at all sites.Vinukollu et al. [32] tested SEBS, PT and PM over 16 FLUXNET sites and concluded that PT outperformed the other models.
Many remote sensing-based ET estimation studies have been performed in different South African landscapes and land uses, as reviewed by Gibson et al. [3].The review stresses the importance of validating the remotely sensed ET estimates to allow for confidence in their use and application in the various biomes.Jarmain et al. [33] used a large discontinuous in situ ET dataset, mainly from agricultural land, to evaluate remote sensing-based ET algorithms (SEBS, SEBAL, METRIC and VIIT), and recorded poor performance of all the models.Evaporative fraction estimation by these models was the main source of error, which led to low accuracy in ET estimation.One of the challenges of this study was the limited data points at each site for substantive statistical analysis, hence, no distinct conclusion could be made.Jovanovic et al. [34] reported an R 2 of 0.72 and 0.75 in their validation of the 30 min and daily MET ET products for the fynbos vegetation of the Riverlands Nature Reserve.Ramoelo et al. [35] validated the PM-derived MOD16 8-day ET product using multiyear eddy covariance-derived ET datasets for two flux towers in savannahs, Skukuza and Malopeni.Inconsistent results were attributed to various factors, including the parameterization of the PM model, flux tower measurement errors, and flux tower footprint vs. MODIS pixel size.Furthermore, Sun et al. [36] used the Skukuza flux tower data to evaluate a remote sensing-based continental ET product.The results were reasonable during the wet season, whereas low coefficients of determination were observed in the dry season.One example of remotely sensed ET applications in South Africa is the use of SEBAL by the eLEAF company, in collaboration with the Water Research Commission (WRC), the South African's Department of Agriculture and academic institutions, to provide information on water use efficiency of irrigated crops, inclusive of grapes, deciduous fruits, sugarcane and grain crops in the Western Cape Province of South Africa [37][38][39].Overall, limited work has been done to assess and compare the performance of different ET models in different SA natural ecosystems.
There is still scope to extensively compare different models for different biomes in semi-arid ecosystems, which would result in the identification of the most accurate and robust model that could be used to map and monitor ET at national scale.Hence, this study intended to evaluate and compare the performance of daily ET estimates derived using TsVI and MOD16 Penman-Monteith-based models, and GLEAM and MET global products, under two different climatic regions, a subtropical, savannah, summer rainfall and a Mediterranean, fynbos, winter rainfall climates.From each climatic region, rainy and dry periods were selected to evaluate the performance of each model.

Materials and Methods
One of the challenges in evaluating evapotranspiration models is the availability of accurate and complete datasets.The comparisons were performed for two eco-regions in South Africa, Skukuza in the northeast of South Africa, an area characterized by summer rainfall, and Elandsberg in the southwest of the country, characterized by winter rainfall (Figure 1).A full description of the sites is given in Section 2.1 below.

Summer Rainfall Savannas and Skukuza Flux Tower
The Skukuza eddy covariance flux tower (25.02 • S, 31.50 • E) (Figure 2) was established early 2000 as part of the SAFARI 2000 experiment, set up to understand the interactions between the atmosphere and the land surface in southern Africa [40,41].
The site is located in the Kruger National Park (South Africa) at 365 m above sea level, and receives 550 ± 160 mm precipitation per annum between November and April, with significant inter-annual variability.The soils are generally shallow, with coarse sandy to sandy loam textures.

Eddy Covariance System
Since 2000, ecosystem-level fluxes of water, heat and carbon dioxide are measured using an eddy covariance system mounted at 16 m of the 22 m high flux tower.The measurements taken and the instruments used are summarized in Table 1.Ancillary meteorological measurements include air temperature and relative humidity, also measured at 16 m, using a Campbell Scientific HMP50 (Loughborough, UK) probe; precipitation at the top of the tower using a Texas TR525M (Dallas, TX, USA) tipping bucket rain gauge; wind speed and direction using a Climatronics Wind Sensor (Bohemia, NY, USA); and soil temperature using Campbell Scientific 107 (Logan, UT, USA) soil temperature probe.
Flux footprint modeling for the Skukuza tower was done using footprint models that are incorporated in the EddyPro software, i.e., the footprint parameterization model of Kljun et al. [42] and footprint model of Kormann and Meixner [43].Depending on a number of factors, including wind direction and friction velocity, atmospheric stability, terrain homogeneity and measurement height, either of the models are selected to compute the footprint.It was estimated that 90% of the fluxes originated from an average of 1.6 km upwind from the flux tower, as reported by Ramoelo et al. [35].The shape of the flux source area is dependent on the wind velocity and direction, and varies throughout the day.With a surface energy balance closure of 1.03, the LE fluxes were then corrected by: where LEcorr is corrected latent heat flux, and Res is the residual (Res = Rn − G − H − LE).Along the 900 m transect of the scintillometer, the average vegetation height is 1 m, with general land cover being shrubland and low fynbos.

Scintillometry System
Prior to selecting the Elandsberg site for LAS installation and ET validation, a study was conducted to ensure that the area selected meets a number of criteria as a suitable flux tower site, including land cover homogeneity (single land cover type within the image pixel, in this case the MODIS pixel), vegetation height, topographic variability, and atmospheric stability [44].The site for the scintillometer study of evapotranspiration was, therefore, selected based on the extent and uniformity of the natural vegetation, and also to ensure to ensure that the match between the selected pixel and the scintillometer transect would be accurate enough for the ground-truthing of ET from that pixel.
A mobile large aperture boundary layer scintillometer (LAS) (BLS 900, Scintec, AG, Rottenburg, Germany) was installed in Elandsberg Nature Reserve in October 2012 and collected data until November 2013.The scintillometer transmitter (Figure 3   Sensible heat flux was calculated from the changes in the refractive index of air between the scintillometer transmitter of monochromatic infrared radiation at 880 nm and the receiver.The net radiation was measured using a net radiometer (Kipp and Zonen, Delft, The Netherlands), installed at 2 m above the vegetation in the middle of the scintillometer transect.The soil heat flux (G) was measured using a cluster of four soil heat flux plates (REBS, Inc., Seattle, WA, USA), installed at a depth of 80 mm at various positions within the MODIS pixel.
An automatic weather station was also installed to record air temperature, relative humidity, wind speed and direction, solar irradiance, and rainfall.Temperature and relative humidity were measured at 1.5 m using a CS500 probe (Vaisala, Helsinki, Finland), while windspeed and direction were measured at 2.5 m height using an RM Young wind sentry (Model 03001-Campbell Scientific Ltd., Logan, UT, USA), and solar irradiance was monitored using a pyranometer (Apogee Instruments, Logan, UT, USA).Rainfall was measured using a tipping rain gauge (Model TE 525WS-Campbell Scientific Ltd., Logan, UT, USA).Soil water content was also measured at 30 min intervals in the depth range 0-20 cm using a CS616 capacitance probe (Campbell Scientific Ltd., Logan, UT, USA).All the sensors were connected to a CR23X data logger (Campbell Scientific Ltd., Logan, UT, USA).

TsVI Method
The concept of the surface temperature-vegetation index triangle (TsVI) method was discovered by Goward et al. [45], and has been used to retrieve soil water content, analyze land use/land cover change and monitor droughts with satellite data [11,[46][47][48].Using the TsVI feature space, Jiang and Islam [49][50][51] adapted the PT equation to estimate regional evaporative fraction (EF) and ET.This method calculates ET as a function of available energy, i.e., net radiation minus soil heat flux.Its main assumption is that ET depends on soil moisture and vegetation cover, the method requires a large heterogeneous area with a varied range of values.The PT parameter ϕ, which accounts for the aerodynamic and canopy resistances in the PT formulation, is replaced with φ in the proposed Sensible heat flux was calculated from the changes in the refractive index of air between the scintillometer transmitter of monochromatic infrared radiation at 880 nm and the receiver.The net radiation was measured using a net radiometer (Kipp and Zonen, Delft, The Netherlands), installed at 2 m above the vegetation in the middle of the scintillometer transect.The soil heat flux (G) was measured using a cluster of four soil heat flux plates (REBS, Inc., Seattle, WA, USA), installed at a depth of 80 mm at various positions within the MODIS pixel.
An automatic weather station was also installed to record air temperature, relative humidity, wind speed and direction, solar irradiance, and rainfall.Temperature and relative humidity were measured at 1.5 m using a CS500 probe (Vaisala, Helsinki, Finland), while windspeed and direction were measured at 2.5 m height using an RM Young wind sentry (Model 03001-Campbell Scientific Ltd., Logan, UT, USA), and solar irradiance was monitored using a pyranometer (Apogee Instruments, Logan, UT, USA).Rainfall was measured using a tipping rain gauge (Model TE 525WS-Campbell Scientific Ltd., Logan, UT, USA).Soil water content was also measured at 30 min intervals in the depth range 0-20 cm using a CS616 capacitance probe (Campbell Scientific Ltd., Logan, UT, USA).All the sensors were connected to a CR23X data logger (Campbell Scientific Ltd., Logan, UT, USA).

TsVI Method
The concept of the surface temperature-vegetation index triangle (TsVI) method was discovered by Goward et al. [45], and has been used to retrieve soil water content, analyze land use/land cover change and monitor droughts with satellite data [11,[46][47][48].Using the TsVI feature space, Jiang and Islam [49][50][51] adapted the PT equation to estimate regional evaporative fraction (EF) and ET.This method calculates ET as a function of available energy, i.e., net radiation minus soil heat flux.Its main assumption is that ET depends on soil moisture and vegetation cover, the method requires a large heterogeneous area with a varied range of values.The PT parameter ϕ, which accounts for the aerodynamic and canopy resistances in the PT formulation, is replaced with φ in the proposed formula by Jiang and Islam [49], which is determined using the triangular shape of the TsVI feature space.ET is estimated using: where φ is a combined-effect parameter accounting for aerodynamic resistance (dimensionless), Rn is surface net radiation (W/m 2 ), G is soil heat flux (W/m 2 ), ∆ is slope of saturated vapor pressure versus air temperature (kPa Defining the φ Parameter The TsVI triangle method proposed by Jiang and Islam [49][50][51] has been applied at different scales to estimate the parameter φ in the PT method, and thus evaporative fraction (EF).Modifications have also been made to this method as shown by Wang et al. [52] who combined the triangle method with thermal inertia, and developed a day-night Ts difference-NDVI triangle using MODIS land surface products to estimate EF.Stisen et al. [53] also combined the triangle method with thermal inertia, and developed a quadratic function of NDVI to determine φ at dry edge using MSG-SERVIRI data.Furthermore, Tang et al. [54] replaced NDVI in the construction of the triangle with fractional cover (Fr), and developed an automatic algorithm to determine the wet and dry edges.
The wet and the dry edges of the two-dimensional triangular space for each vegetation class are determined first, where the wet/cold edge represents the potential evapotranspiration, and the dry/warm edge represents water-stressed conditions.This global φ min and φ max are set at zero for a dry bare soil surface, and 1.26 for a saturated or well vegetated surface, respectively.Assuming that φ increases linearly with a decrease in Ts between φ min and φ max for any given Fr, φ min,i is then linearly interpolated between φ min and φ max , with φ i,max = φ max = 1.26 as: where Then, for each pixel in the triangle, φ i is then derived by using the normalized temperature: where T max and T min are the corresponding maximum and minimum surface temperatures at the dry and wet edges, respectively.The evaporative fraction (EF) is then calculated using:

Penman-Monteith Based MOD16 ET Model
The Penman-Monteith (PM) is a physically-based model that incorporates heat and water vapor mass transfer principles, and is known as the combination equation.Developed and modified by Penman [55] and Monteith [56], this algorithm was initially designed as a single-source model, computing ET from a heterogeneous land surface as a single component.Modified versions of this model now include the estimation of different components of water loss from bare soil and canopy intercepted water (evaporation), and transpiration via the canopy.Mu et al. [20] then modified the model by adding vapor pressure deficit and minimum temperature to constrain stomatal conductance, using enhanced vegetation index instead of NDVI to calculate vegetation cover fraction and included the separate calculation of soil evaporation.Further improvements were made to the algorithm, including calculating ET as a sum of day-and night-time ET, adding soil heat flux calculation, separating dry and wet (interception) canopy surfaces, and soil surface into saturated and moist surface, as well as improving stomatal conductance, aerodynamic resistance and boundary layer resistance estimates [21].Taking into account that relative humidity, and, by extension, vapor pressure deficit were soil moisture stress proxies in Mu et al. [20,21], Sun et al. [57] used the soil moisture index estimated from the TsVI method to constraint soil evaporation.
Latent heat flux is estimated as: where λE is the latent heat flux (W/m 2 ); λ is the heat of vaporization (J/kg); s is the slope of the curve relating saturated vapor pressure to temperature (Pa/K); A is available energy (W/m 2 ); ρ is the air density (kg/m 3 ); Cp is the specific heat capacity of air (J/kg/K); γ (Pa/K) is the psychrometric constant; e sat is the saturation vapor pressure (Pa); e is the actual vapor pressure (Pa), where e sat − e = vapor pressure deficit (VPD); r a (s/m) is the aerodynamic resistance; and r s (s/m) is the canopy resistance, which is the reciprocal of canopy conductance gc (gc = 1/rc).The MOD16 remote sensing-based ET algorithm predicts ET globally at 86% accuracy when compared with eddy measurements of ET over many sites in the AmeriFlux network.Building on previous algorithms, it uses a physically based PM approach driven by MODIS-derived vegetation data.ET is calculated as the sum of daytime and nighttime components using vapor pressure deficit and minimum temperature to control stomatal resistance.Stomatal resistance is scaled up to the canopy level using LAI to calculate canopy resistance for plant transpiration.The algorithm also models soil heat flux and separates evaporation from a wet canopy and transpiration from a dry canopy.Actual soil evaporation is also calculated from a potential evaporation value.

Net Radiation and Soil Heat Flux Estimation
Net radiation (Rn) is the difference between incoming and outgoing long-and shortwave radiation fluxes on the Earth's surface.It plays a very important role in the exchange processes of water and heat over the land surface.It is a critical parameter in the estimation of ET, and all ET models require its estimation [11,12,55,58].
It is expressed in terms of its components: where α is the land surface albedo, R ↓ s is the incoming shortwave radiation (W/m 2 ), R ↓ l is the longwave incoming radiation (W/m 2 ), R ↑ l is the outgoing longwave radiation (W/m 2 ), σ is the Stephan-Botlzmann constant (5.670373 × 10 −8 W/m 2 /K 4 ), ε a is the atmospheric emissivity, ε s is the surface emissivity, and T a and T s are air and surface temperature (K), respectively.
For this study, the method of Shine [59] was used to estimate R ↓ s : where S 0 is the Solar constant (1367 W/m 2 ), θ is the Solar zenith angle, and e is the vapor pressure.
Remote Sens. 2017, 9, 307 9 of 21 Soil heat flux (G) was computed as: where Tc is ratio of G to Rn for full vegetation cover, and Ts is ratio of G to Rn for dry bare soil.The result used in the estimation of ET from Rn and G are estimated using the same set of equations and parameters for the PM and TsVI models, hence the study will also compare these parameters.

Global Evapotranspiration Products
Apart from evaluating the two ET methods, the performance of two global ET products, i.e., the Meteo-sat (MET) and the Global Land Evaporation: the Amsterdam Model (GLEAM) ET products, will also be assessed in this study.

LSA SAF Evapotranspiration Product
The global ET product (MET) by the EUMETSAT Satellite Application Facility on Land Surface Analysis (LSA-SAF) based on the SEVIRI sensor on-board the Meteosat Second Generation geostationary satellites (MSG-SEVIRI), is in near-real time, and also available at 30 min and daily time intervals and 3 km spatial resolution.This product is derived using the physically-based TESSEL SVAT scheme that uses a combination of the European Centre for Medium-range Weather Forecasts (ECMWF) atmospheric model outputs (i.e., air and dew point temperature, humidity, wind speed, atmospheric pressure and soil moisture) and MSG/SERVIRI remote sensing data [22][23][24].This model divides the pixel into different tiles representing different land covers within the pixel, with some parameters defined at pixel level and the data are extracted from the ECOCLIMAP database.The pixel ET is the sum of the weighted contribution of each tile.
It can be downloaded from the Land Surface Analysis Satellite Applications Facility (LSA SAF) website, http://landsaf.meteo.pt/.

GLEAM Evapotranspiration Product
The Global Land Evaporation: the Amsterdam Model (GLEAM) ET product is derived using the semi-empirical PT equation and the Gash analytical model of forest rainfall interception [26].The model uses inputs from different satellites to estimate ET daily at 0.25 • spatial resolution.The use of the simple PT equation means there is no parameterization of stomatal and aerodynamic resistances.It estimates different evaporation components, including transpiration for three land cover types, i.e., tall canopies, short vegetation, interception loss, bare-soil evaporation, snow sublimation and open-water evaporation.The potential ET estimates are then constrained by a multiplicative stress factor computed based on the content of water in vegetation and the root zone.The final ET estimate presented in a grid is then given as a weighted average from the three land covers.

Datasets
The three models that were tested in this study require different inputs and parameterizations.Table 2 lists each model's inputs and their sources, i.e., whether satellite or meteorological based.

In Situ and Meteorological Measurements
To select the ET validation periods, eddy covariance data were filtered based periods having extensive data with minimal gaps per day and per month, and the availability of all required input data used in the different models.For each site, two time periods were selected for the summer and the winter season, i.e., wet and dry periods.For Skukuza, the periods 1-31 January 2012 (wet), and 5 May-5 June 2012 (dry) were selected, whereas for Elandsberg, 9 November-9 December 2012 (dry), and 7 June-6 July 2013 (wet) were selected.
Meteorological data from each site, including air temperature (Ta), precipitation, relative humidity (RH), atmospheric pressure (P) and wind speed (u), were used to calibrate the models.From these data, we calculated daily maximum and minimum air temperature, average daily, daytime and nighttime temperatures, daily atmospheric pressure, wind speed and relative humidity.The input variable at each site were measured using the instruments listed in Table 3.The 30 min LE measurements from the Skukuza eddy covariance were converted to ET, which were summed up to daily ET measurements following Mu et al. [21].
From Elandsberg, LE was estimated using the energy balance equation, i.e., as a residual of the surface energy balance from measurements of Rn, G and H (estimated with the LAS), and then converted to daily ET estimates.

Remote Sensing Data
This study made use of MODIS Terra/Aqua products as data inputs.The following MODIS products were used as inputs for this study: daily 1 km MOD11A1 land surface temperature (LST) and emissivity, 8-day 1 km MOD15A2 leaf area index and FPAR, 16-day 1 km MOD13A2 vegetation indices (NDVI and EVI), 16-day 1 km MCD43A3 albedo, and MOD03 geolocation for solar zenith.These datasets were downloaded from the NASA Land Processes Distributed Active Archive Centre (LP DAAC) website (https://lpdaac.usgs.gov/data_access).
For the eight-day datasets (i.e., LAI/FPAR and surface albedo), images for Julian Days Other methods and techniques that compute ET at daily, monthly and seasonal scale have been developed.These include the integration of the feedback method that uses the complimentary relationship between actual ET and pan ET [60], data fusion methods (i.e., combining different data sources) [61,62], and the backward-averaged iterative two-source surface temperature and energy balance solution (BAITSSS) algorithm [63].

Data Analysis
The models described in Section 2.2 were run using the above data to estimate daily ET during the selected times, i.e., 1 and 31 January 2012 (wet), and 5 May to 5 June 2012 (dry) for Skukuza, and 9 November to 9 December 2012 (dry), and 7 June to 6 July 2013 (wet) for Elandsberg.
The accuracy of the ET estimates depends on a number of factors, including the algorithms and their parameterizations and the accuracy of the input datasets.The performance of each model was assessed using the coefficient of determination (R 2 ), a measure of goodness of fit; bias; mean absolute error (MAE); root mean square error (RMSE); and relative RMSE (rRMSE), the latter three being measurements of error/accuracy.RMSE provides information on the short-term performance of a model by allowing a term by term comparison of the actual difference between the predicted value and the measured value, although it does not differentiate between under-and over-estimation; bias provides information on the long-term performance of a model.A positive value gives the average amount of over-estimation in the estimated values and vice versa.

Ground Measurements of Evapotranspiration and Meteorological Input Variables
The first part of the results describes the meteorological input variations, i.e., daily mean air temperature (Ta), RH, Rn and precipitation in relation to ET, between the two sites (Skukuza located in a rainy summer region, and Elandsberg characterized by dry summers).
Daily time series of Ta, RH, Rn, ET and precipitation for Skukuza and Elandsberg for the validation periods are shown in Figure 4.For Skukuza, the daily ET for the rainy period (DOY 1-31) ranged between 0.96 and 6.24 mm/day, with an average of 4.06 mm/day-and a standard deviation of 1.24 mm/day.Temperature varied between 22.16 • C and 30.93 • C, with a mean of 27.25 ± 2.24 • C, RH ranged between 91.54% and 58.27%, with an average of 69.95%, and average Rn for this period was 144.32 ± 43.7 W/m 2 .A total of 279.69 mm precipitation was recorded within the same period.During the dry period (DOY 128-153), daily ET ranged between 0.53 and 1.47 mm/day, with an average of 0.99 ± 0.26 mm/day, Rn was between 40.15 and 91.92 W/m 2 , and mean air temperature of 18.68 ± 1.44 • C was recorded; no precipitation was recorded during the same period.

Ground Measurements of Evapotranspiration and Meteorological Input Variables
The first part of the results describes the meteorological input variations, i.e., daily mean air temperature (Ta), RH, Rn and precipitation in relation to ET, between the two sites (Skukuza located in a rainy summer region, and Elandsberg characterized by dry summers).
Daily time series of Ta, RH, Rn, ET and precipitation for Skukuza and Elandsberg for the validation periods are shown in Figure 4.For Skukuza, the daily ET for the rainy period (DOY 1-31) ranged between 0.96 and 6.24 mm/day, with an average of 4.06 mm/day-and a standard deviation of 1.24 mm/day.Temperature varied between 22.16 °C and 30.93 °C, with a mean of 27.25 ± 2.24 °C, RH ranged between 91.54% and 58.27%, with an average of 69.95%, and average Rn for this period was 144.32 ± 43.7 W/m 2 .A total of 279.69 mm precipitation was recorded within the same period.During the dry period (DOY 128-153), daily ET ranged between 0.53 and 1.47 mm/day, with an average of 0.99 ± 0.26 mm/day, Rn was between 40.15 and 91.92 W/m 2 , and mean air temperature of 18.68 ± 1.44 °C was recorded; no precipitation was recorded during the same period.
In Elandsberg, in the summer season (DOY 314-346), which is also the dry season, Ta recorded an average of 21 ± 3.45 °C, whereas RH was between 27.76% and 77%, with Rn ranging from 7.83 to 189 W/m 2 ; ET ranged between 1.05 and 4.06 mm/day, with an average of 2.78±0.76mm/day, and a total precipitation of 9.2 mm.The wet, or winter period (DOY 153 and 180) recorded a total of 191 mm rainfall, and daily ET varied between 0.17 and 2.22 mm/day, with a mean of 0.65 mm/day and standard deviation of 0.37 mm/day, Ta ranged from 8.39 to 20.45 °C, RH varied between 8.58% and 100%, whereas Rn recorded daily averages between −14.47 and 46.52 W/m 2 .

ET Models Performance Evaluation
In this subsection, the performance of the evapotranspiration models, TsVI and MOD16-based Penman-Monteith, and MET and GLEAM global products estimates, against in situ measurements is analyzed over two different seasons and ecosystems.This section will look at the performance of the models per location per season, and also compare the different locations per model per season.In Elandsberg, in the summer season (DOY 314-346), which is also the dry season, Ta recorded an average of 21 ± 3.45 • C, whereas RH was between 27.76% and 77%, with Rn ranging from 7.83 to 189 W/m 2 ; ET ranged between 1.05 and 4.06 mm/day, with an average of 2.78±0.76mm/day, and a total precipitation of 9.2 mm.The wet, or winter period (DOY 153 and 180) recorded a total of 191 mm rainfall, and daily ET varied between 0.17 and 2.22 mm/day, with a mean of 0.65 mm/day and standard deviation of 0.37 mm/day, Ta ranged from 8.39 to 20.45 • C, RH varied between 8.58% and 100%, whereas Rn recorded daily averages between −14.47 and 46.52 W/m 2 .

ET Models Performance Evaluation
In this subsection, the performance of the evapotranspiration models, TsVI and MOD16-based Penman-Monteith, and MET and GLEAM global products estimates, against in situ measurements is analyzed over two different seasons and ecosystems.This section will look at the performance of the models per location per season, and also compare the different locations per model per season.The results of the analysis are illustrated in Figures 5 and 6, and the statistical analyses are summarized in Table 4.We also evaluated the modeled Rn and G against values that were measured at the sites.3 highlighting the statistics of the models comparisons.Generally, all the models underestimated ET on both seasons, as shown by the negative biases (Table 3) that ranged between −2.66 mm/day (MET) and −0.79 mm/day (PM) in the wet season, and −0.64 mm/day (GLEAM) and −0.01 mm/day (TsVI) during the dry season.The high underprediction of ET by MET and GLEAM during the wet period is also illustrated in Figures 5a and 6a, whereas PM showed the least underprediction of ET.During the dry season, TsVI (−0.01 mm/day) had the lowest underestimation, and GLEAM (−0.64 mm/day) had the highest underestimation (Figures 5b and 6b).The slopes ranged between 0.19 (PM) and 0.66 (MET), and intercepts were between 1.23 mm/day (MET) and 3.40 mm/day (PM) in the wet season.In the dry season, the slopes ranged from 0.20 (MET) to 3.35 (GLEAM), and the intercepts were between 0.02 mm/day (GLEAM) and 0.66 mm/day (MET).The correlations (R 2 ) of the modeled ET against the measured ET were relatively low, ranging between 0.05 (PM) and 0.45 (MET) during the wet season, and 0.07 (MET) and 0.42 (GLEAM) during the dry season.Although it had the lowest correlation, PM was the best performing model during the wet period, recording the lowest MAE, RMSE and rRMSE of 0.89 mm/day, 1.25 mm/day and 27.34%, respectively, followed by TsVI and GLEAM, which had comparable accuracies of 1.40 and 1.46 mm/day, 1.66 and 1.64 mm/day, and 36.46% and 40.44%,MAE, RMSE and rRMSE, respectively.MET performed the least with MAE, RMSE and rRMSE of 2.66 mm/day, 2.85 mm/day and 67.5%, respectively.In the dry period, TsVI and MET had comparable accuracies of 0.23 and 0.25 mm/day (MAE), 0.29 and 0.32 mm/day (RMSE), 28.34% and 31.09%(rRMSE); GLEAM performed the least with MAE 0.64 mm/day, RMSE 0.67 mm/day, and rRMSE 65.96%.

Elandsberg
Figure 5c,d presents the temporal variation of the LAS-derived ET and the modeled ET estimates, and Figure 6c,d shows the correlations between the measured and estimated ET.In the dry season, there is a general underprediction of ET by the models, as shown in Table 3 and Figure 5c,d and Figure 6c,d, whereas, during the rainy season, GLEAM underestimated ET.The slopes were between 0.012 (GLEAM) and 1.29 (TsVI), and intercepts ranged between −1.51 (MET) and 2.13 (GLEAM) in the dry season; during the wet season the slopes ranged from −0.27 (MET) and 0.57 (MET), and the intercepts were between 0.41 (GLEAM) and 1.41 (PM).The correlations of the modeled ET against the ground ET ranged between 0 (GLEAM) and 0.42 (TsVI) during the dry season, and lowly 0.01 (GLEAM) and 0.12 (MET) during the wet season (Table 3).Although it has higher absolute bias and MAE than MET (−0.09 and 0.65 mm/day), PM (−0.51 and 0.69 mm/day) was best performing with RMSE of 0.85 mm/day, and rRMSE of 28.73%, while MET had RMSE 0.96 mm/day and rRMSE 35

Net Radiation Estimations
Rn and G were computed using the same formulae described above for both ET models at Skukuza and Elandsberg.These were estimated using both meteorological and remote sensing data inputs.In Skukuza, the results show that in the wet season the modeled Rn was more the measured values.During the wet season Rn had exhibited higher accuracies, as shown by R 2 of 0.46, rRMSE of 37.07%, whereas in the dry season the R 2 was extremely low, and rRMSE 22%.On the other hand, the estimation of G was comparably poor at both times.The estimation of Rn and G in Elandsberg was also characterized by low coefficient of determination at both periods, with the wet period giving worse results.To improve the modeling ET using remote sensing based ET models, it is critical to first ensure that intermediate input parameters like Rn and G are accurately estimated.

Discussion
Figure 4 shows high evapotranspiration for both sites during the summer season, even though Elandsberg is located in a winter rainfall region and receives higher rainfall in winter.With higher water availability in Elandsberg in winter, evapotranspiration is expected to be higher, however, low ET is recorded (Figure 4d).In this instance, water availability is not a limiting factor in the process of evapotranspiration in Elandsberg, solar radiation is.Although different times were used for the two sites, this simply illustrates the general climatic trends in the two regions, which experience precipitation at different seasons.
Daily estimates of ET from two models and two products, TsVI, PM, MET and GLEAM, of varying structural complexities, assumptions and parameterizations were compared with ground measurements at two different ecoclimatic regions and times.These models gave differing results, mainly due to algorithm structural errors and input uncertainties and different model sensitivity to the different inputs.Our study shows that no model clearly outperformed the others on both sites and times.In Skukuza, PM gave the least error ET estimates in the wet season, whereas during the dry period TsVI and MET gave relatively comparable results, although TsVI estimates were slightly better.In Elandsberg also, no model distinctly outperformed the others during the dry season, as evidenced by the interchange in hierarchy of the statistics between PM and MET; during the wet season, GLEAM had the least bias, as shown by low MAE and RMSE, despite having the lowest R 2 of 0.01 (Table 3).At both sites, we are looking at periods of high radiation and ET during the summer season, despite limited rainfall in Elandsberg during this period.The poor performance of the models during the wet season in Elandsberg could be related to the weather conditions, particularly since it was raining during these days (Figure 4).The low correlations that were recorded in our study, notably for Elandsberg wet season, are comparable to other studies that have been done in dry ecosystems [30,64].Vinukollu et al. [32] reported low correlations in grasslands and woody savannas with all the models they tested in their study.They reported Kendall's tau of 0.01, 0.32 and 0.33 in woody savanna Tonzi; 0.51, 0.27 and 0.37 in grassland Audubon; and 0.55, 0.54, and 0.59 in grassland Fort Peck, for SEBS, PM-Mu and PT-Fi models, respectively.In their extensive evaluation of SEBS, PT-JPL, PM-Mu and GLEAM models, across different ecoclimatic regions, McCabe et al. [65] also found that the PT-based models, PT-JPL and GLEAM, performed better than the others in most ecoclimatic regions.Low correlations and accuracies of the models were also recorded as the aridity increased, such that they excluded the results of shrublands in their further analysis.Ershadi et al. [30] evaluated the performance of PM, advection-aridity, SEBS and PT-JPL models, and showed that there is no model that consistently performed well across all biomes.Hu et al. [66] also reported low performance of the MOD16 ET and MET on shrublands located in semi-arid climates.
Model performance is affected by different attributes, including structural complexities, model assumptions, parameterizations, and amount of data required.To estimate ET using these models, both meteorological data measured from the study sites and remote sensing data were used as inputs.The error and uncertainties of the data inputs are propagated into the models' outputs.For instance, the remote sensing inputs of LST and surface emissivity are instantaneous in the case, and eight-day averages in the cases of vegetation indices, whereas the ground measurements of air temperature, humidity, atmospheric pressure, wind speed are daily averages.In addition, MODIS products that were used as inputs in this study come with their level of uncertainty and/or error, presented as part of the datasets as reported in their user guides [67][68][69].The MODIS products used as inputs in this study are affected by artifacts caused by clouds, atmospheric aerosols, instrument errors and uncertainties of retrieval algorithms.There is, thus, a need for these data inputs to be improved to improve the ET estimation.The spatial scale difference between MODIS image footprint of 1 km and the flux tower and LAS footprints, confounded by the heterogeneity of the landscape, the wind velocity and direction, and atmospheric stability, within the satellite footprint also contribute error and uncertainty [26,35].As shown by the flux foot print model, and discussed in Ramoelo et al. [35], with a measurement height of 16 m, the Skukuza flux tower footprint is approximately 1.6 km, i.e., based on the general rule of the thumb which suggests the fetch:measurement height ratio of 100,484-m:1 m [70,71].This indicates that the flux tower provides a good match to the MODIS pixel size.Although flux footprint modeling gives a good estimation of the spatial discrepancy between flux tower measurements and the image pixel, it is beyond the scope of this study.The coarser the image resolution, the higher the landscape heterogeneity and spatial mismatch, which introduces errors, as shown by McCabe and Wood [72] who showed that while Landsat and ASTER higher accuracies in ET estimates, MODIS was less accurate.For evapotranspiration models, it is important to have an accurate estimation of Rn which is a critical component of each of these models.A low accuracy in Rn, as shown in our study also contributed to the error and inaccuracies of the ET estimates.
GLEAM and TsVI are PT based, which is a simplified version of the PM, with the TsVI triangle method being a relatively simple formulation that employs empirical means (using to normalize the PT parameter (ϕ) using the LST-vegetation indices triangle feature space, to estimate the evaporative fraction; it requires fewer inputs and parameterizations required (see Table 2) [73,74].The simplification of the model in the estimation of EF also means that the computing complexities of surface and aerodynamic resistances are avoided, thus reducing errors and uncertainty in ET estimation [2].The main challenge of the TsVI method is the subjective determination of the wet and dry edges from the triangle feature space and the neglect of local advection in its formulation.In addition, during the rainy season or in areas of low variability in vegetation cover range, the triangle feature space is hard to establish, as evidenced in Elandsberg, where it was generally outperformed by other models.In other instances, the heterogeneity of the land surface, together with atmospheric forcing complicate the establishment of the TsVI relationship.Since the TsVI feature space is established empirically, the method is site specific.During the rainy days, like in the Elandsberg winter period and Skukuza summer period, the determination of the wet and dry edges presents a challenge, hence the non-performance of the model during this period.GLEAM, on the other hand, is more comprehensive, combining the PT equation, a soil moisture stress computation, and a Gash analytical model to compute ET as a total of transpiration from tall canopy, short vegetation, soil evaporation and canopy interception loss [26].A plus for the GLEAM model is that in water limited regions, atmospheric water demand is constrained by precipitation, surface soil moisture and vegetation optical depth (VOD, which is a proxy for leaf water content [75]).Each of the model components within the GLEAM structure has its own assumptions and complexities with varying levels of error and uncertainty, which are propagated to the final ET estimate.GLEAM, like the PM used in this study also computes interception loss separately, as shown in Miralles et al. [26] and Mu et al. [21].Miralles et al. [26] assessed the GLEAM product across different biomes and showed daily average correlation (R) of 0.83.They, however, stated that the distinct seasonal cycle of evaporation in the dry regions probably balanced out the correlation coefficients positively.In addition to the known challenges of the PT equation, the main challenge in validating GLEAM ET using flux tower derived ET is the coarse spatial resolution of the global product.The MET product is derived from a scheme that is different from the other estimates, based on the energy balance budget.The daily MET product is an aggregation of 30-min ET values obtained.This presents a challenge when, for instance, there are gaps in these 30-min ET data, resulting in underestimation.In addition, the issue of spatial heterogeneity plays a part in the inaccuracies of this course resolution ET product.
PM is the most robust of these models, and theoretically it should present the best performances compared with the other models.The results (Table 3) however, show that this version of PM was only good during periods of high ET.One of the biggest challenges of the PM based models is the parameterization of the aerodynamic and surface resistances, including upscaling stomatal to canopy resistance.LAI is an important input in the parameterization of canopy resistance, as it is linked with the biophysical control of vegetation on ET [20,21,74].Hu et al. [66] explored the relationship between PM-MOD16 estimated ET and LAI, showing that the two are closely related, especially for the savannas and the deciduous broadleaf forests.Because soil water availability plays a key role in ET in semi-arid regions, it is important that models include a soil water constraint function in the model.However, in this version of PM, relative humidity and VPD were used as proxy for soil water in the estimation of soil evaporation, hence low accuracy of the model during low ET periods.Currently, research is focused on incorporating soil water constraint function in the ET modeling, especially in dry regions.Di et al. [76] incorporated two layers of relative soil moisture parameters in the PM model to parameterize the surface resistance, and added a multiplier in the vegetation surface resistance model to cater for the influence of the relative soil moisture in the root zone.Sun et al. [57] also investigated the incorporation of soil moisture in the PM-Mu method to constrain soil evaporation by using actual soil water content and the soil water content at saturation to compute the soil resistance, which they also substituted with a soil moisture index derived from the TsVI triangle method.These studies showed improved estimations of ET by the PM algorithm during water constrained periods.For the estimation of ET to cover different ecosystems, PM would have to be assessed further and modified, especially for semi-arid climatic regions.

Conclusions
Accurate estimates of ET are essential especially in semi-arid and arid regions where there is less water being competed for by different users.Different remote sensing based models and products, of varying complexities and data input requirements, are available, and their applicability at varying ecoclimatic regions and scales is consistently under scrutiny.
This study, thus, presented an evaluation of four ET models and global products, i.e., TsVI, PM and MET and GLEAM global products, in two semi-arid eco-climates.Our results show that there is essentially no model that clearly outperforms others at the two sites.Low coefficients ranging between 0 and 0.45 were recorded on both sites, during both seasons.It was also observed that, during periods of high ET at both sites, i.e., Skukuza in the wet season and Elandsberg in the dry season, PM was relatively more accurate than the other models and products, with rRMSEs of lower than 30%.In Skukuza, TsVI marginally outperformed MET during the dry period, whereas GLEAM gave the least accurate estimates.ET estimation during the rainy season in Elandsberg was quite poor with rRMSEs of over 70%, with GLEAM being most accurate.
The conclusion, therefore, is that none of the models performed well, as shown by low R 2 and high errors in all the models.PM gave the least errors during periods of high ET on both sites, whereas modeling low ET was a challenge.
These results present a prerequisite on the next stage of our study, in which we will investigate the error and uncertainty in the PM model simulations, inclusive input data (both remote sensing and meteorological input data), model structure, and parameterizations.In doing this, Ershadi et al. [30] state that it would allow for the model diagnosis and identification of the main sources of error in ET estimation, in this case in water-scarce regions.
The area is characterized by a catenal pattern of soils and vegetation, with broad-leaved Combretum savanna on the crests dominated by Combretum apiculatum, and fine-leaved Acacia savanna in the valleys dominated by Acacia nigrescens.The vegetation is mainly open woodland, with approximately 30% tree canopy cover of mixed Acacia and Combretum savanna types.Tree canopy height is 5-8 m with occasional trees (mostly Sclerocarya birrea) reaching 10 m.The grassy and herbaceous understory comprises grasses such as Panicum maximum, Digitaria eriantha, Eragrostis rigidor, and Pogonarthria squarrosa.The flux tower was placed on a vegetation transition to measure fluxes from the different types [40].Remote Sens. 2017, 9, 307 4 of 20 squarrosa.The flux tower was placed on a vegetation transition to measure fluxes from the different types [40].

Figure 2 .
Figure 2. Skukuza eddy covariance flux tower station in Kruger National Park.

Figure 2 .
Figure 2. Skukuza eddy covariance flux tower station in Kruger National Park.Figure 2. Skukuza eddy covariance flux tower station in Kruger National Park.

Figure 2 .
Figure 2. Skukuza eddy covariance flux tower station in Kruger National Park.Figure 2. Skukuza eddy covariance flux tower station in Kruger National Park.
left) was located at 33.47404 • S; 19.06239 • E, while the receiver was placed at approximately 900 m from the transmitter (Figure 3 right), at 33.47001 • S; 19.05526 • E.

Figure 3 .
Figure 3. (left) LAS transmitter of infrared radiation; and (right) receiver of the optically modified infrared beam located about 900 m away from the transmitter.

Figure 3 .
Figure 3. (left) LAS transmitter of infrared radiation; and (right) receiver of the optically modified infrared beam located about 900 m away from the transmitter.

20 Figure 5 .
Figure 5.Time series of measured and modeled ET for the: wet (a); and dry (b) periods in Skukuza eddy covariance flux tower site; and the dry (c); and wet (d) periods in Elandsberg LAS site.

Figure 5 .
Figure 5.Time series of measured and modeled ET for the: wet (a); and dry (b) periods in Skukuza eddy covariance flux tower site; and the dry (c); and wet (d) periods in Elandsberg LAS site.

3. 2
Figure 5a,b illustrates the temporal variation of the flux tower based and the modeled ET, Figure 6a,b show the correlations between the flux tower and modeled ET, and Table3highlighting the statistics of the models comparisons.Generally, all the models underestimated ET on both seasons, as shown by the negative biases (Table3) that ranged between −2.66 mm/day (MET) and −0.79 mm/day (PM) in the wet season, and −0.64 mm/day (GLEAM) and −0.01 mm/day (TsVI) during the dry season.The high underprediction of ET by MET and GLEAM during the wet period is also illustrated in Figures5a and 6a, whereas PM showed the least underprediction of ET.During the dry season, TsVI (−0.01 mm/day) had the lowest underestimation, and GLEAM (−0.64 mm/day) had the highest underestimation (Figures5b and 6b).The slopes ranged between 0.19 (PM) and 0.66 (MET), and intercepts were between 1.23 mm/day (MET) and 3.40 mm/day (PM) in the wet season.In the dry season, the slopes ranged from 0.20 (MET) to 3.35 (GLEAM), and the intercepts were between 0.02 mm/day (GLEAM) and 0.66 mm/day (MET).The correlations (R 2 ) of the modeled ET against the measured ET were relatively low, ranging between 0.05 (PM) and 0.45 (MET) during the wet season, and 0.07 (MET) and 0.42 (GLEAM) during the dry season.Although it had the lowest correlation, PM was the best performing model during the wet period, recording the lowest MAE, RMSE and rRMSE of 0.89 mm/day, 1.25 mm/day and 27.34%, respectively, followed by TsVI and GLEAM, which had comparable accuracies of 1.40 and 1.46 mm/day, 1.66 and 1.64 mm/day, and 36.46% and 40.44%,MAE, RMSE and rRMSE, respectively.MET performed the least with MAE, RMSE and rRMSE of 2.66 mm/day, 2.85 mm/day and 67.5%, respectively.In the dry period, TsVI and MET had comparable accuracies of 0.23 and 0.25 mm/day (MAE), 0.29 and 0.32 mm/day (RMSE), 28.34% and 31.09%(rRMSE); GLEAM performed the least with MAE 0.64 mm/day, RMSE 0.67 mm/day, and rRMSE 65.96%.

Table 1 .
Measurements taken at Skukuza eddy covariance system and the instruments used.

Table 2 .
A summary of the different ET algorithms and where they have been applied.

Table 3 .
Summary of meteorological input variables and their measurement instruments for Skukuza and Elandsberg.

Table 4 .
Statistics of estimated daily ET against measured ET for Skukuza flux and Elandsberg LAS.

Table 4 .
Statistics of estimated daily ET against measured ET for Skukuza flux and Elandsberg LAS.
.03%, which was close to TsVI, having recorded RMSE 1.05 mm/day and rRMSE 35.42%.The least accurate in this instance was GLEAM, which also had comparable results of bias −0.65 mm/day, MAE 0.94 mm/day, RMSE 1.15 mm/day, and rRMSE 40.73%.In the wet season, the performance of all the models were very poor, with all of them recording rRMSE of over 70%.The best performing model, GLEAM had rRMSE 73.02%, RMSE 0.42 mm/day, MAE 0.28 mm/day and bias −0.19 mm/day.Although the rRMSE for TsVI and MET were 88.91% and 119.55%, respectively, they generally performed very close to each other, with MAE and RMSE of 0.46 and 0.52 mm/day, 0.60 and 0.63 mm/day, respectively.The lowest performing model, PM, had a bias, MAE, RMSE and rRMSE of 0.61, 0.69, 0.77 mm/day, and 114.64%, respectively.