Uncertainty and Sensitivity Analysis of a Remote-Sensing-Based Penman–Monteith Model to Meteorological and Land Surface Input Variables

: This study analysed the uncertainty and sensitivity of core and intermediate input variables of a remote-sensing-data-based Penman–Monteith (PM-Mu) evapotranspiration (ET) model. We derived absolute and relative uncertainties of core measured meteorological and remote-sensing-based atmospheric and land surface input variables and parameters of the PM-Mu model. Uncertainties of important intermediate data components (i.e., net radiation and aerodynamic and surface resistances) were also assessed. To estimate the instrument measurement uncertainties of the in situ meteorological input variables, we used the reported accuracies of the manufacturers. Observational accuracies of the remote sensing input variables (land surface temperature ( LST ), land surface emissivity ( ε s ), leaf area index ( LAI ), land surface albedo ( α )) were derived from peer-reviewed satellite sensor validation reports to compute their uncertainties. The input uncertainties were propagated to the ﬁnal model’s evapotranspiration estimation uncertainty. Our analysis indicated relatively high uncertainties associated with relative humidity ( RH ), and hence all the intermediate variables associated with RH , like vapour pressure deﬁcit ( VPD ) and the surface and aerodynamic resistances. This is in contrast to other studies, which reported LAI uncertainty as the most inﬂuential. The semi-arid conditions and seasonality of the regional South African climate and high temporal frequency of the variations in VPD , air and land surface temperatures could explain the uncertainties observed in this study. The results also showed the ET algorithm to be most sensitive to the air-land surface temperature difference. An accurate assessment of those in situ and remotely sensed variables is required to achieve reliable evapotranspiration model estimates in these generally dry regions and climates. A signiﬁcant advantage of the remote-sensing-based ET method remains its full area coverage in contrast to classic-point (station)-based ET estimates. ﬁnal A B ﬁnal PM-Mu


Introduction
Evapotranspiration is dependent on meteorological variables such as air temperature (Tair), solar radiation (Rs), humidity (RH) and wind speed (u) and biophysical characteristics of the land surface and vegetation. It is considered the most uncertain component of the hydrological cycle due to its variation both in space and time, and the complex hydrometeorological processes involved. Hence, it is a challenging process to measure of how the model reacts to any change in climatic variables, such as air temperature, net radiation and water availability. Furthermore, through this section, an analysis of how land use/land cover change impacts on ET variation is provided.
The emergence of RS-based ET modelling presents another opportunity to assess not only the reaction of the ET models to climatic input variables but also to land surface parameters, which culminates in how land-use change impacts on ET. Few studies have investigated the sensitivity of ET models to remote sensing input parameters [29,30]. This indicates that more work needs to be done to understand the sensitivity of ET models to remote sensing land surface parameters, especially considering the number of remote sensing-based ET models that have been extensively evaluated across different bioclimates.
A previous study by Majozi, et al. [31] reported a distinct reasonable performance of the PM model as modified by Mu, et al. [32], Mu, et al. [33] (PM-Mu) during high ET periods. Consequently, this study investigated the sensitivity of this model to its input variables, i.e., both meteorological and land surface characteristics, in situ and remote sensing. This study also quantified the uncertainty of the input variables and how these were propagated into the final ET uncertainty. Generally, the PM model is a structurally complex and data-intensive model presenting a combination of the energy balance and aerodynamic components. This makes it important to endeavour to understand and quantify potential errors and uncertainties of the input data and how these impact the final ET output. This process will, therefore, identify the inputs that are most influential and correlated with the dependent variable in a semi-arid environment, in order to improve parameterisations that could eventually improve our results. This study will also highlight the possible implications of land use/land cover and climate change on the hydrology of the study areas.

Site Description
ET was estimated using data from two eddy covariance flux tower sites located in savanna and grassland southern African ecoregions ( Figure 1). The year 2012 daily data were selected for both sites and considered in this study: i.
The Skukuza FLUXNET site, located in a savanna ecosystem in the Kruger National Park, South Africa, sits at 365 m above mean sea level. The site is characterised by low rainfall averaging 550 ± 160 mm per annum between November and April, with notable inter-annual variability, and temperatures ranging between 15.6 and 29.6 • C, with a mean of 22.6 • C. Soils in this part of the park are generally shallow, comprising coarse sandy to sandy-loam texture. The vegetation is mainly open woodland, with approximately 30% tree canopy cover of mixed Acacia and Combretum savanna types. The canopy height is 5-8 m, with occasional trees (mostly Sclerocarya birrea) reaching 10 m. The grassy and herbaceous understorey comprises grasses such as Panicum maximum, Digitaria eriantha, Eragrostis rigidor and Pogonarthria squarrosa. The eddy covariance system, which has been running since 2000, was installed on a vegetation transition characterised 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 [34,35]. ii.
Welgegund flux tower site (26 • 34 10"S, 26 • 56 21"E) is located on a semi-arid, subtropical grazed grassland plain. It is situated approximately 100 km west of Johannesburg, in South Africa. The mean annual rainfall is 540 ± 112 mm, spreading between October and April. Temperature ranges between 0 and 30 • C, with an average of 18

