Sensitivity of Potential Groundwater Recharge to Projected Climate Change Scenarios: A Site-Specific Study in the Nebraska Sand Hills, USA

: Assessing the relationship between climate forcings and groundwater recharge (GR) rates in semi-arid regions is critical for water resources management. This study presents the impact of climate forecasts on GR within a probabilistic framework in a site-speciﬁc study in the Nebraska Sand Hills (NSH), the largest stabilized sand dune region in the USA containing the greatest recharge rates within the High Plains Aquifer. A total of 19 downscaled climate projections were used to evaluate the impact of precipitation and reference evapotranspiration on GR rates simulated by using HYDRUS 1-D. The analysis of the decadal aridity index (AI) indicates that climate class will likely remain similar to the historic average in the RCP2.6, 4.5, and 6.0 emission scenarios but AI will likely decrease signiﬁcantly under the worst-case emission scenario (RCP8.5). However, GR rates will likely decrease in all of the four emission scenarios. The results show that GR generally decreases by ~25% under the business-as-usual scenario and by nearly 50% in the worst-case scenario. Moreover, the most likely GR values are presented with respect to probabilities in AI and the relationship between annual-average precipitation and GR rate were developed in both historic and projected scenarios. Finally, to present results at sub-annual time resolution, three representative climate projections (dry, mean and wet scenarios) were selected from the statistical distribution of cumulative GR. In the dry scenario, the excessive evapotranspiration demand in the spring and precipitation deﬁcit in the summer can cause the occurrence of wilting points and plant withering due to excessive root-water-stress. This may pose signiﬁcant threats to the survival of the native grassland ecology in the NSH and potentially lead to desertiﬁcation processes if climate change is not properly addressed.


Introduction
Socioeconomic drivers have historically determined and will continue to dictate the ever-increasing rate of groundwater depletion around the world. Water consumption in most countries has considerably increased over recent decades due to population and economic growth, with irrigation water use being the leading consumer. In addition to overexploitation, groundwater recharge (GR) regimes are affected by climate change and variability, which threatens environmental sustainability and the livelihoods of billions around the world. The dire water scarcity issues make it imperative to assess the consequences of climate change on water resources at global, continental, and local scales with particular emphasis on GR rates that feed shallow aquifers [1][2][3][4]. Groundwater levels of many aquifers have been showing a decreasing trend over the last few decades due to excessive groundwater extraction for irrigation that surpasses GR and replenishing rates [5,6]. The vulnerability of groundwater resources emphasizes the need for stakeholders and decision-makers to have reliable information regarding the relationship between climate and GR rates, particularly in semi-arid regions where water scarcity issues are well documented [6].
The assessment of climate change can vary in the type of impact, scale, and severity. The negative consequences of climate change in some settings can prove to be pronounced and include severe droughts [7], sea-level rise [8], wildfires [9], extreme heat waves [10], floods and inundation [11], reduction of biodiversity [12], and many others.
The impact of climate change on the water budget in semi-arid regions primarily occurs via changes to potential evapotranspiration (ET) rates and precipitation (P) regime [13]. The rise in temperatures can cause increases in potential ET [14]. Similarly, climate change has been reported to impact the amount and seasonality of P regime [15] which can alter surface runoff and streamflow [16][17][18][19], soil moisture [20], and GR rates [21][22][23].
Climate change can impact GR rates by altering the availability of water through changes in P and temperature, as well as in the composition of vegetation [24]. While on the one hand some studies have found increases in GR rates in areas experiencing permafrost thaws due to temperature rise and increase in P [25,26], on the other hand, several studies have estimated decreased recharge as a result of reduced P or greater temperatures [23,27,28]. Further, Bouraoui et al. (1999) [29] reported considerable reductions in GR near Grenoble, France, predominantly due to increases in evaporation during the recharge season. The extent to which GR is affected by climate warming will also depend on the local soils and vegetation in addition to climate [30].
The objective of this paper is to examine the relationship between historic and projected future climate forcings (in the form of precipitation, reference evapotranspiration, and aridity index) and the GR rates in a site-specific study in the Nebraska Sand Hills (NSH). This groundwater-dependent grassland ecosystem represents the largest stabilized sand dune formation in the western hemisphere. One of the largest grass-stabilized dune regions in the world, the NSH provides the greatest recharge rates within the High Plains Aquifer under semi-arid climate conditions. This paper is a follow-up of four studies in the NSH region [31][32][33][34]. The remainder of this paper is organized as follows: Section 2 presents a brief description of the study area, Nebraska National Forest (NNF), located in the NSH and modeling methodology based on HYDRUS-1D. Section 3 discusses the analysis of climate change impacts on GR rates by comparing the historic and projected water balance at a decadal, annual and sub-annual resolution within a probabilistic framework. Section 4 provides a concluding summary of the analysis presented herein.

