Estimating Irrigation Water Consumption Using Machine Learning and Remote Sensing Data in Kansas High Plains

: Groundwater-based irrigation has dramatically expanded over the past decades. It has important implications for terrestrial water, energy ﬂuxes, and food production, as well as local to regional climates. However, irrigation water use is hard to monitor at large scales due to various constraints, including the high cost of metering equipment installation and maintenance, privacy issues, and the presence of illegal or unregistered wells. This study estimates irrigation water amounts using machine learning to integrate in situ pumping records, remote sensing products, and climate data in the Kansas High Plains. We use a random forest regression to estimate the annual irrigation water amount at a reprojected spatial resolution of 6 km based on various data, including remotely sensed vegetation indices and evapotranspiration (ET), land cover, near-surface meteorological forcing, and a satellite-derived irrigation map. In addition, we assess the value of ECOSTRESS ET products for irrigation water use estimation and compare with the baseline results by using MODIS ET. The random forest regression model can capture the temporal and spatial variability of irrigation amounts with a satisfactory accuracy ( R 2 = 0.82). It performs reasonably well when it is calibrated on the western portion of the study area and tested on the eastern portion that receives more rain than the western one, suggesting its potential transferability to other regions. ECSOTRESS ET and MODIS ET yield a similar irrigation estimation accuracy.