Remote-Sensing-Based Penman-Monteith Model
The Penman-Monteith model as modified by Mu et al. [33,34] was assessed in this study. The latent heat flux was estimated as: where λE is the latent heat flux (Wm −2 ) and λ is the heat of vaporisation (Jkg −1 ), s is the slope of the curve relating saturated vapour pressure to temperature (PaK −1 ), A is available energy (Wm −2 ), is the air density (kgm −3 ), Cp is the specific heat capacity of air (Jkg −1 K −1 ), γ is the psychrometric constant (PaK −1 ), esat is the saturation vapour pressure (kPa), e is the actual vapour pressure (kPa), where esat -e = vapour pressure deficit (VPD, kPa), ra (sm −1 ) is the aerodynamic resistance and rs (sm −1 ) is the canopy resistance, which is the reciprocal of canopy conductance gc (gc = 1rc −1 ). On top of estimating ET as a sum of evaporation from moist soil, interception and transpiration, the daytime and nighttime ET were further computed separately [33,34]. Instead of using NDVI to compute the fraction of vegetation cover (Fc), they used the Fraction of Absorbed Photosynthetically Active Radiation (FPAR) as a surrogate of vegetation cover fraction, with another modification to the derivation of soil heat flux. Moreover, the model was modified by separating (i) dry and wet (interception) canopy surfaces, (ii) soil surface into saturated and moist surface, and (iii) improving stomatal conductance, aerodynamic resistance and boundary layer resistance estimates.

Remote-Sensing-Based Penman-Monteith Model
The Penman-Monteith model as modified by Mu et al. [32,33] was assessed in this study. The latent heat flux was estimated as: where λE is the latent heat flux (Wm −2 ) and λ is the heat of vaporisation (Jkg −1 ), s is the slope of the curve relating saturated vapour pressure to temperature (PaK −1 ), A is available energy (Wm −2 ), is the air density (kgm −3 ), Cp is the specific heat capacity of air (Jkg −1 K −1 ), γ is the psychrometric constant (PaK −1 ), e sat is the saturation vapour pressure (kPa), e is the actual vapour pressure (kPa), where e sat -e = vapour pressure deficit (VPD, kPa), r a (sm −1 ) is the aerodynamic resistance and r s (sm −1 ) is the canopy resistance, which is the reciprocal of canopy conductance gc (gc = 1rc −1 ). On top of estimating ET as a sum of evaporation from moist soil, interception and transpiration, the daytime and nighttime ET were further computed separately [32,33]. Instead of using NDVI to compute the fraction of vegetation cover (Fc), they used the Fraction of Absorbed Photosynthetically Active Radiation (FPAR) as a surrogate of vegetation cover fraction, with another modification to the derivation of soil heat flux. Moreover, the model was modified by separating (i) dry and wet (interception) canopy surfaces, (ii) soil surface into saturated and moist surface, and (iii) improving stomatal conductance, aerodynamic resistance and boundary layer resistance estimates.
The core input variables used in this model are Tair, RH, land surface temperature (LST), surface emissivity (ε s ), leaf area index (LAI) and land surface albedo (α), and were used to derive intermediate inputs like net solar radiation (Rnet), vapour pressure deficit Remote Sens. 2021, 13, 882 5 of 18 (VPD), the slope of the saturated vapour -air temperature curve (∆), the air and saturated air vapour pressures (e, e sat ) and the aerodynamic and surface resistances (r a , r s ).

Uncertainty and Sensitivity Analysis
Uncertainty and sensitivity analyses were performed on the PM-Mu to quantify input uncertainties and how these are propagated to the final ET uncertainties, and identify the inputs and parameters that are most important in modelling ET in a semi-arid environment. We assessed the uncertainty of each input variable, i.e., the direct measured variable uncertainty (core input), derived input variable uncertainty (intermediate variable) and remote-sensing-based input variable uncertainty using Type A and Type B uncertainty methods (see Section 2.3.1). The total uncertainty on the model simulations, i.e., model output uncertainty, was then evaluated by uncertainty propagation using the Gaussian uncertainty analysis method.
Based on the PM-Mu model, ET is defined as a function (f ) of meteorological point measurements of Tair and RH, and spatially explicit remote sensing estimates of LST, and land surface characteristics such as LAI, fraction of green vegetation cover (Fc)/fraction of photosynthetically active radiation (FPAR), NDVI/EVI and surface emissivity (ε s ), which are biome/or land cover characteristics defining parameters. Table 1 below summarises the core and intermediate inputs and parameters used in this study.
The generic model can be presented as: where x 1 to x n represents the n input variables and parameters of the PM-Mu model. The change in ET, i.e., ∆ET, resulting from errors and/or uncertainties in input variables (∆x i ), is then expressed as: The study analysed the uncertainty and sensitivity of the PM-Mu input variables and outputs aimed at: i.
Estimating the uncertainty of model inputs and parameters, i.e., the meteorological and land surface characteristics, representing both point-and remote-sensingbased inputs; ii.
Propagating input uncertainties through to the ET model and computing output uncertainties; iii.
Estimating the sensitivity coefficients of the model inputs. The core inputs of the PM-Mu equation include relative humidity (RH) and air temperature (Tair). We estimated the uncertainty of each core input variable as a combination of Type A and Type B uncertainties. In this study, where each measured input variable was a daily mean of 30-min recordings, we computed the Type A (U A ) standard uncertainty of the meteorological inputs as the standard deviation of the daily mean: where x i is the input value of the variable or parameter under consideration, x i is the average value of the measured values calculated from the n number of independent observations and v is the degrees of freedom equal to n-1. Type B standard uncertainty (U B ) was also computed for the meteorological input variables, based on the instrument manufacturer's published accuracies. The quoted accuracies of the measurement instruments are summarised in Table 2. They were estimated using Equation (5): where a is the quoted accuracy specification from the manufacturer and includes calibration information from calibration certificates. For meteorological data inputs, the combined standard uncertainty was then estimated as: The combined uncertainty was then converted to relative uncertainty for detailed comparison and analysis as follows:

Remote-Sensing-Based Input Uncertainties
In this study, we used the Moderate Resolution Imaging Spectroradiometer (MODIS) Terra/Aqua products as data inputs. The following MODIS products were used as inputs for this study: MOD11A1 LST and εs, MOD15A2 LAI and MCD43A3 α. These datasets were downloaded from the NASA Land Processes Distributed Active Archive Centre (LP DAAC) website using the USGS EarthExplorer platform (https://earthexplorer.usgs.gov/). These downloaded MODIS products were loaded onto ArcGIS together with the two flux tower GPS points. The values for the loaded points were then extracted from the MODIS products using the point extraction feature on ArcGIS software. These remotesensing-derived inputs have uncertainties due to several factors, including model algorithm Remote Sens. 2021, 13, 882 7 of 18 structure and their input variables. The uncertainties of the remotely sensed input variables used in this study (LST, εs, LAI, α) were extracted mainly from their published Algorithm Theoretical Basis Documents (ATBD). Based on the quoted errors (Table 3), the uncertainties were then estimated using Equation (5). We give a short description of each variable and associated remote sensor below.

•
Land surface temperature (LST) and surface emissivity (ε s ): these variables are essential in land surface-atmosphere studies, including the estimation of evapotranspiration and atmospheric water vapour. In our study, we used the MODIS-derived MOD11A1 V006 product, which is generated from the thermal infrared channels 31 (10.78 to 11.28 µm) and 32 (11.77 to 12.27 µm) using the physically-based day-night split-window algorithm by [37]. The uncertainties associated with these products are extensively discussed in the MODIS Land-Surface Temperature ATBD [38][39][40]. They indicate an absolute error of 1 K for LST, which can increase up to 5 K in arid regions. For surface emissivity, the absolute accuracy is reported to be 0.02.

•
Land surface albedo (α): defined as a dimensionless characteristic of the soil-plant canopy system representing the fraction of total solar energy reflected by the surface, it is expressed as the ratio of the radiant energy scattered upward by a surface in all directions to that received from all directions, integrated over the wavelengths of the solar spectrum. Surface albedo is one of the key geophysical parameters that control the surface energy budget. The MODIS bi-directional reflectance distribution function (BRDF) and albedo product (MCD43A3 version VOO6) was used in this study. This product was derived using a kernel-driven semi-empirical BRDF model using the RossThick-LiSparse kernel functions for characterising isotropic, volume and surface scattering [41][42][43]. Studies have given an absolute accuracy of 0.02 to 0.05 as a requirement for climate modelling [44,45], with other validation studies [46,47] reporting errors falling within the 0.02 accuracy.

•
Leaf Area Index (LAI): defined as the total one-sided green leaf area per unit ground surface area, it is also dimensionless. This variable measures the total amount of leaf material in an ecosystem. It is used in the estimation of biogeochemical processes like photosynthesis, evapotranspiration and net primary production. The MOD15A2 V005 product used in this study was derived using the three-dimensional radiative transfer (3D RT) model [48,49]. The product ATBD reports the accuracy of the LAI product at 0.2 [48]. Furthermore, a review by Fang, et al. [50] summarises the uncertainties of MODIS, CYCLOPES and GLOBCARBON LAI products under different biomes, showing a relative uncertainty of 0.26 in the savanna biome for the MODIS product.

Intermediate Input Uncertainty
For intermediate inputs that were derived from the core input variables, like net radiation and surface and aerodynamic resistances, the standard uncertainties were estimated as combined standard uncertainties of their inputs. The Gaussian error propagation method, which describes how the uncertainties in the variables x and y propagate into a function f (x, y) and assumes that the model inputs are uncorrelated, was used, as shown in Equation (8): is the partial derivative of y, which is the output variable, with reference to input variables x i to x n , also called sensitivity coefficient.
Each estimated input variable uncertainty was also propagated to the final ET output uncertainty using Equation (8).

Sensitivity Analysis
One of the aims of SA is to identify and rank input variables according to their importance in modelling a particular phenomenon. This is done to identify the input variables that require a more accurate measurement to reduce model output variance to a minimum. The sensitivity of the PM-Mu estimated daily ET was determined by varying one input variable at a time within ±20% ranges. First, ET was computed with the initial input variables, then one variable was perturbed by 5% within ±20% whilst the rest of the inputs were held constant every day for the whole year of 2012 and the new ET values were recorded. Then, the sensitivity coefficient, S, was computed using Equation (9), after which an overall average was calculated: where Y i is the ET recorded when you vary one variable at a time at each percentage step, and Y 0 is the initial ET value.

Results
Uncertainty analysis gives a range of values likely to enclose the true value, and thus the confidence of the modelled values, and includes all possible sources of error. Meanwhile, sensitivity analysis ranks the input variables according to their sensitivity to errors in a model. In our study, we quantified the uncertainty of the PM-Mu ET model input variables at two FLUXNET sites in semi-arid ecosystems, Skukuza and Welgegund, and analysed how these propagated through to the model's final output uncertainty. Figure 2 illustrates the annual variations of Tair, bound with the absolute uncertainties for the two study areas. With a mean annual Tair of 24 ± 3.57 • C, the absolute standard uncertainty was 1.5 ± 0.7 • C for Skukuza. Welgegund's mean annual Tair was 20 ± 5.2 • C, with the absolute standard uncertainty similar with that of Skukuza, with an average of 1.5 ± 0.42 • C. Meanwhile, this translated to Tair's relative uncertainty, between 0.5 and 7.6%, with an average of 3.1 ± 1.5% for Skukuza, and it varied between 0.9 and 10.6% with a mean of 4.0 ± 1.7% for Welgegund.

Core Input Variables Uncertainty
As shown in Figure 3, the annual daily average RH varied between 34 and 94% during the same period for Skukuza, and the mean absolute uncertainty was 5 ± 1.5%. This result was converted to an average relative uncertainty of 10 ± 3.5%. On the other hand, Welgegund's annual daily average RH ranged between 13 and 98%, and the mean absolute uncertainty was 6 ± 1.9%, while the relative uncertainty was 14.2 ± 5.4%.
Relative uncertainties of both Tair and RH have strong seasonal variability, with relative uncertainties being higher during the drier months of the year, i.e., between April and September, compared to the wet months. During this period, daily temperatures tend to be highly variable throughout the day, hence the high Type A standard uncertainty. Furthermore, there was much less variation in Tair's relative uncertainty between the two sites compared to RH's relative uncertainty. Remote Sens. 2021, 13, x FOR PEER REVIEW 9 of 20 As shown in Figure 3, the annual daily average RH varied between 34 and 94% during the same period for Skukuza, and the mean absolute uncertainty was 5 ± 1.5%. This result was converted to an average relative uncertainty of 10 ± 3.5%. On the other hand, Welgegund's annual daily average RH ranged between 13 and 98%, and the mean absolute uncertainty was 6 ± 1.9%, while the relative uncertainty was 14.2 ± 5.4%. Relative uncertainties of both Tair and RH have strong seasonal variability, with relative uncertainties being higher during the drier months of the year, i.e., between April and September, compared to the wet months. During this period, daily temperatures tend to be highly variable throughout the day, hence the high Type A standard uncertainty. Furthermore, there was much less variation in Tair's relative uncertainty between the two sites compared to RH's relative uncertainty.

Intermediate Input Uncertainty
This subsection reports how the core input uncertainties estimated above were propagated onto the ET intermediate inputs, i.e., Rnet and the aerodynamic and surface resistances. As shown in Figure 3, the annual daily average RH varied between 34 and 94% during the same period for Skukuza, and the mean absolute uncertainty was 5 ± 1.5%. This result was converted to an average relative uncertainty of 10 ± 3.5%. On the other hand, Welgegund's annual daily average RH ranged between 13 and 98%, and the mean absolute uncertainty was 6 ± 1.9%, while the relative uncertainty was 14.2 ± 5.4%. Relative uncertainties of both Tair and RH have strong seasonal variability, with relative uncertainties being higher during the drier months of the year, i.e., between April and September, compared to the wet months. During this period, daily temperatures tend to be highly variable throughout the day, hence the high Type A standard uncertainty. Furthermore, there was much less variation in Tair's relative uncertainty between the two sites compared to RH's relative uncertainty.

Intermediate Input Uncertainty
This subsection reports how the core input uncertainties estimated above were propagated onto the ET intermediate inputs, i.e., Rnet and the aerodynamic and surface resistances.

Intermediate Input Uncertainty
This subsection reports how the core input uncertainties estimated above were propagated onto the ET intermediate inputs, i.e., Rnet and the aerodynamic and surface resistances.

Net Radiation Uncertainty
Net radiation (Rnet) estimation depends on a number of atmospheric and land surface variables, including α, ε s , ε a , LST and Tair. The overall relative uncertainty of Rnet was 4.0 ± 0.6% of the estimated 558.0 ± 105.2 Wm −2 daytime Rnet in Skukuza, whereas for Welgegund, a 2.8 ± 0.8% relative uncertainty was reported for the derived 556.4 ± 127.7 Wm −2 Rnet. As summarised in Table 4, for Skukuza, a mean relative air temperature (Tair) uncertainty of 3.1% was associated with relative Rnet uncertainties of 23.28 ± 9.85% of the total Rnet uncertainty, whereas a land surface temperature (LST) error of 3.5 K contributed 59.31 ± 12.87% to the total uncertainty, which was the highest contributor to Rnet uncertainty. The surface emissivity (ε s ) error of 0.02 contributed an Rnet uncertainty between 4 and 7 Wm −2 , while the albedo (α) uncertainty contributed to an average 2.8 ± 0.42 Wm −2 . Similar results were found for Welgegund, where a mean relative Tair uncertainty of 4% resulted in an Rnet relative uncertainty of 30.25 ± 10.23%, whereas 89.42 ± 22.07% of the total Rnet uncertainty was attributed to an LST error of 3.5 K, showing that it was the highest source of uncertainty in Rnet estimation. ε s and α contributed mean relative uncertainties of 38.81 ± 9.43 and 22.06 ± 10.31%, respectively.

Aerodynamic and Surface Resistances
In the estimation of wet canopy evaporation, the aerodynamic resistance to evaporated water on the wet canopy surface (r a wc ) is a function of Tair, LAI and RH (in the form of wet surface fraction (Fwet)), whereas the surface resistance to evaporated water on the wet canopy (r s wc ) is a function of LAI and RH. Further, in plant transpiration estimation, the aerodynamic resistance to water vapour from a dry canopy surface (r a t ) is a function of Tair only; and the canopy resistance to transpired water (r s t ) is estimated from LST, LAI, Tmin and RH. In the computation of soil evaporation, both the surface (r s tot ) and aerodynamic resistances (r a s ) to water vapour from the soil surface are a function of LST and VPD (which is indirectly RH).
Our results, as illustrated in Table 5 (only the standard uncertainties for resistances are shown here), show that the mean standard uncertainties for r a wc were 0.0011 ms −1 ± 6.25% and 0.001 sm −1 ± 0.17% for Skukuza and Welgegund, respectively. Of the total standard r a wc uncertainty, Tair contributed the highest uncertainty of average 91 ± 5.0% and 96.02 ± 16.54%, with low contributions from the LAI and RH uncertainties for Skukuza and Welgegund, respectively. Meanwhile, r s wc standard uncertainties were an average 10.34 ± 10.0 sm −1 and 18.02 ± 19.0 sm −1 , respectively, with RH uncertainty contributing most to the total r s wc uncertainty (approximately 80%), on both sites. In the estimation of r a t standard uncertainty, the values ranged between 0.00019 and 0.0038 sm −1 , and 0.00031 and 0.0032 sm −1 (low average relative uncertainties of 0.81 ± 0.36% and 0.84 ± 0.22%) for Skukuza and Welgegund, respectively. These low values indicate that Tair uncertainties have the least effect in the estimation of r a t uncertainty. Total standard r s t uncertainty ranged from 0 to 90 sm −1 (mean relative uncertainty of 8.82 ± 2.71%) for Skukuza; whereas for Welgegund, it was between 20 and 146 sm −1 , (average relative uncertainty of 8.4 ± 1.7%). The r a s relative uncertainties were on average around 2% for both sites. Finally, r s tot standard uncertainty ranged from 0.45 to 0.73 sm −1 , and 0.35 to 0.77 sm −1 , for Skukuza and Welgegund, respectively, an average 1% relative uncertainty for both sites. Of the total uncertainty, LST uncertainty contributed the most of the two input variables with an average of 58% and 63%, whereas 5.25% of the total r s tot uncertainty was attributed to VPD uncertainty, for Skukuza and Welgegund, respectively.

Uncertainty in Evapotranspiration
The final estimate of ET uncertainty is a result of uncertainties propagated from the measured and remote sensing input variables, through intermediate parameters, up to the final ET uncertainty. The standard uncertainty was computed for each ET component, i.e., evaporation from intercepted rainfall (wet canopy), transpiration and soil evaporation, and ultimately combined to give the total ET uncertainty.
In Skukuza (Table 6), of the 0.038 mmday −1 wet canopy evaporation standard uncertainty, r s wc uncertainty contributed the highest with 21.46 ± 5.97%, with VPD also having a considerable impact of 6.37 ± 1.45% while the rest of the inputs (r a wc , Fc and Fwet) contributed very little. Also, of the 0.33 mmday −1 transpiration uncertainty, 14.63 ± 8.71% of it was attributed to VPD uncertainty and 9.64 ± % to Fwet uncertainty. Wet soil evaporation uncertainty of 0.11 mm day −2 was made of 12.3% of Fwet uncertainty, 9.7% of VPD uncertainty, and very low contributions from the rest of the inputs. Lastly, VPD uncertainty contributed the highest to the potential soil evaporation uncertainty of 21.5%. In Welgegund (Table 7), ET estimation was only the sum of the potential soil evaporation and transpiration, and so were the uncertainties, since wet soil evaporation and wet canopy evaporation gave zero values and did not contribute to the final ET. The potential soil evaporation standard uncertainty of 1.05 mmday −1 was mainly a result of VPD, which contributed 46.11%, while other inputs had very low contributions. The transpiration standard uncertainty of 0.13 mmday −1 mostly resulted from the VPD, that contributed 26.93 ± 15.13%, while other inputs had relatively low contributions.
The total ET mean relative uncertainty for Skukuza was 76.19 ± 30.82%. The total uncertainty for Welgegund was similar to that of Skukuza, with a mean relative uncertainty of 81.1 ± 17.57%. In both sites, the highest uncertainty was attributed to soil evaporation, which contributed 76.74 ± 19.13% of the 1.38 ± 0.51 mmday −1 in Skukuza, and 90.93 ± 32.46% of the 1.62 ± 0.36 mmday −1 in Welgegund; subsequently, plant transpiration uncertainty presented a mean of 23.06 ± 18.83% and 18.21 ± 18.62% for Skukuza and Welgegund, respectively. On both sites, the wet canopy evaporation uncertainty was very low, which corresponded with this portion of evapotranspiration.

Sensitivity of PM-Mu Model to Core Input Variables
A sensitivity assessment of the ET output to input variables was done on the PM-Mu model to determine which input variable contributes the most to ET output variation. The percentage change in ET with respect to a percentage change in input variables at the study sites is summarised in Table 8 and illustrated in Figure 4. In the savanna biome, the PM-Mu model was mainly sensitive to LST and Tair. A change of −10% and −5% in Tair resulted in a 92% and 65% ET decrease, respectively, whereas an increase of 5% increased ET by 39%. In contrast, an LST decrease by the same values resulted in ET increasing by 55% and 40%, respectively. On the other hand, an LST increase of 5% and 10% resulted in an ET decrease of 51% and 77%, respectively. ε a changes from −20% to 0% gave an ET increase of +10%, and generated an ET increase and decrease of 12%, respectively, and α changes of −20% to +20% decreased ET by between 6.16 and −6.16%. ET had the lowest sensitivities to RH and LAI, with the computed parameter variations producing ET variations mostly inferior to 3%.
The grassland biome results were quite comparable to the savanna biome results, and ET was again mostly sensitive to Tair and LST. The percentage change in ET in relation to Tair decreases of −10% and −5% were −84% and −48%, while an increase of the same magnitudes showed an increase of 51% and 93%, respectively. Similarly, an LST decrease of the same magnitudes showed that ET increased by 85% and 44%, respectively, and an increase of 5% and 10% resulted in an ET decrease of 58% and 63%, respectively. ET was least sensitive to RH, LAI and α, with ET variations generally below 2.5%. Remote Sens. 2021, 13, x FOR PEER REVIEW 14 of 20 The grassland biome results were quite comparable to the savanna biome results, and ET was again mostly sensitive to Tair and LST. The percentage change in ET in relation to Tair decreases of −10% and −5% were −84% and −48%, while an increase of the same magnitudes showed an increase of 51% and 93%, respectively. Similarly, an LST decrease of the same magnitudes showed that ET increased by 85% and 44%, respectively, and an increase of 5% and 10% resulted in an ET decrease of 58% and 63%, respectively. ET was least sensitive to RH, LAI and α, with ET variations generally below 2.5%.

Discussion
In this study, we evaluated how input variable uncertainty was propagated in the PM-Mu algorithm to the final ET uncertainty, along with the analysis of the sensitivity of the ET output to the different input variable uncertainties. Our study only concentrated on the uncertainty associated with the input variables including uncertainty propagation, and not on the model algorithms used to compute the intermediate variables and the final ET product.
The measured meteorological input variable uncertainties were estimated as a combination of Type A and B uncertainties, whereas for the remote-sensing-based inputs, values from the literature were used to compute Type B uncertainties. Another essential assumption made was that there is no correlation between the inputs variables. This, in reality, might not be fully true, an example being the possible relationship between Tair and RH. A sensitivity analysis of ET to both the measured and remote sensing estimated input variables was also conducted.

Core Inputs
The mean absolute standard uncertainties for air temperature (Tair) were within similar ranges for the two sites, i.e., 1.5 ± 0.7 °C and 1.5 ± 0.42 °C, for Skukuza and Welgegund, respectively, with the variable showing slight seasonal uncertainty variation during the year. The cooler, drier season exhibited higher relative uncertainties compared with the hotter, wetter season on both sites. This is explained by high diurnal temperature ranges, and hence the high Type A standard uncertainties during the dry, cooler months. In contrast, relative humidity's (RH) absolute standard uncertainty was 5 ± 1.5% for Skukuza, and 6 ± 1.9% for Welgegund, showing a rather lower variability throughout the year compared to Tair, especially for Skukuza. This was probably due to more stable RH readings throughout the day, resulting in less variation in the estimated Type A uncertainty.

Discussion
In this study, we evaluated how input variable uncertainty was propagated in the PM-Mu algorithm to the final ET uncertainty, along with the analysis of the sensitivity of the ET output to the different input variable uncertainties. Our study only concentrated on the uncertainty associated with the input variables including uncertainty propagation, and not on the model algorithms used to compute the intermediate variables and the final ET product.
The measured meteorological input variable uncertainties were estimated as a combination of Type A and B uncertainties, whereas for the remote-sensing-based inputs, values from the literature were used to compute Type B uncertainties. Another essential assumption made was that there is no correlation between the inputs variables. This, in reality, might not be fully true, an example being the possible relationship between Tair and RH. A sensitivity analysis of ET to both the measured and remote sensing estimated input variables was also conducted.

Core Inputs
The mean absolute standard uncertainties for air temperature (Tair) were within similar ranges for the two sites, i.e., 1.5 ± 0.7 • C and 1.5 ± 0.42 • C, for Skukuza and Welgegund, respectively, with the variable showing slight seasonal uncertainty variation during the year. The cooler, drier season exhibited higher relative uncertainties compared with the hotter, wetter season on both sites. This is explained by high diurnal temperature ranges, and hence the high Type A standard uncertainties during the dry, cooler months. In contrast, relative humidity's (RH) absolute standard uncertainty was 5 ± 1.5% for Skukuza, and 6 ± 1.9% for Welgegund, showing a rather lower variability throughout the year compared to Tair, especially for Skukuza. This was probably due to more stable RH readings throughout the day, resulting in less variation in the estimated Type A uncertainty. Welgegund's RH relative uncertainties were substantially higher than Skukuza's uncertainties, indicating a higher diurnal variation of RH measurements at this site compared to the Skukuza site.
Our results are consistent with the ranges reported in other studies that have been conducted, albeit for different purposes. In most cases, Tair and RH uncertainties have been evaluated simultaneously. For instance, ascertained Muniz, et al. [52] the standard uncertainty of air temperature and relative humidity measured by thermography and found a standard Tair uncertainty of between 0 and 2 • C, and 0 and 5% for RH, in their study to ascertain the uncertainty of air temperature measured by thermography. In their case, though, they only considered Type B uncertainties, whereas we took into account Type A uncertainty, which is the variation of temperature over time. In addition, Lin and Hubbard [53] reported that the uncertainty of derived dewpoint temperature ranged from 1.8 to 3.3 • C. Comparing our results with other studies, like the ones discussed above, is a challenge because of the different metrics used in each study.
It is important to have an understanding of the uncertainties associated with the remote sensing products that are used in simulating ET, as shown in this study. These uncertainties are normally estimated and included in the ATB documents, with further research being carried out per biome. For example, the absolute quoted accuracy for LST is 1 K and 0.02 for εs, in the Modis ATBD document [39]. However, these accuracies vary with land cover type and the type of uncertainties included in the estimations [54] and should be investigated in detail.

Intermediate Data Components
The core input variable uncertainties had varying effects on the uncertainties of intermediate parameters like net solar radiation, and the aerodynamic and surface resistances. There was little variation in Rnet uncertainty between the grassland and savanna sites, with the relative uncertainties modelled at 4% and 3% for Skukuza and Welgegund, respectively. Our results showed that on both sites, LST uncertainty contributed the most to the Rnet uncertainty, with ε s and Tair uncertainties also contributing meaningfully to the Rnet uncertainty. Contrasting with other studies, our total Rnet uncertainties are much different. For example, we recorded much lower Rnet uncertainties than those reported by Mira, et al. [55], who reported overall uncertainties of between 15% and 20% in varying sites of rainfed to irrigated meadows and crops in the Mediterranean region of the Rhône Valley, in Southeastern France. They also found that the main contribution to the total uncertainty was equally distributed between the measured incoming short and longwave radiations (at 5% and 8%, respectively), with LST contributing the least uncertainty.
The aerodynamic resistance (r a wc , r a t , r a s ) relative uncertainties were consistent at an average of 0.8, 0.5 and 2%, respectively, throughout the year, at both sites. Considering these low uncertainties, their contribution to the ET uncertainties was also low. These results concur with the findings of Ershadi, et al. [56], who also showed that aerodynamic resistances play a relatively minor role in ET estimation in the PM model. Furthermore, it has been shown that changes in the parameterisation of aerodynamic resistance in the PM model produce minor improvements to the model output [57,58]. Our results show that in comparison to the aerodynamic resistance uncertainties, relative surface resistance uncertainties were quite high. Given that surface resistances (r s wc , r s t , r s tot ) are a critical parameterisation in ET estimation, the corresponding uncertainties also have a noteworthy impact on standard uncertainties of ET. This was also reported by Ershadi, et al. [56], who determined that surface resistance parameterisation substantially affects PM model performance.

Uncertainty in PM-Mu Evapotranspiration Estimation
The final ET uncertainty is a product of all the input variable uncertainties that were propagated through the PM-Mu model. In our study, we only investigated the uncertainty associated with the input variables, and not with the algorithms used to compute the intermediate inputs and the final ET model. Soil evaporation uncertainty contributed the most to the final ET uncertainty in our study areas, with wet canopy evaporation uncertainty contributing slightly less.
In both biomes, our results show that RH uncertainty, including input variables and parameters derived from RH, like VPD, Fwet and the different resistances contributed the most to the uncertainties of all the ET components. These results are in agreement with a study by Langensiepen, et al. [59], who reported that VPD is one of the principal meteorological variables in ET estimation using the PM model. Consequently, any error in the humidity and temperature measurement significantly impacts VPD uncertainty, and hence increases the overall uncertainty.
The overall mean relative ET uncertainty in our study was around 80% for both biomes, as measured from the propagation of input uncertainties. A similar approach was used to quantify the propagated ET uncertainty of several ET estimation methods, including the Penman method, in a riparian area of the Middle Rio Grande Basin, in New Mexico [5]. They used different values of input variable errors and obtained a much lower relative ET uncertainty, of only 10%, on the Penman method. This notable difference is due to a number of issues, including the methods and ecoclimates under investigation, the difference in the determination of the individual input uncertainties and the number of input variables and parameters considered in the propagation. Our uncertainties were a product of the U A and U B estimates propagated from the core inputs, through intermediate parameters, through to the final uncertainty product. In another study, predefined input uncertainty values from the manufacturers were used, for example, an Rnet uncertainty of 15%, while we reported a low Rnet uncertainty of 3 and 4% for our sites [5]. Ferguson, et al. [60] evaluated the uncertainty contributions of input variables to the final ET output. They show that the overall contributions of α and εs to ET uncertainties were minor, whereas LAI contributed quite substantially.

Sensitivity of PM-Mu ET Estimates to Input Variables
It is always a challenge to compare results on sensitivity analysis with other studies because of the difference between models, input variables and datasets, and the procedures being used to estimate the sensitivity coefficients. Moreover, these are applied under different ecoclimatic conditions under investigation. In our study, we used the simple one-at-a-time local sensitivity analysis method to estimate the sensitivity coefficients. Based on the maximum value of each input, we perturbed each of the input variables within the ±20% range. Our results show that PM-Mu is most sensitive to Tair and LST, thus making them the most influential input variables in ET estimation in southern Africa using the PM-Mu method. This result is explained by the fact that in dry environments the air has a high capacity to hold water vapour, which can then transfer energy to the land surface. Meanwhile, our results also show that the land surface parameters, like εa, LAI and α, have little effect on the PM method in these semi-arid ecoclimates. These results are consistent with other studies in similar semi-arid/arid regions, where the PM model was reported to be most sensitive to Tair [61,62]. A comprehensive sensitivity analysis assessment of the PM and Priestley-Taylor models to various inputs in different climates in Australia, Guo, et al. [26] showed that these models were most sensitive to Tair.
Other studies in similar dry climates contradict, however, our finding that PM-Muestimated ET is most sensitive to air and land surface temperatures. Tabari and Hosseinzadeh Talaee [63] observed that ET o was more sensitive to wind speed in a semi-arid climate, with less sensitivity to Tair and sunshine hours. Garcia, et al. [64] also reported that wind speed is a critical variable in ET o estimation in arid and semi-arid climates. They reasoned that this is because of the importance of the aerodynamic term under dry and high wind speed conditions. Gong, et al. [65] reported that the reference ET was most sensitive to RH variations, followed by solar radiation, in the Changjiang basin, in China.

Conclusions
In this study, we conducted a comprehensive uncertainty and sensitivity assessment of the PM-Mu model with regard to both meteorological and land surface characteristics in situ, i.e., air temperature and relative humidity, and remote sensing input variables, i.e., LAI, LST, ε s and α, in both savanna and grassland biomes of southern Africa. We only investigated the input variable uncertainties and quantified how these were propagated to the final ET uncertainty, and not the uncertainties related to the model algorithms used. With average relative uncertainties of 10 ± 3.5% and 6 ± 1.9% for Skukuza and Welgegund, respectively, RH, including its derivative, VPD, contributed the highest uncertainty to the final ET uncertainty. On the other hand, these results highlight that land surface and air temperatures and surface emissivity contributed the most to the solar net radiation uncertainty. Our results showed an overall ET uncertainty of approximately 80% in both biomes. This final propagated uncertainty is considerably larger than those reported in other studies. A number of reasons have been highlighted, including the number of input variables assessed for their uncertainty contribution, the assumption that there is no correlation between the input variables and the uncertainty analysis method used, that gives the total uncertainty as a combination of Type A and Type B uncertainties. This highlights the importance of accurate input data collection in ET estimation, as any errors are propagated to the final product. In contrast, the PM model was most sensitive to air and land surface temperatures, indicating the importance of these input variables to ET estimation using the PM-Mu model in our study areas. However, the sensitivity of the PM-Mu model to land surface parameters was quite low.