Study Area and Modeling Approach
The study site is located in a century-old experimental forest (Nebraska National Forest, near the town of Halsey, Nebraska, USA) surrounded by mixed-grass prairie in the south-central part of the NSH ( Figure 1). The NSH landscape is comprised mainly of eolian sand dunes that were deposited as recently as a few thousand years ago [35]. The soil is approximately 92%-97% sand [36], which contributes to the greatest recharge rates in the High Plains Aquifer [37]. The native vegetation of the NSH region consists of mixed-prairie grassland including little bluestem (Schizachyrium scoparium), switchgrass (Panicum virgatum), sand dropseed (Sporobolus cryptandrus), and Kentucky bluegrass (Poa pratensis), and is suitable to the historic land uses of ranching and cattle grazing. A weather station (Longitude −100 • 15 00 ', Latitude 41 • 9 00 ', elevation 824 m) recorded daily values of precipitation, wind speed, air temperature, relative humidity and solar radiation from 1990 to 2012. The climate is semi-arid continental with annual precipitation (P) ranging from 36.0 to 82.0 cm yr −1 , and with reference crop evapotranspiration (ET 0 , referred to a well-watered hypothetical reference alfalfa-crop) ranging from 128.5 to 166.7 cm yr −1 . We remind that with the availability of all meteorological data, we used the well-known standard Penmann-Monteith equation for computing ET 0 as recommended by the Food and Agriculture Organization (FAO), in Paper no. 56 [38]. Overland flow is negligible due to the high infiltration capacity of the sandy soils. Other details are described by Adane and Gates (2015) [39]. The native grassland plot (Longitude −100 • 21 16 ', Latitude 41 • 50 41 ', elevation 860 m), denoted with the acronym G, represents the natural grassland ecosystem in the NSH [32,33].  [38]. Overland flow is negligible due to the high infiltration capacity of the sandy soils. Other details are described by Adane and Gates (2015) [39]. The native grassland plot (Longitude −100°21′16′', Latitude 41°50′41′', elevation 860 m), denoted with the acronym G, represents the natural grassland ecosystem in the NSH [32,33]. The effect of climate change on potential GR rates was assessed using a scenario-based approach by considering historic  and plausible climate projections (2000-2100). Both historic and projected climate data contained daily values of P, and minimum and maximum temperatures obtained from the High Plains Regional Climate Center (HPRCC). In this case, la ack of data for wind speed, relative humidity, and solar radiation made it a necessity to use the temperature-based Hargreaves equation [40] to estimate reference crop evapotranspiration (ET0). The formula only requires minimum and maximum temperature data while the extraterrestrial radiation term is estimated by using study site latitude and day of the year. Nonetheless, the performance of the Hargreaves equation to estimate ET0 has considerable limitations given the limited input dataset [41]. For this reason, the performance of Hargreaves formula was tested against the Penman-Monteith equation, by using the meteorological data at Halsey weather station. In doing so, ET0-values based on solely temperature data were compared with those based on a complete set of meteorological data from 1990 till 2012. The results indicate that the Hargreaves equation underestimates annual average ET0 (119.2 cm yr −1 ) compared to ET0 values estimated by Penman-Monteith equation (145.5 cm yr −1 ) [42]. The root mean square error (RMSE) and R 2 -value calculated on the daily ET0-values are 0.11 cm d −1 and 0.88, respectively. Archibald and Walters (2014) [43] found similar and acceptable performance at Mead station in eastern Nebraska (Latitude 41°17′00′'). A total of 19 climate models were used in this study (Table 1). The effect of climate change on potential GR rates was assessed using a scenario-based approach by considering historic  and plausible climate projections (2000-2100). Both historic and projected climate data contained daily values of P, and minimum and maximum temperatures obtained from the High Plains Regional Climate Center (HPRCC). In this case, la ack of data for wind speed, relative humidity, and solar radiation made it a necessity to use the temperature-based Hargreaves equation [40] to estimate reference crop evapotranspiration (ET 0 ). The formula only requires minimum and maximum temperature data while the extraterrestrial radiation term is estimated by using study site latitude and day of the year. Nonetheless, the performance of the Hargreaves equation to estimate ET 0 has considerable limitations given the limited input dataset [41]. For this reason, the performance of Hargreaves formula was tested against the Penman-Monteith equation, by using the meteorological data at Halsey weather station. In doing so, ET 0 -values based on solely temperature data were compared with those based on a complete set of meteorological data from 1990 till 2012. The results indicate that the Hargreaves equation underestimates annual average ET 0 (119.2 cm yr −1 ) compared to ET 0 values estimated by Penman-Monteith equation (145.5 cm yr −1 ) [42]. The root mean square error (RMSE) and R 2 -value calculated on the daily ET 0 -values are 0.11 cm d −1 and 0.88, respectively.
Archibald and Walters (2014) [43] found similar and acceptable performance at Mead station in eastern Nebraska (Latitude 41 • 17 00 '). A total of 19 climate models were used in this study (Table 1). Each climate model considered four plausible scenarios based on future greenhouse gas concentration trajectories (Representative Concentration Pathways, RCPs = 2.6, 4.5, 6.0, and 8.5). The four RCPs represent possible ranges of radiative forcing values in the year 2100 relative to pre-industrial values (+2.6, +4.5, +6.0, and +8.5 W/m 2 , respectively). RCP 2.6 assumes that global annual greenhouse gas emissions peak between 2010 and 2020 and decline substantially thereafter (optimistic scenario). Emissions in RCP 4.5 and RCP 6.0 peak around 2040 and 2080 then decline, respectively while the RCP 8.5 scenarios assume emissions continue to rise throughout the 21st Century (worst-case scenario). A total of 76 climate projections (four emission scenarios over 19 climate models) were expected, but only 65 were available because some climate models did not provide time series of daily P, or minimum and maximum temperatures under RCPs 4.5 and 6.0 emission scenarios. The aridity index, AI is described as the annual P over the annual ET 0 (AI = P/ET 0 ) and is commonly used for climate classification [44][45][46][47]. This climate indicator was widely used for evaluating the effects of climate change on overall water balance [48][49][50].
The water balance in the soil-plant-atmosphere system was numerically evaluated using HYDRUS 1-D [51]. Numerical simulations at daily temporal resolution utilized available P data and the calculated ET 0 from the Hargreaves equation.
The model setup in HYDRUS 1-D requires vegetation parameters, which were determined by direct measurements or were obtained from tabulated values in the scientific literature [33]. It is worth noting that considering net P (P net ) instead of total precipitation accounts for precipitation that was intercepted by foliage cover and does not reach the soil surface. The leaf area index (LAI) determines the amount of precipitation interception that is subtracted from P in order to obtain P net falling on the soil surface [52].
The crop coefficient (K c ) of grass converts ET 0 into crop potential evapotranspiration (ET p ). The LAI is used for partitioning ET p into potential evaporation (E p ) and potential transpiration (T p ) [33]. Therefore, P net and E p represent the system-dependent time-variable water fluxes of the upper boundary condition whilst T p determines the potential root water uptake and depends on maximum root depth (z r ) and root distribution (assumed to be linear along the soil profile by varying from maximum at the soil surface to minimum at z r ). Both E p and T p are reduced by water stress conditions into actual evaporation (E a ) and actual transpiration (T a ) depending on soil moisture storage status. A set of five unknown soil hydraulic parameters (SHPs) belonging to the water retention curve and the hydraulic conductivity curve [53], denoting residual soil water content, saturated soil water content, two water retention shape parameters, and saturated hydraulic conductivity, were calibrated and validated in the study presented by Adane et al. (2018) [33]. The depth of the soil profile was assumed to be 220 cm following Adane et al. (2018) [33], and the lower boundary condition is set to free drainage representing potential GR, as the depth to the water table is at least 4 m in the NNF [31]. The assessment of the lag time will provide an estimate of the actual GR to the water table, rather than just draining through the root zone [31]. Initial conditions were set in terms of pressure head values along the soil profile to impose hydraulic equilibrium.
The total number of numerical simulations in HYDRUS-1D included: Actual transpiration (T a ), actual evaporation (E a ), and GR rates represent the model outputs that were stored at daily time resolution and subsequently aggregated to monthly and annual sums.
In this study, the Student's t-distribution is used as an alternative to the normal distribution when sample sizes are small to estimate the expected value of P, ET 0 , ET p , AI, E a , T a and GR. The Student's t-distribution is used to construct the 95% confidence interval (α = 0.05) around the sample mean x in each emission scenario to infer the expected value (µ) of a continuously distributed population: where n is the number of observations (degrees of freedom are therefore defined as n − 1), and s is the sample standard deviation, and t α,n−1 is student's t-value at 95% confidence interval and n − 1 degrees of freedom.

Results and Discussion
This Section is divided into three sub-sections. The first sub-section presents the comparison between historic and future decadal-average GR rates in response to varying climate forcings. The second part presents the annual average GR rates. Finally, in the third sub-section, three representative scenarios were selected out of 65 available climate projections to compare results at sub-annual time resolution.

Multi-Decadal-Average Results of Historical and Projected ET a and GR
Multi-decadal-average values (referring to the average of annual values of 2000-2100) of climate forcings (P, ET 0 , ET p ), aridity index (AI) and model output (ET a and GR) for historic and climate projections (16 projections in RCP2.6, 18 projections in RCP4.5, 12 projections in RCP6.0, 19 projections in RCP8.5) are listed in Table A1 (in the Appendix A). Given that the climate models were relatively few, the observations were assumed to be distributed normally (P, ET 0 , AI, ET p , ET a ) or lognormally (GR).
The historical average values of all 50 years (1950-2000) are presented in the first row of Table A1. The historical aridity index value was AI = 0.49 corresponding to the semi-arid class (0.2 < AI ≤ 0.5) based on aridity index classification reported in Spinoni et al. (2015) [54]. The simulated historical data indicate that ET p and ET a are 69% and 40% of the ET 0 , respectively. Considering negligible losses to precipitation interception and overland flow, baseline GR was approximately 16% of the historic P. This result agrees with previous GR estimates in the Nebraska Sand Hills [33,39,55,56] that on average estimated approximately 10 cm yr −1 under native grassland conditions. However, spatial differences may be important across the Sand Hills as estimated net GR rates based on remote-sensing and other ancillary climate data. For instance, data presented by Szilagyi et al. (2011) [57] for the period 2000-2009 showed that GR rates varied spatially from approximately −0.3 to 12.5 cm yr −1 where the negative values indicate that in some parts of the NNF annual ET was greater than P. Detailed comparisons of historical GR rates in the Sand Hills in general and NNF, in particular, are available in Adane and Gates (2015) [39].
In all emission scenarios evaluated, the projected rise in temperature generated a corresponding rise in ET 0 , resulting in projected multi-decadal-averages of 123.8 ± 1.8 cm yr −1 , 126.5 ± 1.9 cm yr −1 , 127.3 ± 2.3 cm yr −1 , and 149.8 ± 8.7 cm yr −1 , for RCP2.6, RCP4.5, RCP6.0 and RCP8.5, respectively (see Equation (1)). The mean values correspond to an 8.5%, 10.9%, 11.6%, and 31.4% increase, respectively, when compared to the historic value (114.1 cm yr −1 ). It is also worth noting that the uncertainty (given by the confidence interval around the mean quantified by Equation (1)) in ET 0 rise also increased with increasing CO 2 concentration projected, with the greatest uncertainty observed under the worst-case climate scenario (RCP8.5). Further, P is characterized by a slight increase (about 1% for all four emission scenarios) when compared to the decadal-average historic value (55.7 cm yr −1 ). This result agrees with the conclusion by Lauffenburger et al. (2018) [58] that the projected change in P between 1990 and 2050 in western Nebraska (Northern High Plains) is a statistically non-significant increase.
Due to the substantial differences in the rates of change between ET 0 and P, the use of a synthetic climate measure such as the aridity index (AI) is effective to evaluate recharge rates between the historic and climate change projections. Figure 2a shows the box plots of the projected multi-decadal-average AI-values obtained from the 19 climate models for each emission scenario. The AI values resulting from the analysis were 0.47 ± 0.014, 0.45 ± 0.015, 0.46 ± 0.021, and 0.40 ± 0.029, for RCP2.6, RCP4.5, RCP6.0 and RCP8.5, respectively (see Equation (1)). These results suggest that the region is projected to become more arid (with decreasing AI-values) compared to the historical baseline (AI = 0.49) but still remaining in the same semi-arid class [59]. The greatest change in aridity is likely to occur only under the RCP8.5 climate projection scenario.
The ANOVA test at 5% significance level rejected the hypothesis that the mean AI values of the four groups were equal to each other given that F-value = 12.75 was larger than the critical F-value. However, if the AI-values of the RCP8.5 were removed, the ANOVA test would accept the hypothesis that the mean values of the RCP2.6, RCP4.5 and RCP6.0 group were significantly equal. This statistical analysis suggests that the change in aridity index under the climate change scenarios was only significant under the worst-case scenario.
The analysis also indicates that all climate scenarios caused a decline in GR rates compared to the historical baseline (9.02 cm yr −1 ). Projections of climate change induced a general decrease in the simulated potential GR, resulting in future period multi-decadal-averages of 7.15 ± 0.60 cm yr −1 , 6.87 ± 0.71 cm yr −1 , 6.95 ± 0.74 cm yr −1 , and 5.51 ± 0.78 cm yr −1 , for RCP2.6, RCP4.5, RCP6.0 and RCP8.5, respectively (see Equation (1)). The mean values represent decreases ranging from 21% in RCP2.6 to 39% in RCP8.5 in terms of percentage values (see box plots in Figure 2b). Nevertheless, box plots denote high uncertainty around the median value, with a coefficient of variation spanning from 18% to 30% for RCP2.6 and RCP8.5, respectively. The decrease in potential GR rates is due to a substantial increase in ET losses. This climatic setting also poses questions on probabilities related to plant survival, covered in Sub-Section 3.3.

Annual-Average Results of Historical and Projected AI, ETa, and GR
The model simulations comprise 50 years of historical and 100 years of 65 climate projections of climate forcings (P, ET0, ETp), aridity index (AI), and model output ETa and GR. The resulting annual-averages covering the projected 100-year simulation period are presented within a probabilistic framework shown in empirical histograms of AI values and corresponding GR ( Figure  3). The normal and lognormal distributions (red lines) were fit on annual AI values and GR rates (blue histograms).
The probabilities for each climate class and corresponding median potential GR rates are reported in Table 2. Approximately 69% of the total aridity index values fall in the semi-arid climate class and are associated with the median potential GR rate (GR = 5.19 cm) that is significantly lower than the historic GR (9.02 cm yr −1 , see Table A1). The probabilities of the climate shifting from semiarid to dry (0.5 < AI ≤ 0.65) are approximately 26% and would produce a median potential GR of 9.64 cm yr −1 , which is slightly greater (4.6%) than the historic value (9.02 cm yr −1 ). The combined probabilities attributed to extreme conditions such as arid, sub-humid, and humid climates were very low (~5%). The median-annual GR rates for the extreme climate conditions were far from the baseline rates.  [58] reported that median annual recharge for the low (+1.0 • C, 478 ppm CO 2 ) and high (+2.4 • C, 567 ppm CO 2 ) global warming scenarios in western Nebraska rangeland are projected to decrease by about 9% and 94% relative to historical values, respectively.
The decrease in potential GR rates is due to a substantial increase in ET losses. This climatic setting also poses questions on probabilities related to plant survival, covered in Section 3.3.   The coarse sandy soils in the Nebraska Sand Hills are not conducive to overland flow, thus the water balance of the soil profile is likely controlled by vertical hydrological processes, namely infiltration and ET. Moreover, precipitation is one of the main drivers of GR under semi-arid climate conditions [60]. The relationship between annual P and GR was quantified through regression lines and the corresponding coefficient of determination, namely the R-squared (R 2 ), for evaluating its goodness of fit for both the historic and projected data group. The data are presented in Table 3 and visualized through a scatter plot in Figure 4, where the circles are colored according to the aridity index values.  The probabilities for each climate class and corresponding median potential GR rates are reported in Table 2. Approximately 69% of the total aridity index values fall in the semi-arid climate class and are associated with the median potential GR rate (GR = 5.19 cm) that is significantly lower than the historic GR (9.02 cm yr −1 , see Table A1). The probabilities of the climate shifting from semi-arid to dry (0.5 < AI ≤ 0.65) are approximately 26% and would produce a median potential GR of 9.64 cm yr −1 , which is slightly greater (4.6%) than the historic value (9.02 cm yr −1 ). The combined probabilities attributed to extreme conditions such as arid, sub-humid, and humid climates were very low (~5%). The median-annual GR rates for the extreme climate conditions were far from the baseline rates. The coarse sandy soils in the Nebraska Sand Hills are not conducive to overland flow, thus the water balance of the soil profile is likely controlled by vertical hydrological processes, namely infiltration and ET. Moreover, precipitation is one of the main drivers of GR under semi-arid climate conditions [60]. The relationship between annual P and GR was quantified through regression lines and the corresponding coefficient of determination, namely the R-squared (R 2 ), for evaluating its goodness of fit for both the historic and projected data group. The data are presented in Table 3 and visualized through a scatter plot in Figure 4, where the circles are colored according to the aridity index values. The coarse sandy soils in the Nebraska Sand Hills are not conducive to overland flow, thus the water balance of the soil profile is likely controlled by vertical hydrological processes, namely infiltration and ET. Moreover, precipitation is one of the main drivers of GR under semi-arid climate conditions [60]. The relationship between annual P and GR was quantified through regression lines and the corresponding coefficient of determination, namely the R-squared (R 2 ), for evaluating its goodness of fit for both the historic and projected data group. The data are presented in Table 3 and visualized through a scatter plot in Figure 4, where the circles are colored according to the aridity index values.  Even though the regressions are clearly affected by heteroscedasticity (with higher spread around the regression line for higher AI-values), the results presented in Figure 4 and Table 3 represent a good opportunity to roughly estimate potential GR rates from projected P with annual temporal resolution. Given the R 2 -values listed in Table 3, P is able to explain about 30%-40% of GR temporal variability. It is apparent that all projected regression lines lie below the historic line (dotted black line), which means that GR rates will likely decrease, especially under the worst-case scenario (RCP8.5, dashed red line), as can be observed from the slopes of the regression lines reported in Table 3.

Selection of Three Representative Climate Projections for Assessing the Relation between Precipitation and Groundwater Recharge Rates at Sub-Annual Time Scale
Detailed results can be explored at sub-annual resolution. Nonetheless, it was not practical to analyze simulations pertaining to all climate projections. Therefore, for the purpose of encompassing the full spectrum of most probable future climate projections, three representative climate scenarios were selected following the procedure recommended by Rossman et al. (2018) [34]. The temporal dynamics of all cumulative GR rates, colored according to RCP, are shown in Figure 5a. The corresponding frequency distribution of final cumulative GR rates recorded for the year 2100 with sample mean, x, and standard deviation, s, is displayed in Figure 5b. The dry, mean, and wet GR projections selected based on this analysis were CSIRO (RCP8.5), GFDL_M1 (RCP2.6), and MPI_MR (RCP4.5), respectively (Table A1), based on the closest GR rates identified using the x − s, x and x + s distribution criteria, respectively (Figure 5b). Figure 5a. The corresponding frequency distribution of final cumulative GR rates recorded for the year 2100 with sample mean, ̅ , and standard deviation, s, is displayed in Figure 5b. The dry, mean, and wet GR projections selected based on this analysis were CSIRO (RCP8.5), GFDL_M1 (RCP2.6), and MPI_MR (RCP4.5), respectively (Table A1), based on the closest GR rates identified using the ̅ , ̅ and distribution criteria, respectively (Figure 5b). The annual time series of ET a , P, ET 0 , and potential GR rates of the three representative scenarios are shown in Figure 6. While the projected multi-decadal-averages of P are 54.0 cm yr −1 , 59.4 cm yr −1 , 58.9 cm yr −1 , in dry, mean, wet scenarios, respectively, the three ET 0 time series are characterized by time-variant increasing rates as predicted in other semi-arid environments [61]. Thus, the corresponding AI time series (not shown) suggests that the climate might suffer from an increasing aridity even though the progressive shift does not change the projected climate class from the historical designation (semi-arid). The results of the analyses are consistent with Oglesby et al. (2015) [62] that suggested that while the form and seasonality of P will change in Nebraska, the annual amount of P may not vary substantially. Similarly, the report forecasts that Nebraska will likely experience temperature increases ranging from 4 to 5 • F under low emission scenario and 8 to 9 • F in the high emission scenario. Both emission scenarios would result in greater aridity index given that P would remain relatively stable. from the historical designation (semi-arid). The results of the analyses are consistent with Oglesby et al. (2015) [62] that suggested that while the form and seasonality of P will change in Nebraska, the annual amount of P may not vary substantially. Similarly, the report forecasts that Nebraska will likely experience temperature increases ranging from 4 to 5 °F under low emission scenario and 8 to 9 °F in the high emission scenario. Both emission scenarios would result in greater aridity index given that P would remain relatively stable.  [34] also found an increase in projected GR under both the median and wet projections at a rate of 0.27 cm yr −1 (+5%) and 2.21 cm yr −1 (+30%), and a decrease under the dry projection by 1.50 cm yr −1 (-40%) averaged over 2010-2099 relative to the baseline averaged over 2000-2009. However, it is worth noting the differences in the baseline historical recharge rates used between this analysis (9.02 cm yr −1 ) for NNF and Rossman et al. 2018 [34] (7.0 cm yr −1 ) for NSH. The difference in projected recharge rates may also have occurred due to simulation scales (spatial and temporal) and differences in the modeling approach. Spatial resolution differences can create discrepancies in the historical averages of climate forcings and the resulting parameters such as ETa and GR. While the focus of this study is on the NNF site, Crosbie et al. (2013) [56] and Rossman et al. (2018) [34] provide averaged values over Northern High Plains and the NSH, and make use of WAVES, a one-dimensional point-scale model that simulates the fluxes of mass and energy between the atmosphere, vegetation, and soil Annual-average groundwater recharge in the wet, mean, and dry projections are 8.40 cm yr −1 , 6.74 cm yr −1 , and 3.94 cm yr −1 , representing −7%, −25% and −56% changes relative to the historical baseline (9.02 cm yr −1 ). Crosbie et al. (2013) [56] estimated a baseline historical recharge of 7.8 cm yr −1 for the Northern High Plains. The selected three future climate scenarios, namely dry, median, and wet correspond to 10%, 50%, and 90% exceedances for low (+1.0 • C, 478 ppm CO 2 ), medium (+1.7 • C, 523 ppm CO 2 ), and high (+2.4 • C, 567 ppm CO 2 ) global warming scenarios. The study reported that GR was projected to be 5.93 cm yr −1 (−24%), 8.42 cm yr −1 (+8%), and 10.30 cm yr −1 (+32%) under the dry, median, wet scenarios, respectively. Rossman et al. (2018) [34] also found an increase in projected GR under both the median and wet projections at a rate of 0.27 cm yr −1 (+5%) and 2.21 cm yr −1 (+30%), and a decrease under the dry projection by 1.50 cm yr −1 (-40%) averaged over 2010-2099 relative to the baseline averaged over 2000-2009. However, it is worth noting the differences in the baseline historical recharge rates used between this analysis (9.02 cm yr −1 ) for NNF and Rossman et al. 2018 [34] (7.0 cm yr −1 ) for NSH. The difference in projected recharge rates may also have occurred due to simulation scales (spatial and temporal) and differences in the modeling approach. Spatial resolution differences can create discrepancies in the historical averages of climate forcings and the resulting parameters such as ET a and GR. While the focus of this study is on the NNF site, Crosbie et al. (2013) [56] and Rossman et al. (2018) [34] provide averaged values over Northern High Plains and the NSH, and make use of WAVES, a one-dimensional point-scale model that simulates the fluxes of mass and energy between the atmosphere, vegetation, and soil systems [63] and the variable infiltration capacity (VIC) macroscale hydrology model [64], respectively. Further, differences can also arise from the increased number of climate models used in this analysis, as opposed to the number of climate models used in Crosbie et al. (2013) [56] and Rossman et al. (2018) [34]. Despite these differences, GR rates are nearly equivalent to those reported in Crosbie et al. (2013) [56] when only looking at results derived from the same climate models. Discrepancies will continue to exist when comparing GR rates derived from averaging all climate projections over the entire NHP and NSH spatial extents. To describe and test the entirety of reasons for differences among previous studies was not an objective of this investigation, but it could be a focus for future research. In semi-arid regions, the impact of precipitation on GR is more important when analyzing the seasonal distribution rather than the annual sums [65]. The distribution of monthly sums of precipitation, actual ET, and potential GR highlights interesting aspects. Historic precipitation, actual ET and potential GR follow similar temporal trends with maxima in May (P = 9.07 cm), June (ET a = 9.82 cm), and June (GR = 1.30 cm). Minima occur in winter months likely related to negligible precipitation, low temperatures, and freezing conditions. The projected monthly distributions under the mean and wet scenarios roughly reproduce similar trends for P, while the dry scenario generates precipitation lower than the historic baseline during the summer seasons and higher than the historic baseline during the winter seasons. The dry climate projection indicates that summer months produce 17.4 cm yr −1 of P, which is a 25% reduction in summer month P compared to the historic value (23.1 cm yr −1 ). This significant water deficit is reflected in a 60% loss in GR rates in the summer, which is the major season for groundwater replenishment in this region. Spring is the next most important season for GR. However, GR was also substantially reduced in April and May months due to the projected excessive actual ET demand (+25%) in the dry scenario. Outside of the summer and spring months of the year, the values for the water balance components observed under the historic time period and the dry scenario projection were similar.
The water deficit detected by using monthly sums (Figure 7) poses the risk of threatening plant survival in this site. In this study, we applied time-invariant parameters for root water uptake as we assumed no change in grassland cover in HYDRUS-1D [17]. Prolonged root water stress is known to cause plant wilting and subsequent withering. Increasing chances of root water stress induced by a change in climate will have consequences to the native grassland ecology of the Nebraska Sand Hills. In this study, daily values of root water stress were explored in the historical data and the mean, wet, and dry scenario projections. Frequency distribution of decimal logarithmic values (18,262 daily values) of pressure head at root depth was used to explore the impact of climate change on grassland root water stress. The results were fitted by using a nonparametric kernel probability smoothing distribution. The root water stress analysis was performed within a probabilistic framework, where the probabilities associated with no-stress, stress, and wilting point were considered (Table 4). The analysis indicates that the historic conditions resulted in equivalent percentages of root water stress and no stress conditions (about 49% and 47%, respectively) as well as low probabilities (4.3%) to reach wilting points. While the root water stress probabilities under the mean and wet climate projections were equivalent to the historic condition, the wilting point probabilities under the projections nearly double from 4.3% to 8.5% (Table 4). In contrast, the dry scenario analysis indicates that while the stress probability remains equivalent, the no stress probability is considerably reduced from 48.8% to 35.1%. Of particular interest is also a significant increase in the probability of wilting point (Figure 8). The wilting point probability under the dry climate scenario increased about four times compared to the historic condition. This analysis indicates that the conditions that would be induced by the dry scenario will increase the risk for grass to die and increase the chance of desertification of this region substantially. The root water stress analysis was performed within a probabilistic framework, where the probabilities associated with no-stress, stress, and wilting point were considered (Table 4). The analysis indicates that the historic conditions resulted in equivalent percentages of root water stress and no stress conditions (about 49% and 47%, respectively) as well as low probabilities (4.3%) to reach wilting points. While the root water stress probabilities under the mean and wet climate projections were equivalent to the historic condition, the wilting point probabilities under the projections nearly double from 4.3% to 8.5% (Table 4). In contrast, the dry scenario analysis indicates that while the stress probability remains equivalent, the no stress probability is considerably reduced from 48.8% to 35.1%. Of particular interest is also a significant increase in the probability of wilting point (Figure 8).
The wilting point probability under the dry climate scenario increased about four times compared to the historic condition. This analysis indicates that the conditions that would be induced by the dry scenario will increase the risk for grass to die and increase the chance of desertification of this region substantially.  Figure 7. Monthly averages of precipitation (blue line), actual evapotranspiration (green), and potential groundwater recharge (GR; black line) under (a) historical records, (b) dry, (c) mean and (d) wet climate projections. P and ETa are referred to in the primary y-axis while GR is referred to in the secondary y-axis.
The root water stress analysis was performed within a probabilistic framework, where the probabilities associated with no-stress, stress, and wilting point were considered (Table 4). The analysis indicates that the historic conditions resulted in equivalent percentages of root water stress and no stress conditions (about 49% and 47%, respectively) as well as low probabilities (4.3%) to reach wilting points. While the root water stress probabilities under the mean and wet climate projections were equivalent to the historic condition, the wilting point probabilities under the projections nearly double from 4.3% to 8.5% (Table 4). In contrast, the dry scenario analysis indicates that while the stress probability remains equivalent, the no stress probability is considerably reduced from 48.8% to 35.1%. Of particular interest is also a significant increase in the probability of wilting point (Figure 8). The wilting point probability under the dry climate scenario increased about four times compared to the historic condition. This analysis indicates that the conditions that would be induced by the dry scenario will increase the risk for grass to die and increase the chance of desertification of this region substantially.

Conclusions
The study examined the relationship between historic (1950-2000) and projected future (2000-2100) climate forcings and the GR rates in a site-specific study in the Nebraska National Forest of the Nebraska Sand Hills region. Toward this objective, we considered 65 projections of GR that included four greenhouse gas emissions scenarios. These climate inputs were used in the local-scale numerical model of soil-atmosphere-vegetation dynamics based on HYDRUS-1D. The results have been averaged over decadal, annual, and sub-annual time scales within a probabilistic framework accounting for uncertainty in climate projections. Dry, mean, and wet GR recharge scenarios were defined from the statistical parameters of cumulative GR from daily model runs. The historic climate in the NSH is classified as semi-arid with annual-average potential GR rate of 9.02 cm yr −1 . Projections of decadal precipitation appear to remain quasi-steady with negligible (1%) increase while increasing forecasts of ET 0 range from 8.6% to 31.4%. These main factors are shown to cause a decrease in decadal GR rates ranging from 21% to 39% in extreme cases. It is also worth noting that the Hargreaves approach to estimate ET 0 can inject uncertainty in recharge estimation. In this analysis, the results indicate that GR will likely be overestimated if ET 0 is estimated using the Hargreaves equation. In contrast, GR is likely to be even lower if ET 0 is estimated by using the Penman-Monteith equation.
The results highlight the need to consider adaptation and mitigation measures to respond to projected decreases in GR, which can exacerbate groundwater depletion in the High Plains Aquifer. The outcomes also underscore the need for exploring measures to combat the projected potential increases in wilting probability (~8.3% to 18.5% compared to 4.3% historical baseline), to protect the native grassland ecosystem, and to prevent desertification processes.

Funding:
The authors acknowledge the financial support provided by the NSF IGERT program (DGE-0903469), the Nebraska Geological Society's Yatkola-Edwards scholarship and the American Association of Petroleum Geologists.

Conflicts of Interest:
The authors declare no conflict of interest.