A Hybridized NGBoost-XGBoost Framework for Robust Evaporation and Evapotranspiration Prediction

We analyze the relationship between potential evapotranspiration (ETo), actual evapotranspiration (ETa), and surface water evaporation (Esw) in the semi-arid south-central Texas, using hourly climate data, daily lake evaporation measurements, and daily actual evapotranspiration measurements from an eddy covariance (EC) tower. The deterministic analysis reveals that ETo set the upper bound for ETa, but the lower bound for Esw in the study area. Unprecedentedly, we demonstrate that a newly developed probabilistic machine learning (ML) model, using a hybridized NGBoost-XGBoost framework, 5 can accurately predict the daily ETo, Esw, & ETa from local climate data. The probabilistic approach exhibits great potential in overcoming data uncertainties, in which 99% of the ETo, 90% of the Esw, and 91% of the ETa test data at three watersheds were within the model’s 95% prediction interval. The probabilistic ML model results suggest that the proposed framework can serve as a robust and computationally more efficient tool than the hourly Penman-Monteith equation to predict the ETo while avoiding computationally-involved net solar radiation calculations. Additionally, the performance analysis of the probabilistic 10 ML model indicates that it can be successfully implemented in practice to overcome the uncertainties associated with pan evaporation & pan coefficients in Esw estimates, and to offset the high capital & operational costs of EC towers used for Ea measurements. Finally, we demonstrate, for the first time, a coalition game theory approach to identify the order of importance, dependencies & interactions of climatic variables on the ML-based ETo, Esw, and ETa predictions. New knowledge gained through the game theory approach is beneficial to strategically locate weather stations for enhanced evapo(transpi)ration 15 predictions, and plan out sustainability and resilience efforts, as part of water management and habitat conservation plans.


Introduction
Evapo(transpi)ration is one of the key components of a groundwater budget in drought-prone regions with scarce water supplies (Heilman et al., 2009;Gokmen et al., 2013;Glenn et al., 2015), facing challenges of sustainable development and climate resilience. Reliable prediction of evapo(transpi)ration is useful in such regions to determine aquifer recharge (Hauwert and 20 Sharp, 2014; Xie et al., 2018), and subsequently, evaluate groundwater sustainability to meet municipal, agricultural, ranching, industrial, and recreational water demands, while sustaining quality & quantity of environmental flows to protect and maintain a healthy ecologic environment for endemic groundwater-obligated species.
the aquifer region to collect the most relevant data for enhanced evapotranspiration predictions and/or assess the suitability of simplified evapotranspiration prediction models for the watersheds with scarce data. Moreover, the interpretability of our ML model in combination with the game theory approach is capable of revealing new knowledge that may not be immediately apparent. For example, although soil moisture content was not included in our ML-model as a feature, but its effect was captured in ET a measurements at the EC tower, the ML model predicted low ET a despite high ET o & low RH in certain 70 times, which could be an indication of critical water deficiency in the soil. Such new knowledge from the ML model is essential for the current and future "well-informed" groundwater management and habitat conservation plans.