Introduction
Groundwater-based irrigation accounts for more than 70% of global freshwater use [1], especially in the arid and semi-arid regions [2,3]. The High Plains is one such region where groundwater-based irrigation has dramatically increased over the last decades and led to groundwater and streamflow depletions [4][5][6]. It is estimated that the cumulative groundwater depletion in the High Plains Aquifer (HPA) in the last century is over 353 km 3 , which is equivalent to a 1 mm sea-level rise [7]. Groundwater irrigation introduces additional water to the land surface and reduces the surface air temperature (the cooling effect) by shifting sensible heat to latent heat [8][9][10]. The latent heat, or evapotranspiration (ET), is further enhanced by the augmented available energy as a result of increased land greenness and decreased albedo of irrigated vegetation [10][11][12][13]. The extra atmospheric water vapor from the groundwater stimulates cloud cover and precipitation [9]. Nearly 80% of the extracted groundwater ends up in the ocean by surface runoff and increased precipitation caused by atmospheric vapor transport from the land to the ocean [14], which ultimately results in sea-level rise. On the other hand, groundwater depletion can also modify the monsoon circulation by changing the thermal contrast between the land and the ocean [15].
Despite the important implications of agricultural irrigation, discrepancies exist between the observations and numerical models in quantifying this human intervention in water cycles [16]. For example, Liu et al. [17] suggest a high bias (i.e., overestimation) in the temperature in the central US simulated by the Weather Research and Forecasting coupled to the Noah-Multiparameterization model (WRF/Noah-MP), especially during summer, likely due to the lack of representation of irrigation activity in the model. Scanlon et al. [18] indicate that the modeled terrestrial water storage (TWS) of seven global hydrologic models shows an opposite trend to that derived from the Gravity Recovery and Climate Experiment (GRACE), mainly due to the inability of the models to capture the human-induced water storage change. Adding an irrigation scheme to Noah-MP reduces the gap between the model simulation and the remotely sensed ET, GRACE TWS, and groundwater level measurements for the southern HPA, highlighting the importance of the enhanced representation of human water use [19].
The growing repository of in situ and remote sensing observations and data fusion techniques provide opportunities to develop a retrospective dataset of groundwater irrigations [20]. Two commonly used types of remote sensing data for irrigation estimation are soil moisture (SM) and ET. SM products from Soil Moisture Active Passive (SMAP), Soil Moisture and Ocean Salinity (SMOS), Advanced SCATterometer (ASCAT), and Advanced Microwave Scanning Radiometer (AMSR2) have been used to detect irrigation areas [21][22][23][24][25]. Methods have also been developed to quantify irrigation amounts using SM data. For example, Brocca et al. [26] calculate irrigation as the difference between storage change (from remotely sensed SM) plus outflow (drainage and ET calculated using an SM balance model) and in situ rainfall measurements at pilot sites located in the U.S., Europe, Africa, and Australia. Other studies extend this method with various estimated and remote sensing ET products [27][28][29]. Zaussinger et al. [30] infer the monthly irrigation water use over the contiguous United States (CONUS) based on the gap between remote sensing SM and Modern-Era Retrospective Analysis for Research and Applications Version 2 (MERRA-2) reanalysis SM that does not include irrigation. Zappa et al. [31] further incorporate calculated ET and drainage loss for their study sites in Germany into this method. Studies have also estimated the irrigation amount by comparing precipitation to remotely sensed ET, such as the MODIS (Moderate Resolution Imaging Spectroradiometer) product (MOD16A). More specifically, MODIS ET, which is treated as the ground truth crop ET rate, surplus over the effective precipitation or Land Surface Model (LSM) simulated ET, is treated as a crop consumptive water use supplied by irrigation [27,[32][33][34][35][36].
The above approaches based on water balance rely heavily on the accuracy of the remote sensing and reanalysis products. Due to the coarse resolution and uncertainty of remote sensing SM [26,27], these approaches often perform less well for low-intensity irrigation signals, such as supplementary irrigation, irrigation in semi-humid and humid areas, and irrigation with high-efficiency technology [26,27,30]. In parallel, various irrigation schemes have been developed within LSMs to derive irrigation driven by crop demands [19,37]. However, these LSMs' capability of predicting SM is questionable considering the uncertainties in the upper (e.g., infiltration of ponding water and runoff) and lower boundary conditions (e.g., groundwater exchange).
Irrigation water management is heavily dependent on socioeconomic-controlled irrigation behavior, including local policy, water rights, and farm producer irrigation preferences [37][38][39][40]. Irrigation schemes implemented in LSMs often fall short of representing the socioeconomic controls on irrigation behavior. Machine learning algorithms and remote sensing data provide opportunities to capture human-water interactions that are hard to represent using physically based methods [39,[41][42][43][44][45]. Several studies have applied machine learning approaches in the High Plains region to map irrigated farmlands [46], identify irrigation driving factors [38], and estimate groundwater withdrawal from MODIS ET [47].
We present a machine learning approach for estimating spatial irrigation amounts by fusing hydrometeorological, in situ, and multi-source remote sensing data that reflect dynamic crop conditions. Our approach is applied over the Kansas HPA, where in situ pumping records [48] are available for training and validating the machine learning model. Through a spatial prediction experiment, we demonstrate the potential transferability of our approach to other regions with scarce or no in situ pumping records to create an annual, gridded irrigation amount dataset, which can then be used to improve the representation of anthropogenic irrigation within ecohydrological and atmospheric models and enhance the prediction capability of regional land-atmosphere feedbacks. In addition, we compare the irrigation estimation accuracy using the recently available ECOsystem Spaceborne Thermal Radiometer Experiment on Space Station (ECOSTRESS) ET products and MODIS ET.

Study Area
Our study area is a portion of the central HPA in Kansas (Kansas HPA, Figure 1a). In recognition of the "lead from local need" perspective, five groundwater management districts (GMDs) were established in the 1970s (Figure 1b) for water use administration and planning [49]. The primary crop types in this region are winter wheat, corn, soybeans, sorghum, and alfalfa [50] (Figure 1c). The study area experiences a climatic gradient from semi-arid in the west to sub-humid in the east and a decreasing temperature from south to north (Figure 1d,e; [51]). The agriculture of Kansas HPA heavily relies on irrigation to meet crop water demands. The Kansas Geological Survey Water Information Management and Analysis System (WIMAS) [48] dataset provides a continuous annual irrigation water use record reported by individual water right holders since 1990. For over 95% of the irrigation wells in Kansas HPA, the reported water use volume is measured by totalizing flowmeters [51]. The non-metered wells are required to provide pumping hours and pumping rates to ensure the reported data quality [52].

Data Collection
We compiled various remote sensing products, climate data, and ancillary information (Table 1, Figure 2) to estimate irrigation using a machine learning model. All data were aggregated into a 6 km × 6 km spatial resolution, which can be readily adapted to

Data Collection
We compiled various remote sensing products, climate data, and ancillary information (Table 1, Figure 2) to estimate irrigation using a machine learning model. All data were aggregated into a 6 km × 6 km spatial resolution, which can be readily adapted to the model's needs. We further aggregated time-varying remote sensing and climatic data to an annual scale. For these data sources, we chose to focus on 1 May-31 August (MJJA), which covered most of the growing season, according to the planting and harvesting schedule of primary crop types in the study area [53,54]. The remote sensing data include the Normalized Difference Vegetation Index (NDVI), Normalized Difference Water Index (NDWI), and Greenness Index (GI), and ET products. The vegetation indices (NDVI, NDWI, GI) were calculated from 30 m-resolution Landsat 5, Landsat 7, and Landsat 8 surface reflectance images, respectively. Due to the differences in band designations and sensors [55,56], we normalized the vegetation indices from each Landsat image by linearly scaling them to the range of [0, 1] before aggregating to the annual scale. We used the following three ET products: MODIS 8-day 500 m resolution ET (MOD16) [57], ECOSTRESS PT-JPL instantaneous ET [58], and ECOSTRESS DisALEXI-JPL daily ET [59]. Both DisALEXI-JPL and PT-JPL products are at a 70 m resolution using the ECOSTRESS Level 1 geolocation information, Level 2 Land Surface Temperature (LST) and Emissivity product, near surface temperature, wind speed and relative humidity from National Centers for Environmental Prediction (NCEP), and data from MODIS and LANDSAT 8. The ECOSTRESS DisALEXI-JPL daily ET is based on an Atmosphere-Land Exchange Inverse (ALEXI) surface-energy balance model. The ECOSTRESS PT-JPL is based on the Priestley-Taylor (PT) method [60] and available both as instantaneous (i.e., at the overpass time) and daily rates. We used the instantaneous product because the daily product is interpolated from the instantaneous ET assuming a clear sky diurnal solar cycle. As a result, the daily ECOSTRESS PT-JPL ET product tends to overestimate, especially in frequently cloudcovered regions [61].

Methods
We used random forest regression (RFR) to estimate the irrigation amount from the input datasets (Table 1). The irrigation amount was represented as depth by dividing the total volume of the irrigation water use within each 6 km × 6 km grid by the area of the 2 Figure 2. Flowchart of a machine learning-based estimation of the spatial irrigation amount by fusing hydrometeorological data, in situ irrigation water use records, multi-source remote sensing data, and derived irrigation maps and crop fractions. The second group of inputs was calculated from GRIDMET climate data [62], including daily maximum temperature, total precipitation, and daily mean water vapor deficit (VPD). The third group of inputs were derived products of land use and land cover, including annual irrigation maps for HPA (AIM-HPA) [46] and crop fractions. The AIM-HPA dataset identifies irrigated fields at a 30 m resolution and was derived based on Landsat imagery, hydroclimatic data, and ancillary information, such as the USDA Soil Survey Geographic Database (SSURGO). The irrigation maps were resampled to the 6 km resolution used in this study by calculating the percentage of the irrigated area within one 6 km × 6 km grid. Lastly, we calculated crop fractions based on the USDA National Agricultural Statistics Service Cropland Data Layers [63] (available since 2006) as the percentage of the cultivated area of each primary crop type in a 6 km × 6 km grid.

Methods
We used random forest regression (RFR) to estimate the irrigation amount from the input datasets (Table 1). The irrigation amount was represented as depth by dividing the total volume of the irrigation water use within each 6 km × 6 km grid by the area of the grid (36 km 2 ). A random forest is an ensemble of decision trees with bootstrap sampling [68]. By averaging the results across all trees and randomly choosing a subset of variables at each split point, the RFR is less prone to overfitting than a single tree [68]. RFR has been used in many remote sensing applications as it can manage high dimensional, correlated data with easy hyperparameter tuning and a relatively fast process [69]. In this study, we trained and validated RFR models using the in situ irrigation records as targets. We also used variable importance scores to examine the contribution of each variable to the estimation accuracy. The variable importance was calculated based on the Gini impurity, measuring how splitting with an input variable decreases the impurity (variance) of the data [70]. The performance of the model was evaluated by root-mean-square error (RMSE), coefficient of determination (R 2 ), and percentage bias (PBIAS) statistics calculated using in situ irrigation records reserved for the validation. The PBIAS measures the average trend of the predicted value to be larger or smaller than the actual value [71].
We performed three experiments to test the prediction skill of the RFR model, probe its potential for transferability, and investigate the utility of the ECOSTRESS ET for an annual irrigation amount estimation ( Figure 2). In Experiment A, we randomly split the data across the study area from 2006 to 2020 into 70% for training and 30% for test. To examine the capability of the RFR model to capture the overall spatial and temporal variability of the irrigation amount, we trained the RFR model with all years' training data and assessed the accuracy of the trained model on the validation data points each year.
In Experiment B, we split the data from 2006 to 2017 by region. More specifically, we divided the model domain into the western and eastern portions for training and test, respectively. Given the differences in climate conditions and crop types between the training and test regions (Figure 1c-e), this experiment represents a challenging case and was designed to investigate the spatial transferability of our approach. Since the ECOSTRESS ET product was not available until late 2018, Experiments A and B were performed using MODIS ET.
In Experiment C, we performed the same random split test for years 2019 and 2020 using the ECOSTRESS in one run and the MODIS ET in another run. We compared the accuracy resulting from the use of the ECOSTRESS and MODIS ET products.

Experiment A: Random Split
Overall, the RFR achieves satisfactory performance, with an average RMSE, PBIAS, and R 2 of 22.9 mm, −0.73%, and 0.82, respectively ( Figure 3). Despite interannual variabilities, a declining trend in R 2 is noticeable primarily due to the decreasing variability of the irrigation amount among the grids in the study area ( Figure S2). The decrease in variability suggests that the irrigation behavior became more homogeneous. Possible causes may be water use being constrained by local water conservation programs and the expansion of water-saving irrigation technology during the study period, which will be further discussed in Section 4. Figure 3b shows the variable importance scores calculated by the RFR for the inputs listed in Table 1. The irrigation area from the AIM-HPA dataset [46] receives the highest score. This is not surprising, since a higher irrigation area is directly related to higher irrigation amounts. It is noteworthy that the AIM-HPA dataset is available for the entire HPA and does not rely on irrigated acreage records. This enables our approach to be extended to other parts of HPA, where in situ records such as WIMAS are not available. The irrigation area is correlated with remote sensing vegetation indices, thus, overshadows the vegetation indices in terms of variable importance scores (Figure 3b).
In Experiment C, we performed the same random split test for years 2019 and using the ECOSTRESS in one run and the MODIS ET in another run. We compared accuracy resulting from the use of the ECOSTRESS and MODIS ET products.

Experiment A: Random Split
Overall, the RFR achieves satisfactory performance, with an average RMSE, PB and 2 of 22.9 mm, −0.73%, and 0.82, respectively ( Figure 3). Despite interannual v bilities, a declining trend in 2 is noticeable primarily due to the decreasing variabili the irrigation amount among the grids in the study area ( Figure S2). The decrease in iability suggests that the irrigation behavior became more homogeneous. Possible ca may be water use being constrained by local water conservation programs and the ex sion of water-saving irrigation technology during the study period, which will be fur discussed in Section 4. Figure 3b shows the variable importance scores calculated by RFR for the inputs listed in Table 1. The irrigation area from the AIM-HPA dataset receives the highest score. This is not surprising, since a higher irrigation area is dir related to higher irrigation amounts. It is noteworthy that the AIM-HPA dataset is a able for the entire HPA and does not rely on irrigated acreage records. This enables approach to be extended to other parts of HPA, where in situ records such as WIMAS not available. The irrigation area is correlated with remote sensing vegetation ind thus, overshadows the vegetation indices in terms of variable importance scores (Fi 3b).
The irrigation volume of the Kansas HPA predicted by the RFR overall follows interannual variability in the in situ irrigation records from WIMAS ( Figure 4). Both th situ and estimated irrigations peaked in 2011 and 2012, which were drought years lower-than-usual growing season precipitation ( Figure 4). The RFR model tends to form better in dry years than in wet years (e.g., 2015, 2019) ( Figure S1). The RFR cont ously overpredicts the annual irrigation amount over Kansan HPA since 2016. This be due to irrigation technology and water management policy changes, as discusse Section 4.1.  The irrigation volume of the Kansas HPA predicted by the RFR overall follows the interannual variability in the in situ irrigation records from WIMAS ( Figure 4). Both the in situ and estimated irrigations peaked in 2011 and 2012, which were drought years with lower-than-usual growing season precipitation (Figure 4). The RFR model tends to perform better in dry years than in wet years (e.g., 2015, 2019) ( Figure S1). The RFR continuously overpredicts the annual irrigation amount over Kansan HPA since 2016. This may be due to irrigation technology and water management policy changes, as discussed in Section 4.1.

Experiment B: Spatial Split
The RFR model achieved a lower but still satisfactory accuracy in the Experiment B spatial split test (Figures 5 and 6) than in Experiment A (random split). Noteworthy, the annual precipitation during the study period was in the range of 431-635 mm for the training region and 584-863 mm for the test region. The RFR predicted irrigation amount for the test region is systematically higher than in situ records, suggesting an overestimation.
Nevertheless, the RFR prediction is overall following the interannual variability of the in situ records. For example, the peak irrigation amount from the RFR estimation and the in situ records both occurred in 2011 and 2012. This is the same as in Experiment A and consistent with the drought conditions in the two years.

Experiment B: Spatial Split
The RFR model achieved a lower but still satisfactory accuracy in the Experiment B spatial split test (Figures 5 and 6) than in Experiment A (random split). Noteworthy, the annual precipitation during the study period was in the range of 431-635 mm for the training region and 584-863 mm for the test region. The RFR predicted irrigation amount for the test region is systematically higher than in situ records, suggesting an overestimation. Nevertheless, the RFR prediction is overall following the interannual variability of the in situ records. For example, the peak irrigation amount from the RFR estimation and the in situ records both occurred in 2011 and 2012. This is the same as in Experiment A and consistent with the drought conditions in the two years.
We also examined the spatial distribution of the irrigation amounts of the in situ records and RFR estimates, respectively, for a moderate model performance year (Figure 6) and other years ( Figures S3 and S4). The dashed line separates the training (western) and test (eastern) portions of the study area. When trained over the western Kansas HPA, the RFR model is able to capture the spatial variability of the irrigation amounts in the eastern Kansas HPA reasonably well despite a slight overestimation.   . Annual irrigation amounts from the RFR (red triangles) and in situ records (green squares) for the test data (eastern portion of Kansas HPA) in Experiment B (spatial split). The growing season (May, June, July, August, MJJA) total precipitations for training (grey bars) and test (purple bars) portions, respectively, are also shown.  Table 2 lists the PBIAS, RMSE, and 2 metrics resulting from Experiment C (see Section 2.2). To constrain the uncertainty from random sampling in the random forest algorithm, the RFR was repeated 200 times, and the mean and standard deviations of the metrics are reported. Based on the three metrics, ECOSTRESS PT-JPL, dis-ALEXI, and MOD16 ET products yield irrigation volume estimates of similar accuracy. The average variable importance scores calculated by repeating the RFR 200 times using each ET product show that MOD16 ET receives a slightly higher importance score than ECOSTRESS ET products ( Figure 7). However, all ET products appear to be less important than the irrigation area, hydrometeorological data, and vegetation indices. It is also noteworthy that the vegetation indices receive a higher importance score in Experiment C than in Experiment A. This may be because these variables became more important for capturing spatial irrigation patterns during a shorter time period (2 years). We also examined the spatial distribution of the irrigation amounts of the in situ records and RFR estimates, respectively, for a moderate model performance year ( Figure 6) and other years ( Figures S3 and S4). The dashed line separates the training (western) and test (eastern) portions of the study area. When trained over the western Kansas HPA, the RFR model is able to capture the spatial variability of the irrigation amounts in the eastern Kansas HPA reasonably well despite a slight overestimation. Table 2 lists the PBIAS, RMSE, and R 2 metrics resulting from Experiment C (see Section 2.2). To constrain the uncertainty from random sampling in the random forest algorithm, the RFR was repeated 200 times, and the mean and standard deviations of the metrics are reported. Based on the three metrics, ECOSTRESS PT-JPL, dis-ALEXI, and MOD16 ET products yield irrigation volume estimates of similar accuracy. The average variable importance scores calculated by repeating the RFR 200 times using each ET product show that MOD16 ET receives a slightly higher importance score than ECOSTRESS ET products ( Figure 7). However, all ET products appear to be less important than the irrigation area, hydrometeorological data, and vegetation indices. It is also noteworthy that the vegetation indices receive a higher importance score in Experiment C than in Experiment A. This may be because these variables became more important for capturing spatial irrigation patterns during a shorter time period (2 years).

Groundwater Irrigation Trend
Despite the strong climatic control on the irrigation amount interannual variability (e.g., higher irrigation in drought years), irrigation water use has overall decreased in the Kansas HPA during the study period. The declining trend is primarily due to irrigation management policy and irrigation technology changes in recent years [40,72]. Since 2006, most Kansas HPA irrigation fields have adapted central pivot irrigation to central pivot with Low Energy Precise Application (LEPA). The LEPA technology can reduce wind drift-caused ET, therefore, improving the irrigation efficiency and reducing the irrigation intensity (or irrigation volume per irrigated area) required for crop growing [73]. Figure  8 shows the spatially averaged irrigation intensity during the study period for the training (western) and test (eastern) portions of the study area. It was calculated by dividing the total irrigation amount of each grid with the irrigated acreage in that grid (as opposed to the grid size, 36 km 2 , used to calculate the irrigation depth) and taking the average across the Kansas HPA. Both the irrigation amount and irrigated acreage are from the WIMAS dataset.
The declining trends in irrigation amount and intensity can also be attributed to the Local Enhanced Management Area (LEMA) groundwater conservation program. The LEMA was initially introduced in the GMD 4 Sheridan County 6 from 2012 and has been expanding within GMD 4 and to other GMDs since then [40,74,75]. Irrigators participating in the LEMA project either switched to water conservation crops or reduced the irrigation intensity, resulting in decreased irrigation amounts [73].

Groundwater Irrigation Trend
Despite the strong climatic control on the irrigation amount interannual variability (e.g., higher irrigation in drought years), irrigation water use has overall decreased in the Kansas HPA during the study period. The declining trend is primarily due to irrigation management policy and irrigation technology changes in recent years [40,72]. Since 2006, most Kansas HPA irrigation fields have adapted central pivot irrigation to central pivot with Low Energy Precise Application (LEPA). The LEPA technology can reduce wind drift-caused ET, therefore, improving the irrigation efficiency and reducing the irrigation intensity (or irrigation volume per irrigated area) required for crop growing [73]. Figure 8 shows the spatially averaged irrigation intensity during the study period for the training (western) and test (eastern) portions of the study area. It was calculated by dividing the total irrigation amount of each grid with the irrigated acreage in that grid (as opposed to the grid size, 36 km 2 , used to calculate the irrigation depth) and taking the average across the Kansas HPA. Both the irrigation amount and irrigated acreage are from the WIMAS dataset. Despite the declining irrigation amount, the study period has experienced an increase in the irrigated area, according to the AIM-HPA dataset ( Figure S5). Given the high importance score received by the irrigation area, the accuracy of the irrigation amount The declining trends in irrigation amount and intensity can also be attributed to the Local Enhanced Management Area (LEMA) groundwater conservation program. The LEMA was initially introduced in the GMD 4 Sheridan County 6 from 2012 and has been expanding within GMD 4 and to other GMDs since then [40,74,75]. Irrigators participating in the LEMA project either switched to water conservation crops or reduced the irrigation intensity, resulting in decreased irrigation amounts [73].
Despite the declining irrigation amount, the study period has experienced an increase in the irrigated area, according to the AIM-HPA dataset ( Figure S5). Given the high importance score received by the irrigation area, the accuracy of the irrigation amount estimation is sensitive to the accuracy of the AIM-HPA. A higher irrigation area in later years likely contributed to the overestimation of the irrigation amounts after 2015. While AIM-HPA has been validated using the USDA National Agricultural Statistics Service (NASS) records at county level and achieved an R 2 = 0.86 to 0.93 for counties within the central High Plains region, the map tends to have commission (false positive) errors in humid climate areas with rainfed crops [46]. Therefore, it might overestimate the irrigation area within our study area.
Declining trends were found in the GI and NDVI during the study period ( Figure S6). However, NDWI exhibited an increasing trend, seemingly contradictory to the expectation that a reduced irrigation intensity would result in a lower NDWI. Given the opposite trends and low importance scores received by these indices, they are likely overshadowed by the irrigation areas that have increased over years. Therefore, the declining trend in the irrigation amounts caused by irrigation policy and technology changes is not adequately captured by the input variables used in this study. In Experiments A and B, the RFR model was trained and tested using data throughout the study period. The trained RFR model sought an overall goodness-of-fit; thus, it learned the "average" irrigation behavior. Therefore, the RFR slightly underestimated the annual total irrigation amount in earlier years, while it overestimated it in later years of the study period ( Figure 4). Including irrigation technology and policy information in the model inputs (e.g., LEMA regulated the cultivated area fraction of each grid) might improve the estimation accuracy, although limiting the transferability of our approach to other regions without such information.

Spatial Transferability
In addition to the temporal shift in irrigation technology and practices, the study area exhibits a spatial variability of climatic conditions, crop types, and irrigation practices that are typical for HPA. Experiment B is a challenging case where the training and test datasets differ in climate, crop types, and irrigation practices to examine the spatial transferability of our approach. The western part (GMDs 1, 3, and 4) is dryer; thus, it needs an overall higher irrigation amount than the eastern part (GMDs 2 and 5). Additionally, the training and test areas differ in fractions of the major crops. In particular, the test region has substantially more soybeans and alfalfa than the training region ( Figure 9a). Overall, the major crops had a higher irrigation intensity in the dryer training region (Figure 9b). As a result, the training and test data in Experiment B have different input distributions. Predicting at a data point different from the training data (out-of-distribution) is a well-known challenge for machine learning techniques [76]. By using a variety of input variables, such as remote sensingderived data, our approach partially overcomes the out-of-distribution challenge. The RFR is shown to be relatively skillful in capturing the inter-annual and spatial variability in irrigation amounts in Experiment B. This suggests that our approach of estimating irrigation amounts using remote sensing and hydroclimatic data can be potentially extended to other parts of the HPA.
Predicting at a data point different from the training data (out-of-distribution) is a wellknown challenge for machine learning techniques [76]. By using a variety of input variables, such as remote sensing-derived data, our approach partially overcomes the out-ofdistribution challenge. The RFR is shown to be relatively skillful in capturing the interannual and spatial variability in irrigation amounts in Experiment B. This suggests that our approach of estimating irrigation amounts using remote sensing and hydroclimatic data can be potentially extended to other parts of the HPA.

ECOSTRESS ET Utility for Irrigation Amount Estimation
Despite their higher spatial resolution and accuracy, the ECOSTRESS ET products (dis-ALEXI daily and PT-JPL instantaneous) generated irrigation estimations similarly accurate as the results from MODIS ET. This is likely due to the low variable importance scores received by the ET products. Since ET is highly correlated to VPD, its importance as an explanatory variable for irrigation amounts may have been overshadowed by VPD. Uncertainty associated with remote sensing ET products further reduces their information content for estimating irrigation amounts. Another reason for the similar accuracy of using ECOSTRESS and MODIS ET may be the low temporal and spatial coverage of the ECOSTRESS images. As a result, important irrigation information during the growing season may be missed. Recent work on fusing the ECOSTRESS ET with spatially coarser yet temporally continuous ET products, such as MODIS [77], to generate a high spatial and temporal resolution ET holds great potential to further improve irrigation amount estimations.

Conclusions
In this research, we develop a random forest regression model to fuse multiple remote sensing products, hydrometeorological information, and in situ annual irrigation records in the Kansas portion of the HPA to estimate annual irrigation amounts at a 6 km spatial resolution. Specifically, we use remote sensing ET products, vegetation indices (GI, NDVI, NDWI), crop types, climate data (precipitation, temperature, water vapor deficit), and a remote sensing-derived irrigation map dataset as input variables and train a random forest model using in situ irrigation records. In all three experiments, the RFR model

ECOSTRESS ET Utility for Irrigation Amount Estimation
Despite their higher spatial resolution and accuracy, the ECOSTRESS ET products (dis-ALEXI daily and PT-JPL instantaneous) generated irrigation estimations similarly accurate as the results from MODIS ET. This is likely due to the low variable importance scores received by the ET products. Since ET is highly correlated to VPD, its importance as an explanatory variable for irrigation amounts may have been overshadowed by VPD. Uncertainty associated with remote sensing ET products further reduces their information content for estimating irrigation amounts. Another reason for the similar accuracy of using ECOSTRESS and MODIS ET may be the low temporal and spatial coverage of the ECOSTRESS images. As a result, important irrigation information during the growing season may be missed. Recent work on fusing the ECOSTRESS ET with spatially coarser yet temporally continuous ET products, such as MODIS [77], to generate a high spatial and temporal resolution ET holds great potential to further improve irrigation amount estimations.

Conclusions
In this research, we develop a random forest regression model to fuse multiple remote sensing products, hydrometeorological information, and in situ annual irrigation records in the Kansas portion of the HPA to estimate annual irrigation amounts at a 6 km spatial resolution. Specifically, we use remote sensing ET products, vegetation indices (GI, NDVI, NDWI), crop types, climate data (precipitation, temperature, water vapor deficit), and a remote sensing-derived irrigation map dataset as input variables and train a random forest model using in situ irrigation records. In all three experiments, the RFR model achieves a satisfactory performance and generates annual irrigation estimates in good agreement with in situ records. The spatial split test with spatially different climates and irrigation patterns between the training and test regions, which appears to be the most challenging, indicates that the RFR can still capture the irrigation patterns relatively well, thus, suggesting the model's potential for spatial transferability. Finally, the irrigation amounts estimated at the 6 km resolution using two ECOSTRESS ET products (dis-ALEXI daily and PT-JPL instantaneous) have a similar accuracy compared with the MODIS ET results. This is likely due to the low importance scores received by the ET products and their correlation with the VPD, as well as the ECOSTRESS images may miss key irrigation periods due to their low spatial and temporal coverage. However, the high spatial resolution of the ECOSTRESS ET is vital for estimating irrigation amounts at the field scale. In addition, recent advances in developing a fused ET product that has high spatial and temporal resolutions could further improve the accuracy of irrigation amount estimations.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/rs14133004/s1, Figure S1: Experiment A (random split) test R 2 varies with the growing season (MJJA) precipitation amounts for 2006-2020; A statistically significant (p = 0.02) linear declining trend is also shown; Figure S2: Annual variance of in situ irrigation depth (mm 2 ) of the test data in Experiment A (random split). Also shown is a statistically significant declining trend (p = 0.03); Figure S3: (a) In situ and (b) RFR estimated irrigation amounts in Experiment B (spatial split) for the low performance year 2020 with an R 2 = 0.37. The red dash line divides the training (western) and test (eastern) regions; Figure S4: (a) In situ and (b) RFR estimated irrigation amounts in Experiment B (spatial split) for the good performance year 2006 with an R 2 = 0.77. The red dash line divides the training (western) and test (eastern) regions; Figure S5: Annual total irrigation acreage within the study area from in situ records (WIMAS, purple) and a remote sensing derived irrigation map (AIM-HPA, grey); Figure S6: Normalized annual GI (green, triangles), NDVI (red, squares), and annual NDWI (blue, circle) averaged across the study area during 2006-2020. Linear regression lines of GI, NDVI, and NDWI are also shown with p-values.