Methods
Description of Evapo(transpi)ration Measures. Different evapo(transpi)ration measures, including pan and lake evaporation, potential evapotranspiration, and actual evapotranspiration considered in this paper are briefly described here prior to 75 associated calculations and ML methods are introduced. For more comprehensive discussion, the reader may refer to the paper by McMahon et al. (2013).
Evaporation pans are used to determine evaporation from water surface at the pan-scale (E p ), which are then scaled-up to estimate evaporation from open water bodies (E sw ) such as lakes (Dingman, 1992). Therefore, lake evaporation can be interpreted as hybrid measured-estimated evaporation. E sw was viewed to represent regional potential evaporation (Vercauteren 80 et al., 2009) and has been used in terrestrial water balance calculations (Roderick et al., 2009). Daily or monthly empirical Meyer's formulas (MF) have been used to calculate E sw , based on surface water temperature, relative humidity, and wind speed measurements (Penman, 1948;Xu and Sing, 2002;Burn and Hesch, 2006). Potential evapotranspiration (ET o ), on the other hand, accounts for climate-driven watershed-scale evapotranspiration from a hypothetical reference crop in a saturated soil, which reflects the evaporation power of the atmosphere. The Penman-Monteith 85 equation (PME), based on the energy-balance, is used to calculate ET o (Allen et al., 1998). PME calculations require time series of shortwave solar radiation, air temperature, atmospheric pressure, relative humidity, and wind speed data. PME can be used for hourly to monthly ET o estimates, depending on the temporal resolution of the input climate data. ET o was used to estimate E sw (Vercauteren et al., 2009), actual evaporation (Boughton, 2004), vegetation potential evapotranspiration (Jia et al., 2009) or aridity index (Nash et al., 1997). The PME can be coupled with surface conductance models and leaf area 90 indices for evapotranspiration predictions. Using this approach, Zhang et al. (2008) reported good agreement between 5-year average evaporation rates predicted through the PME and using water balances at 120-gauged catchments in Australia. Leuning et al. (2008) noted reliable estimates of daily evapotranspiration rates at the kilometer-scale using the PME. PME calculations, however, are more complicated and involve more climate variables than empirical MFs. Therefore, simplified versions of the PME with fewer climate variables have been explored and tested for watersheds with scarce data Irmak et al., 95 2003; Peng and Feng, 2017).
However, neither ET o nor E sw provides a direct estimate for actual evapotranspiration (ET a ), which is the sum of evaporation from soil and transpiration from vegetation. As compared to ET o , ET a is more site-specific and spatially-variable, depending on soil and vegetation types. Reportedly, an increase in transpiration from vegetation could result in a two-fold decrease in soil evaporation (Yongqiang et al., 2016). Eddy covariance (EC) is the most direct method of measuring land surface 100 water vapor flux (Burba, 2013) without disturbing the water-air interface (Vesala et al., 2006), and hence, provides accurate site-scale ET a measurements (Wang et al., 2015). When coupled with the energy balance method, the EC technique provides an alternative measure of latent heat flux equivalent to ET a (Wilson et al., 2001;Zitouna-Chebbi et al., 2018). Shi et al. (2008) noted that PME resulted in higher latent flux than the EC method in estimating ET a of dry forest canopy. The relation between ET o and ET a is to be explored for the semi-arid region in this paper. 105 Penman-Monteith Equation (PME). A detailed description of underlying physical processes and calculation steps of the PME for hourly ET o estimates can be found in the FAO by Allen et al. (1998). This section provides the main equations and critical implementations for the solution of the PME for hourly ET o given by where σ is the Stefan-Boltzmann constant (4.903 × 10 −9 MJ / (K 4 m 2 d), T 4 a,max and T 4 a,min are the maximum and minimum 120 absolute air temperatures during the 24-hour period [K]. R so is the clear-sky radiation [MJ/(m 2 d)]. Linearized Beer's radiation law leads to R so = 0.75 + 2 × 10 −5 z R a , in which z is the elevation of the weather station above the sea level [m] and R a is the extraterrestrial radiation [MJ/(m 2 d)]. In other words, R so ∼ 0.75 of R a , which accounts for 25% reduction in Ra due to the interaction of R a with atmospheric gases Raza and Mahmood, 2018). (R s /R so ) is the relative shortwave radiation, representing the cloud cover, defined as in which the lower bound of 0.33 and the upper bound of 1.0 represent the dense cloud cover and clear sky on a particular day, respectively. The first, second, and third terms in Eq. 2 account for the effect of air temperature, air humidity, and cloudiness on R nl . R a depends on the geographic location of the weather station and time of the day, and is computed as In hourly ET o calculations, R a = 0 when the sun is below the horizon at ω < −ω s or ω > ω s . To keep the cloudiness, R s /R so in Eq. 3, and hence, R nl in Eq. 2 finite, R s /R so at night times (i.e., when the sun is below the horizon) is set to R s /R so 140 value 2-3 hours prior to sunset. The sunset time in each day of the year can be identified by (ω s − 0.79) ≤ ω ≤ (ω s − 0.52).
When the sun is above the horizon (R a > 0), G = 0.1R n corresponds to smaller heat outfluxes, promoting soil warming during day times. In contrast, when the sun is below the horizon (R a = 0), G = 0.5R n corresponds to larger heat outfluxes, promoting soil cooling at nights. Moreover, wind speed, u 2 ≥ 0.5 m/s in ET o calculations to account for the effects of boundary layer instability and buoyancy of air in promoting exchange of vapour at the surface when air calm. 145 Meyer's Formula (MF). Meyer's formula (MF) is a mass transfer-based, empirically constructed formula (Meyer, 1915), whose different versions have been used to calculate daily or monthly E sw . It is typically expressed in the form of E sw = β (e o − e a ), in which β is the empirically determined constant, e o and e a are defined in terms of surface water temperature, unlike in the PME. Reportedly, the best form of the MF to predict daily E sw from free water surface constructed using data from England (Penman, 1948;Xu and Sing, 2002) 150 where u 2 is expressed in [mm/d], and e o and e a are expressed in [mm-Hg] in Eq.6. For monthly evaporation estimates from surface of a water body using data from Canadian prairies (Burn and Hesch, 2006), where u 7.56 is the monthly-averaged wind speed measured at 7.56 m above the ground surface, expressed in [km/hr], and e o 155 and e a are expressed in [mbar] in Eq.7. C = 1 in the original work by Burn and Hesch (2006). C = 1 is introduced here to adjust the magnitude of monthly E sw . When compared to PME in Eq. 1, MFs involve fewer climate variables and do not involve computationally-involved net solar radiation calculations.  XGBoost, proposed by Chen and Guestrin (2016), is a tree-based ensemble learning algorithm that follows the principle of boosting. Boosting is a general technique in ML, where multiple weak learners such as Classification and Regression Trees (CART) are organized to produce a strong learning model (Marsland, 2014). The fundamental concept behind this technique is to produce new learners that are sequentially fitted to the residuals from the previous learner, which are then added to the model to update the residuals. Gradient boosting enhances the flexibility of the boosting algorithm by generating the new learners that 170 are maximally correlated to the negative of the gradient of the loss function. This process enables the convergence of the loss function and allows arbitrary differentiable loss functions to be used in the model building process (Chen, 2014;Chen and He, 2015). From the computational standpoint, XGBoost is built with a multiprocessing OpenMP API (Chandra et al., 2001), which enables XGBoost to use all the CPU cores in parallel during while training, making it computationally efficient and scalable. Moreover, XGBoost presorts the independent variables at the beginning of the training process, which further reduces 175 the training complexity and computational time.
NGBoost, proposed by Duan et al. (2019), is a supervised learning algorithm with generic probabilistic prediction capability.
A probabilistic prediction produces a full probability distribution over the entire outcome space; thus, enabling the users to quantify the uncertainties of the evapotranspiration predictions produced by the model. In standard point prediction settings, the object of interest is an estimate of the scalar function E(y|x), where x is the feature vector and y is the prediction target, without accommodating uncertainty estimates. In contrast, under a probabilistic prediction setting, a probabilistic forecast with probability distribution P θ (y|x) is produced by predicting the parameters θ. NGBoost can perform probabilistic forecast with flexible tree-based models, given that NGBoost is designed to be scalable and modular with respect to the base learner (e.g. decision trees), probability distribution parameter (e.g. normal, Laplace), scoring rule (e.g. Maximum Likelihood Estimation).
We utilized NGBoost's modular design to hybridize it with XGBoost base learners to enhance the resulting model's pre-185 dictive capability. As shown in Fig. 1, the input feature vector x in the hybrid NGBoost-XGBoost model is passed on to the XGBoost base learners to produce a probability distribution of the predictions P θ (y|x) over the entire outcome space y (i.e., ET o ). The models are then optimized by scoring rule S(P θ , y) using a maximum likelihood estimation function that yields calibrated uncertainty and point predictions. The feature vector x for ET o and ET a predictions consists of T a , P , RH, u 2 , and R s ; and the feature vector x for E sw predictions consists of T sw , T a , P , RH, u 2 , and R s .  Local climate data at the NDR weather station. For hourly-ET o calculations, hourly-averaged T a , P , RH, and u 2 and hourly-summed R s at the NDR station, shown in Fig. 3, were used as input in Eq. 1. The total number of missing hourly records was 2, which were filled in by linear interpolation. The NDR weather station was selected in the analysis due to its 200 proximity to Uvalde County, TX, where monthly representative cloud cover data was available, which were used to test the model accuracy in Section 3. Local climate data at the BCRA and CBS stations, along with ET a measurements at Savanna, Well 10 are provided in Appendix A.
Surface Water Data. Daily and monthly surface water evaporation data closest to the NDR site were obtained from Ingram Lake in Texas. Daily pan evaporation measurements (E p ) from 9/1/2015 to 12/31/2019 were taken by the Texas Water De- velopment Board (TWDB). 1.9% of these measurements were missing, which were filled in by linear interpolation. These measurements were upscaled to daily lake evaporation totals (E sw ) using monthly-varying pan coefficients developed by the TWDB. However, sporadically extremely high and low E sw values, shown in Fig. 4a, were found to be quantitatively inconsistent with the climatic data (T a , R s , and RH) trends at the BCRA station, provided in Appendix A. Therefore, this time series is regarded as anomalous. Such anomalies are quite common in E p measurements due to birds drinking from the pan, 210 debris falling in, or water splashing out (Thompson, 1999). Subsequently, these anomalies are carried into the daily E sw data, but largely smoothed out in monthly-averaged E sw . Because the ML model was run with daily E sw data here, a 7-day rolling median function was used to reduce the noise and outliers in the daily E sw data (Fig. 4a). Monthly E sw , derived from daily E sw (Fig. 4) were then used to determine the suitability of the MFs to predict the monthly E sw at Ingram Lake.
The daily and monthly MFs rely on surface water temperature, T sw , rather than T a , in e o and e a calculations (Eqs. 6-7). The 215 closest gauging station, with the surface water temperature data from 9/1/2015 to 12/31/2019 at the 15-min (or 1-hr) intervals, to Ingram Lake is the U.S. Geological Survey Station (USGS 08195000) located at the Frio River in Concan, TX. The Frio River at the USGS 08195000 and Ingram Lake are small-size surface water bodies fed by groundwater from the Trinity aquifer. Therefore, T sw from the Frio river were used in MF-based E sw calculations at Ingram Lake. Daily-averaged T sw are shown in Fig. 4b. Because the Frio river is a groundwater-fed river, T sw ≥ T a in winter; whereas, T sw ≤ T a in summer. 0.26% of daily 220 T sw were missing, which were filled in by linear interpolation. Hourly ET o were calculated via Eq.1, using local climate data at the NDR weather station from 9/1/2015 to 12/31/2019. In Fig. 5a, R a , computed using Eq. 4, relies on information on the geographic location of the weather station and hourly-varying solar time angle. R a was subsequently used to calculate R so .
The ratio of the measured R s to the computed R so provides an estimate for cloudiness, defined as the fraction of the number of 230 cloudy-sky hours in a day. In Fig. 5b, PME-computed monthly-averaged cloudiness from 2016 through 2019 agrees well with the monthly-averaged representative cloudiness for the Uvalde city.  10 https://doi.org/10.5194/hess-2020-260 Preprint. Discussion started: 13 July 2020 c Author(s) 2020. CC BY 4.0 License.
Measured R s was the input to the PME. R nl in Fig. 5c was computed by Eq. 2, based on daily extremes of air temperature, and hourly cloud cover and actual vapor pressure. R nl was then used to calculate the net solar radiation, R n = R nl − R s , in Eq.1. Hourly ET o were aggregated to daily ET o , which is the time interval at which both E sw and ET a measurements 235 were available. Although hourly ET o contains negative values at the humid and rainy hours, daily ET o 's were persistently non-negative (Fig. 5d), as expected for the Texas climate.
Lake Evaporation Using Meyer's Formula. Daily and monthly E sw data from Ingram Lake, the closest water body to the NDR site, have been reported by the TWDB. Applicability of Eqs. 6 and 7 to predict monthly E sw were tested here using monthly E sw data. In this analysis, E sw at Ingram Lake computed using Eq. 6 averaged over each month to obtain monthly 240 E sw . Local climate data from two weather stations, the NDR station (∼ 90 km away) and the BCRA station (∼ 40 km away), were used in calculations.
In Fig. 6a, monthly-averaged daily E sw computed by Eq. 6 matched the E sw data almost perfectly (R 2 = 0.99) when the climate data at the NDR station was used, but underpredicted the E sw when the climate data at the BCRA station was used. As shown In Fig. 6b, the original form of Eq. 7 with C = 1 matched the overall monthly trend of the E sw data, but underestimated 245 the magnitude of E sw irrespective of climate data from the NDR or BCRA station. When C = 1.6 with the BCRA data (and C = 1.5 with the NDR data), Eq. 7 matched the monthly E sw data almost perfectly (R 2 = 0.99). Although the empirical relations in Eqs. 6 and 7 were independently derived using site-specific data at two markedly different geographic locations in Canada and England, these equations matched monthly E sw in south-central Texas surprisingly closely in Fig. 6.  Computed Potential Evapotranspiration vs. Lake Evaporation. Fig. 7 shows that ET o , in general, set the lower bound for 250 in the summer of the following years, with the largest difference in the summer of 2018, ET o appears to be a reliable predictor for E sw especially in spring and winter months. Fig. 7 reveals that ET o computed by the PME using the climate data from the NDR or BCRA stations can be used to estimate the minimum monthly E sw , computed by MF and determined by TWDB, for Ingram Lake. This is rather an interesting result, 255 given that the empirical MF does not involve complex solar radiation (e.g., extra-terrestrial, clear-sky, outgoing longwave) calculations as in the PME and relies on surface water temperatures rather than air temperatures. It can be argued that surface water temperature may implicitly account for the solar radiation effect on E sw . Nevertheless, the main conclusion from Fig. 7 is that ET o at the NDR station set the lower bound for E sw at nearby Ingram Lake (i.e., ET o ≤ E sw ). Fig. 7 further suggests that if uncertainty in local climate measurements are higher than lake evaporation measurements, mathematically simpler MF, 260 after being validated with historical lake water evaporation data, can be used to predict potential evapotranspiration from new lake evaporation data. Fig. 7 also shows that monthly-aggregated ET o near Ingram Lake is slightly higher in summers than at the NDR or BCGA sites. However, if the climate data near Ingram Lake is not available, data from a farther weather station to the west can be used to predict the ET o , and hence, the lower bound for E sw . Based on the existing data, the results also suggest that no additional weather station is needed between the NDR and BCGA stations to predict ET o and/or E sw from 265 other surface water bodies fed by the same groundwater system between these two stations. Potential Evapotranspiration vs. Actual Evapotranspiration. Fig. 8 compares daily or monthly Bowen-ratio-corrected ET a measurements from the EC tower against ET o computed by the PME, using local climate data from the CBS station. According to this plot, ET o ≥ ET a during the monitoring period, as expected. In some summer months, ET o was about three times higher than ET a (e.g., July 2017), indicating that the soil was too dry in summer times to contribute to the evapotranspiration at the CBS site. As compared to ET o at the NDR site, ET o is typically higher at the CBS site in summer months, revealing spatial variability in ET o . The results in Figs. 7 and 8 lead to  ET a is the most critical evapotranspiration estimate, especially for irrigation, agricultural, and water resources management practices. As discussed previously, the EC method provides the most accurate prediction for ET a ; however, the associated 275 capital and maintenance costs are high (e.g., the capital cost for the EC tower at the CBS site was about $40,000 and required frequent maintenance The ML-based E sw prediction model was trained by using the first 90% of the daily climate data & the month of the year as features, and the measured E sw data as the target. Subsequently, the model was tested on the remaining 10% of the data. The comparison between ML-predicted daily E sw and the measured E sw on the testing data is shown in Fig. 9b. The ML-305 predicted E sw matched the measured E sw very closely, and ∼ 90% of the actual E sw were within the model's 95% prediction interval in the testing dataset. Based on the statistical measures in Table. 1, probabilistic prediction of E sw by the hybrid NGBoost-XGBoost ML model can be used as a reliable method for E sw projections. The total training time for the E sw hybrid model was ∼ 6 minutes, including choosing the optimum model from 690 model fits. The main advantage of the ML-based E sw prediction model is that E sw predictions are not affected by anomalies in E p measurements (Fig. 4a)   The ML-based ET a prediction model was trained by using the first 90% of the daily R s , PME-computed ET o , & the month of the year as features, and the measured ET a data as the target. Subsequently, the model was tested on the remaining 10% of the data. The comparison between ML-predicted daily ET a and the actual ET a measurements -on the testing data -is shown in Fig. 9c, in which ∼ 91% of the actual ET a values were found within the model's 95% prediction interval. ML-based

315
ET o predictions was more accurate than the ET a predictions due, in part, to the availability of less data from the EC tower at the CBS site than at Ingram Lake or the NDR site for the ML-model training. However, based on the statistical measures in Table. 1, probabilistic prediction of ET a by the hybrid NGBoost-XGBoost ML model, with an R 2 of 0.804 on testing data, can still be used as a reliable method to estimate the future ET a . The total training time for the ET a hybrid model was ∼ 9 minutes, including the 690 model fits to choose the optimum model. The main advantage of the ML-based ET a prediction 320 model is that it offsets the high capital and maintenance costs for the installation and operation of EC towers to acquire ET a measurements. Obviously, if more E sw and ET a measurements are available for ML-model training, the predictive accuracy of the respective ML models would improve. The EAA is planning to construct additional EC towers in other parts of the aquifer region to collect more ET a data, which would enhance the training and predictive performance of the ML-model. Similarly, as the TWDB continues to collect E sw data by effectively filtering out observed anomalies, availability of longer E sw data with 325 less noise for ML-model training would enhance the predictive accuracy of the model.  This finding is important to evaluate the suitability of the simplified versions of the PME proposed for semi-arid watersheds with scarce climate data. Irmak et al. (2003) proposed two simplified PMEs that require less number of climate variables to 340 calculate the net radiation (R n in Eq. 1). The first equation relied on the measured T a and R s , whereas the second equation relied on predicted R s , and measured T a and RH. Although the simplified equations were used to estimate R n only, the second equation, built on the three most important climate variables identified in Fig. 10a   R s was recently used as a surrogate variable to reduce the uncertainty of ET o projection data (Yoo et al., 2020). This can be justified by the findings from Fig. 10a, in which R s displayed a more profound impact on ET o than the other forcing variables.
On the other hand, the mean annual temperature was used by Hartmann et al. (2017) as a proxy for ET o in assessing aquifer recharge sensitivity to climate variability based on the argument that R n is temperature-dependent and temperature is the best-understood and most common climatic variable for large-scale hydrological models. Similarly, a computationally simple 350 method of Berti et al. (2017) that relies only on T a was reported to be the best alternative method to the PME in describing spatiotemporal characteristics of ET o in different sub-regions of mainland China (Peng and Feng, 2017). Such assumptions (Hartmann et al., 2017) and conclusions (Peng and Feng, 2017) Figure 11. Correlation map between daily climatic variables and (a) the potential evapotranspiration at the NDR site, (b) lake evaporation at Ingram Lake, and (c) and actual evapotranspiration at the Camp Bullis site. Gong et al. (2006) noted that although the order of importance of climate variables on ET o estimates through the PME varied with season and region in their study, ET o in general was most sensitive to RH, followed by R s , T a and u 2 . The authors used time-histories of daily T a , u 2 , RH, and daily sunshine duration. In our analysis, however, measured climate variables were available at the 15-min intervals, including also R s and P . Unlike the general conclusion by Gong et al., our ML-based feature 360 importance calculations in Fig. 10a revealed that both R s and T a were more critical than RH on ET o estimates. Fig. 10b shows that the E sw is largely impacted by the T sw followed by RH. T a is given lower importance because of the high correlation (R 2 = 0.95) between T sw and T a (Fig. 11b), and thus, the model considers T a as redundant. Fig. 10b also highlights the model's understanding of the underlying hydrological process. For example, we see that the model tries to push the E sw predictions upward -represented by higher Shapley value on the x-axis -when the T sw feature values are high -(a) RH dependence plot with ETo interaction.
(b) ETo dependence plot with RH interaction. Eddy covariance (EC) towers provide accurate estimates for site-specific ET a , but at the expense of high capital, operational, 400 and maintenance costs, which often limit the use of multiple EC towers and/or their operational periods in resources management projects. Pan evaporation method, on the other hand, is a simple, inexpensive, and widely-used data acquisition method to predict E sw at open water bodies, but suffers from uncertainties in pan evaporation measurements and in pan coefficients for the existing or projected climatic conditions when water evaporation is upscaled from the pan-scale to the large openwater body-scale. ET o is often computed by energy-balance models, such as Penman-Monteith equation (PME) that relies on 405 time-series of more local climate variables and includes rather complex calculations for net solar radiation computations. In brief, ET a measurements are often challenged by the project budget; whereas, E sw measurements are affected by uncertainties in and upscaling of pan-evaporation measurements. ET o calculations, on the other hand, require computationally-involved calculations.
To eliminate these shortcomings in ET a , E sw , and ET o predictions, we proposed a hybrid ML probabilistic prediction 410 model of ET o , E sw , and ET a using the local climate data and the month of the year as the only independent feature. Different from other ML models, the proposed hybrid ML model is able to produce point predictions as well as a probability distribution over the entire outcome space for quantifying the uncertainties related to hydrological predictions. The proposed hybrid model could provide practitioners with a better understanding of the uncertainty in the ET o , E sw , and ET a predictions without compromising the accuracy of the predictions. Our results showed that the hybrid (NGBoost-XGBoost) ML model successfully 415 predicted the PME-computed ET o , and the measured E sw and ET a , in which ≥ 90% of the target data points were within the 95% prediction interval of the model, and R 2 values for the point predictions were 0.99, 0.75, and 0.8, respectively, using data from the model testing period. These results exhibit that the proposed hybrid ML model is a reliable and robust alternative method to predict ET o , E sw , and ET a from local climate data, without implementing computationally-intensive PME calculations, or coping with uncertainties in E sw estimates using evaporation pans, or having expensive EC tower setups 420 for ET a measurements.
We also demonstrated that the hybrid ML model, based on a game theory approach, generated new knowledge, different from sensitivity analyses findings in the existing literature, about the importance of features (variables) on the ET o , E sw , and ET a predictions. The underlying idea behind this analysis was to explain the prediction of an instance by computing the contribution of each feature to the prediction. Our analysis revealed that the shortwave solar radiation, air temperature, and relative humidity 425 are the most critical features for the ET o predictions, whereas the surface water temperature, relative humidity, and the month are the most critical features for the E sw predictions, and the shortwave solar radiation, month, and relative humidity are the most critical features for the ET a predictions in the semi-arid climate.
The EAA has 12 active weather stations across the Edwards aquifer region. The proposed hybrid ML models would allow continuous ET a predictions, without the need for expensive EC tower setups, from continuously streaming climate data at these 430 weather stations and through their interpolation between the stations. Moreover, the ML model would allow E sw predictions from surface water bodies without equipped with sensors or tools for E sw measurements, if they are located closer to the weather stations. Thus, the ML-model would curtail data acquisition costs and ML-based ET a predictions would be particularly useful for real-time aquifer recharge estimates, and irrigation and agricultural water management.
The proposed hybrid ML model would also be a useful tool to project local evapotranspiration and its impacts on aquifer 435 recharge when projected local climate data over the next 30-50 years is statistically-downscaled from global climate/regional climate data and reinterpreted as a time series of local climate data. Such projections are crucial to predict potential drought of records ahead of time before it results in irreversible damages on the sustainability of groundwater resources for diverse consumptive uses and habitats of groundwater-bound endangered or threatened aquatic species.
Data availability. Data are available from the authors.