Evaluating the Risks of Groundwater Extraction in an Agricultural Landscape under Di ﬀ erent Climate Projections

: Groundwater resources worldwide are being depleted at alarming rates since 1960 to support agriculture, industry, and domestic water demand. Water harvesting and the implementation of reduced application or more e ﬃ cient irrigation technologies were identiﬁed as two of the most e ﬃ cient practices to mitigate the declining patterns on groundwater resources. However, prior to implementing these practices, understanding how groundwater interacts with surface water and responds to natural and anthropogenic stressors is crucial. Integrated modeling tools that are able to exchange ﬂuxes in both domains are needed to assess how conservation practices will a ﬀ ect our water resources under di ﬀ erent projected climate and land use scenarios. This study aimed to evaluate the most likely impacts of current land management practices under the most severe projections of future climate and quantify the potential mitigation e ﬀ ects of three conservation scenarios on the water resources of the Fort Cobb Reservoir Experimental Watershed (FCREW) in western Oklahoma. The semicoupled SWAT-MODFLOW (SWATmf) model was used to simulate the hydrologic responses of the FCREW to a 50% reduction in the irrigation depths and the transition of 50% and 75% of croplands to rangelands under 32 distinct climate projections. Results showed that future climate can drive a reduction in the streamﬂow ( − 18%) and an increase to the depth of the water table (99%–120%) in the western part of the FCREW by the end of the century. The Fort Cobb Reservoir was expected to reduce its release after the mid-2060s to maintain its current target level. All the scenarios, aimed at decreasing groundwater extractions or implementing conservation measures, signaled a full recovery response in the groundwater levels 7–10 years after the year the conservation practices were implemented. The 50% reduction in the irrigation depths was found to elicit faster hydrologic systemic responses than the two that implemented conservation measures, which contravene tradition and would imply cessation of agricultural activities. This study can enable stakeholders to formulate timely adaptation and mitigating strategies to adopt to land use changes.


Introduction
Water stress, the ratio of groundwater withdrawal relative to its recharge rate over a given aquifer, is a potential indicator of regional water insecurity [1]. One of the main drivers of water stress is groundwater abstraction for irrigation that can lead to aquifer depletion [2]. In fact, groundwater resources worldwide are being depleted at alarming rates because of the overwhelming groundwater exploration in the second half of the 20th century to support agriculture, industry, and domestic water demand [3]. Hydrologic studies have shown that this phenomenon has resulted in increased global groundwater depletion of approximately 25% since 1960, which corresponds to almost 39% of the global annual groundwater extractions in 2000 [4]. As a result, excessive groundwater depletion and water stress have been reported in North and South Africa, Central U.S., Northeast China, Australia, and South and Central Asia, among others. Consequently, groundwater depletion is expected to continue as agricultural activities intensify to enhance food production. This condition is further exacerbated by expected increase in water demand due to population growth and climate variability and change. In order to mitigate this issue, more sustainable management practices that can alleviate the declining conditions of groundwater resources are being implemented. The most efficient practices have been water harvesting, aiming to increase groundwater recharge while decrease groundwater extractions by implementing changes in land use and/or transition to more efficient irrigation technologies and crops [5].
Prior to implementing these practices, understanding how groundwater responds to natural and anthropogenic stressors is crucial. Specifically, understanding the interactions between surface water and groundwater and how these interactions are affected by changes in land use and climate are essential in assessing the viability of these practices. For example, changing the land use from row crops to pasture is expected to reduce irrigation usage and, hence, the withdrawal from the groundwater, leading to a recovery in groundwater levels. Climate variability and change, on the other hand, can increase the pressure on groundwater resources because of the reduction in the reliability of surface water supplies driven by increasing temperatures and projected prolonged droughts [6]. However, the impacts of increased evaporation associated with projected increased temperatures can be mitigated by storing surface water excess in groundwater banks [7]. Therefore, an efficient reallocation of water will require a conjunctive management of groundwater and surface water resources because of their vital connection [8].
As the analysis of groundwater-surface water interactions involves a larger spatial domain, the influence of climatic and anthropogenic stressors becomes more relevant. In the catchment and subcatchment scale, for example, these interactions are characterized by the physiographic and climatic configuration of many diverse landscapes [9], making their monitoring and assessment laborious and very costly. For this reason, the lack of data has been the major barrier that impedes the transition from the groundwater development (i.e., allocation of new explorations) to the implementation of a sustainable management of the current allocated groundwater resources [5]. Counteracting these facts, the diversification of monitoring systems at experimental fields, the use of upscaling techniques from small-scale studies, and the acquisition of data by remote sensing have allowed the development of various hydrologic modeling tools. Such models have led to a better understanding of the water circulation in the biosphere at the surface and subsurface domains and provide an alternative to the lack of field data at a large scale. Although this may be true, most of the specialized groundwater models are typically not connected to the surface; thus, the impacts of surface alteration are not properly propagated to the groundwater domain. Similarly, hydrologic models that are capable of representing changes brought about by agricultural activities do not have a realistic representation of the groundwater domain. Computational efforts have been made towards coupling two or more models that are able to partially exchange fluxes in both domains [10][11][12][13][14][15][16][17]. These integrated tools have been tested to be reliable models but have not widely been used to project possible future conditions of water resources under integrated climate and land use scenarios.
To propose comprehensive management practices for the sustainable use of our water resources, there is a critical need to understand the dynamics between groundwater and surface water [12]. The main objective of this study was to assess the changes in the groundwater under future climate and land management schemes. Specifically, this study was geared to (1) develop a hydrologic distributed Water 2020, 12, 400 3 of 20 modeling framework to simulate groundwater-surface water interactions at the watershed scale, (2) determine the effects of future climate scenarios on fluctuation and probability of occurrence of streamflow, groundwater level, and reservoir storage, and (3) quantify the potential effects of land management scenarios on the changes in the temporal patterns of these resources. The semicoupled Soil and Water Assessment Tool-Modular Finite-difference Flow (SWATmf; [15]) model was used to simulate the hydrologic responses of the Fort Cobb Reservoir Experimental Watershed (FCREW) in central western Oklahoma to different management schemes intended to mitigate surface groundwater depletion. Groundwater extraction in the FCREW is mainly used for irrigation of row crops, but the projected climate changes can further aggravate the decreasing groundwater levels. To determine the feasibility of groundwater level recovery for FCREW, we developed tools and techniques to evaluate the projected mid-and long-term states of the water resources in the watershed. Furthermore, this approach can be used by stakeholders in other areas to assess the impacts of proposed land management practices on surface-groundwater interactions under projected climates.

SWATmf Implementation
The changes in water resources under future climate and land use scenarios in the FCREW ( Figure 1) were assessed using the hydrologic modeling framework SWATmf [15]. This framework coupled the models SWAT rev. 664 [18] and MODFLOW-NWT v. 1.1.3 [19] allowing exchange of fluxes between the surface and the subsurface model domains.
Water 2020, 12, x FOR PEER REVIEW 3 of 21 modeling framework to simulate groundwater-surface water interactions at the watershed scale, (2) determine the effects of future climate scenarios on fluctuation and probability of occurrence of streamflow, groundwater level, and reservoir storage, and (3) quantify the potential effects of land management scenarios on the changes in the temporal patterns of these resources. The semicoupled Soil and Water Assessment Tool-Modular Finite-difference Flow (SWATmf; [15]) model was used to simulate the hydrologic responses of the Fort Cobb Reservoir Experimental Watershed (FCREW) in central western Oklahoma to different management schemes intended to mitigate surface groundwater depletion. Groundwater extraction in the FCREW is mainly used for irrigation of row crops, but the projected climate changes can further aggravate the decreasing groundwater levels. To determine the feasibility of groundwater level recovery for FCREW, we developed tools and techniques to evaluate the projected mid-and long-term states of the water resources in the watershed. Furthermore, this approach can be used by stakeholders in other areas to assess the impacts of proposed land management practices on surface-groundwater interactions under projected climates.

SWATmf Implementation
The changes in water resources under future climate and land use scenarios in the FCREW ( Figure 1) were assessed using the hydrologic modeling framework SWATmf [15]. This framework coupled the models SWAT rev. 664 [18] and MODFLOW-NWT v. 1.1.3 [19] allowing exchange of fluxes between the surface and the subsurface model domains.   The FCREW watershed was used as the study area based on its predominantly agricultural landscape, dependency on groundwater resources to sustain productivity, its reservoir playing a key role for eco-hydrological regulation, and the availability of a robust monitoring network and long-term data ( Figure 1) [20]. In accordance with the 2011 National Land Cover Database [21], cultivated crops and grassland/herbaceous comprise approximately 60% and 29% of the FCREW area, respectively. The watershed's main water source for irrigation is groundwater extracted from the Rush Springs Aquifer (RSA), on which the FCREW rests. The Fort Cobb Reservoir (FCR) was constructed in 1959 as a multipurpose structure, serving not only as source of public water supply for Anadarko and Chickaska, Oklahoma, but also for flood control, recreation, and to sustain aquatic biodiversity [20,22]. However, in 2012 the Fort Cobb Master Conservancy District was fully allocated, which means it would not be able to meet future additional demands [23]. Moreover, the 2012 Oklahoma Comprehensive Water Plan [24] projected a significant groundwater depletion by 2060 in the RSA. Additionally, in some areas of the nearby High Plains Aquifer (HPA), mainly devoted for irrigation, depletion in the groundwater levels, and a reduction in the saturated thickness of more than 30 m and 50%, respectively, has been observed [3]. Despite these declines, McGuire et al. [25] and McGuire [26] have reported that the introduction of agricultural conservation programs together with the improvement of irrigation efficiency during the last decades have produced a recovery of the water levels in portions of the HPA.
SWATmf is a unified modeling framework in which SWAT and MODFLOW representing the surface and groundwater domains, respectively, are independently discretized and parametrized. A detailed description of these two models' conceptualization, evaluation, and limitations can be found in Arnold et al. [27], Gassman et al. [28], and Harbaugh [19]. In the case of SWAT, the model is typically set up using ArcSWAT, which first discretizes the watershed into sub-basins using elevation-based flow direction and accumulation maps. Then, unique combinations of geospatial data of land use, soil type, and slope are generated to describe the spatial heterogeneity of the study area. These combinations, represented as irregular fractions of sub-basin area, are commonly known as Hydrologic Response Units (HRUs). On the other hand, MODFLOW requires the groundwater domain to be discretized into stratified structured grids in accordance with the properties of the geological formation that describe the three-dimensional groundwater flow equation. These include the vertical and horizontal hydraulic conductivities, specific storage, and specific yield as well as some specific boundary and initial conditions that constrain the system. For the surface model, SWAT, the watershed was subdivided into 12 sub-basins and then into 1000 HRUs. The RSA aquifer (groundwater domain), simulated by MODFLOW, was discretized into two layers [15] using 280 m cells. In the MODFLOW model, the following processes were represented: (1) water fluxes between the reservoir and the RSA using the reservoir package (RES) including the daily reservoir water elevation and the bathymetry, sediment thickness, and hydraulic conductivity of the reservoir bed; (2) well water extractions using the well package (WEL). Note that because of the ungauged characteristics of the well extractions in this watershed, the withdrawn water volume was estimated using the SWAT auto-irrigation model and scaled to a well influence area based on central pivot coverage. This was implemented to avoid the possible inclusion of HRUs with the same properties but were outside the well influence area; (3) water fluxes within cells using the Newton solver and layer properties (e.g., vertical and horizontal hydraulic conductivities, specific yield, etc.) specified in the Upstream Weighing Package (UPW); and (4) recharge fluxes on the top of the domain using the recharge package (RCH) based on the SWAT percolation depths to the deep aquifer at the HRU level. Thus, the modeling framework used in this research, SWATmf, obeyed a semicoupled scheme in which the deep percolation fluxes and irrigation volumes computed in SWAT were used to define the recharge and well extraction rates in MODFLOW. It is important to note that feedbacks from MODFLOW to SWAT were not considered.
Geologic data reported by Becker and Runkle [29], Penderson [30], and Guzman et al. [15] were initially used to describe the physical properties of the groundwater domain (i.e., hydraulic conductivities, specific yield, and specific storage) required to setup the MODFLOW UPW package.
The soil properties, river conductance, and reservoir-bed hydraulic conductivity were estimated based on the Soil Survey Geographic Database (SSURGO; [31]). Likewise, a total of 655 extraction wells primarily used for irrigation were identified within FCREW from the OWRB Groundwater wells database [32] (Figure 1).

Reservoir Model
Considering the important ecosystem services (e.g., water supply, fish and wildlife propagation, recreation) that the FCR supports, this study aimed to assess the combined effects of the changes in climate and irrigation demands on the reservoir storage. Some routines were modified in SWATmf to better represent the reservoir storage and release under future scenarios based on historic records of pool elevation, current release regulations, and projected withdrawals for drinking water supply systems in the surroundings of the study area. The reservoir water balance was determined using a controlled release scheme with a target pool elevation as follows: where S k+1 and S k are the reservoir storage volume at time steps k+1 and k, F k in is the reservoir inflow from streams at time step k, O k v is the reservoir inflow from the overland flow at time step k, P k is the precipitation volume falling over the reservoir surface area at time step k, E k v is the evaporated volume from the reservoir surface area at time step k, W k is the water withdrawals at time step k, S k e is the seepage at time step k, F out_max is the maximum allowed outflow, F k out is the outflow at time step k, and a and b are calibration parameters that incorporate abstractions of other kind of losses. Reservoir releases through the spillway were allowed when the pool elevation at the end of each simulation day (h k ) was greater than the defined target pool elevation but less than the maximum allowed discharge (h ta and F out_max ; 409 m above mean sea level (mamsl) and 1300 cfs by the United States Army Corps of Engineers; [33]). In this case, the released volume corresponds to the difference between the storage volume at the end of the simulation day (S k ) and the storage volume at the target pool elevation (S ta ). Otherwise, the released volume equals the maximum allowed discharge.
The elevation-area-storage curves for the FCR were constructed from the 1993 sediment survey [34]. Similarly, records and projections of population growth and future water demand of Anadarko and Chickasha were obtained for the period 2000-2060 [35] and extended until 2100 using nonlinear regression models. These were employed to compute the withdrawals from the reservoir. The aforementioned elevation-area curves combined with the precipitation records and the SWAT-generated potential evapotranspiration depths were used to determine the precipitation and evaporation volumes at each simulation day. Additionally, the contribution from streams and overland flow was assumed to be the sum of discharges coming from the Cobb, Lake, and Willow Creeks and the surface runoff produced at the sub-basins located at the FCR vicinity, respectively.

SWATmf Calibration
Calibration of the SWATmf model was conducted using a two-step process [36]: (1) calibration of the surface domain using daily streamflow data at three gauging stations and pool elevation in the Fort Cobb reservoir, and (2) calibration of the groundwater domain using three continuous piezometric heads on a daily basis. In both steps, metrics as well as the three-dimensional groundwater level representation were used to assess the infiltration fluxes at the boundary condition between the two domains. SWAT was first calibrated until acceptable model metrics were achieved. The combination of parameters that resulted in the highest metrics (e.g., Nash-Sutcliffe Efficiency coefficient, NSE; [37]) Water 2020, 12, 400 6 of 20 and lowest errors (e.g., percent bias, PBIAS; [38]) were used to set up SWATmf. Next, the parameters describing the groundwater flow in MODFLOW were adjusted to obtain the best possible agreement between the observed and simulated groundwater levels at the given monitoring points. Also, the spatial patterns of groundwater heads were verified to ensure that constraints, such as groundwater levels not reaching the ground surface, were simulated properly. If unrealistic conditions were detected, SWAT was recalibrated by locally adjusting the groundwater parameters (e.g., RCHRG_DP, GWQMIN) in specific HRUs. The performance of SWATmf was then reassessed in a new iteration until the model performance was validated in both domains.
The performance of the SWATmf model for the FCREW was evaluated using a ten-year period (2008 to 2017) of daily streamflow and groundwater level from the monitoring stations. Three years of this period were used to set up the initial conditions of the model. The performance of SWATmf in representing the response of FCREW was assessed using NSE and PBIAS. The observed streamflow at the Lake, Cobb, and Willow Creek gauges and groundwater levels at the Eakly, Core2, and Alfalfa monitoring wells were used to evaluate model performance. The simulated reservoir volume storage was also calibrated against daily time series provided by the USACE Tulsa District [39] and applying the same metrics.

Climate Projections
The impacts of future climate on water resources of the FCREW were evaluated using the CMIP5 dataset downscaled using the Localized Constructed Analogs (LOCA) technique [40,41]. These datasets resulted from a multimodel ensemble consisting of 32 Integrated Atmosphere-Ocean Global Climate Model (AOGCM) projections for the most aggressive representative concentration pathway (i.e., RCP8.5). The RCP8.5 depicts a scenario with the highest greenhouse gas emissions based on assumptions about a high global population growth in absence of mitigation policies [42]. The dataset covers the period 2006-2099, and it is available for the contiguous U.S. territory and portions of southern Canada and northern Mexico [41]. This dataset was selected because it aggregates daily projections of precipitation and minimum/maximum temperature at a resolution of 1/16 • latitude-longitude (i.e., approx. 7 × 7 km at the equator) for a significant ensemble of models. Therefore, it is expected to be sufficient to generate climate-driven scenarios at the watershed scale for hydrologic models as SWATmf that run on a daily temporal scale.
In accordance with the climate ensemble data for the FCREW domain, the projected trend of total annual precipitation was nearly flat, describing a small decrease of just 2.1 mm per decade from 2018 until the end of the century with values ranging between 640 and 798 mm ( Figure 2). On the contrary, the mean annual temperature presented an upward linear trend at a rate of 0.62 • C per decade for the projection period. This represents a rise of approximately 6.3 • C in the mean annual temperature at the end of the century in comparison with the historic mean (15.9 • C; Figure 2). Previous studies have reported similar projected trends at the watershed and state scales [43,44].

Integrated Climate and Land Use Scenarios.
Complex human-environmental interactions define the functioning of many ecosystems, where certain anthropogenic factors can result in long-term impacts that degrade the environment [45]. Ecologists and hydrologists have attempted to evaluate the effects of changes in the land management practices on water availability and circulation and ecosystem biodiversity [46]. The development and evaluation of future land use change scenarios is vital not only to understand their impacts but also to provide a strong scientific-based framework for mitigation measures and policy making.
Three land use scenarios were considered in this study ( Figure 3). The first scenario was aimed at decreasing groundwater extractions by increasing the water use efficiency of irrigation, for instance, the transition from sprinkler to drip irrigation or the introduction of deficit irrigation techniques. Scenarios 2 and 3 were aimed at implementing conservation measures intended to lessen water demand, erosion, and sediment yield, among others. For this study, conservation was focused on reducing irrigation application by changing two different percentages (50% and 75%) of the agricultural crops to pasture.

Integrated Climate and Land Use Scenarios.
Complex human-environmental interactions define the functioning of many ecosystems, where certain anthropogenic factors can result in long-term impacts that degrade the environment [45]. Ecologists and hydrologists have attempted to evaluate the effects of changes in the land management practices on water availability and circulation and ecosystem biodiversity [46]. The development and evaluation of future land use change scenarios is vital not only to understand their impacts but also to provide a strong scientific-based framework for mitigation measures and policy making.
Three land use scenarios were considered in this study ( Figure 3). The first scenario was aimed at decreasing groundwater extractions by increasing the water use efficiency of irrigation, for instance, the transition from sprinkler to drip irrigation or the introduction of deficit irrigation techniques. Scenarios 2 and 3 were aimed at implementing conservation measures intended to lessen water demand, erosion, and sediment yield, among others. For this study, conservation was focused on reducing irrigation application by changing two different percentages (50% and 75%) of the agricultural crops to pasture.
The three land management scenarios were combined with the 32 climate projections to form a set of 96 scenarios. All the scenarios were applied in 2043 considering a grace period of 25 years from 2018, the year in which the climate projections began. This grace period accounted for transition protocols and activities (e.g., feasibility studies, government regulations) that may be required before implementing changes in the land use. The land management scenarios were assumed to elicit positive impacts on the groundwater resources in terms of depletion and recovery of the groundwater levels. Scenario 1 considered a 50% decrease in the irrigation depth (i.e., from 25 to 12.5 mm) of all HRUs defined as row crops (i.e., AGRR), while Scenarios 2 and 3 exemplified the transition of 50% and 75% of the agricultural lands area (irrigated) into rangelands (i.e., RNGE; no irrigation), respectively. The three land management scenarios were combined with the 32 climate projections to form a set of 96 scenarios. All the scenarios were applied in 2043 considering a grace period of 25 years from 2018, the year in which the climate projections began. This grace period accounted for transition protocols and activities (e.g., feasibility studies, government regulations) that may be required before implementing changes in the land use. The land management scenarios were assumed to elicit positive impacts on the groundwater resources in terms of depletion and recovery of the groundwater levels. Scenario 1 considered a 50% decrease in the irrigation depth (i.e., from 25 to 12.5 mm) of all HRUs defined as row crops (i.e., AGRR), while Scenarios 2 and 3 exemplified the transition of 50% and 75% of the agricultural lands area (irrigated) into rangelands (i.e., RNGE; no irrigation), respectively.

Interpretation of Results
The 32 climate projections were assumed to be equally likely, thus providing a representative sample of possible future patterns of climate from global circulation models. For the streamflow, seasonal cumulative density functions were constructed for three gauging locations to evaluate the relative changes in the seasonal means across scenarios over the two 25-year periods (2043-2067 and 2068-2092) that followed the transition to the new land management practices. For the groundwater levels and reservoir storage, the relative frequency distribution was computed from daily variables in the 32 climate projections for every month from 2018 to 2097 ( Figure 4; steps 1-2). This means that for every month there were at least 960 values (sample size; 30 days × 32 projections) which were then distributed into 11 class intervals to construct the histogram (Figure 4

Interpretation of Results
The 32 climate projections were assumed to be equally likely, thus providing a representative sample of possible future patterns of climate from global circulation models. For the streamflow, seasonal cumulative density functions were constructed for three gauging locations to evaluate the relative changes in the seasonal means across scenarios over the two 25-  Since the historic baseline model was validated in a relatively short period, the analyses of results were performed using the future climate projections. The differences between scenarios were determined relative to the baseline simulated under future climate projection. To simplify the reporting, the likelihood of occurrence was further classified into three main groups: less likely (0.0 ≤ P ≤ 0.1), likely (0.1 < P ≤ 0.2), and most likely (P > 0.2). In the probability maps, the most probable Since the historic baseline model was validated in a relatively short period, the analyses of results were performed using the future climate projections. The differences between scenarios were determined relative to the baseline simulated under future climate projection. To simplify the reporting, the likelihood of occurrence was further classified into three main groups: less likely (0.0 ≤ P ≤ 0.1), likely (0.1 < P ≤ 0.2), and most likely (P > 0.2). In the probability maps, the most probable pathway (highest probability in a month) can be traced for the entire simulation period. This represents the most likely pattern that the variable may take as a result of changes in climate and land management. Statistical tests were also conducted to determine if there were significant differences between the mean of the baseline and that of different scenarios as the first filter to quantify the impacts of land management practices on water resources of the FCREW.

Model Calibration
Overall, the performance of SWATmf in simulating the stream flow and the groundwater levels showed acceptable model performance at all the monitoring stations (Figures 1 and 5). For stream flow, the Lake Creek and Willow Creek presented the same goodness-of-fit (i.e., NSE = 0.64) at daily time scale, while monthly streamflow patterns were better described at the Willow Creek with values of NSE greater than 0.92. On the other hand, the NSE ratings for the Cobb Creek showed daily streamflow NSE = 0.54 and monthly NSE = 0.68. The PBIAS on daily time step was less than 5% at the three monitoring gauges. This indicated an acceptable averaged representation of the water balance by the model during the calibration period (Figure 5a). However, note that the model showed a systematic deficit in the Willow Creek streamflow estimations at the beginning of the simulation and a systematic surplus at the end of the calibration period. A reasonable explanation to this model behavior was the misrepresentation of baseflow patterns that led to an underestimation of peak flows greater than 1.7 m 3 /s.
For the groundwater domain, daily NSE values were equal to 0.50, 0.63, and 0.69 for the Eakly, Core2, and Alfalfa monitoring wells, respectively, while the PBIAS ranged between −0.1% and 0.2%. Overall, the model was able to simulate the general trends in groundwater levels especially the recovery in 2015 after the shift in the climatic regime from dry to wet period. However, short-term fluctuations (interannual) were not properly represented by the model. This was explained as the smoothing effect of recharge patterns occurring across the watershed that resulted from the spatial discretization of clustered, spatially disconnected HRUs derived from ArcSWAT [36]. HRUs are created based on unique combinations of soil, land use, and topography that can result in spatially disconnected areas clustered to a single SWAT process representation. Consequently, in a cluster of HRUs one or more of these unique combinations may represent contiguous or discontiguous locations within the same sub-basin with dissimilar hydrological responses treated as equal by the SWAT model.
The calibration of the reservoir storage showed a good agreement between the measured and simulated volume (NSE = 0.84; PBIAS = −0.27%) except for the underestimation of a peak in May 2014. This disagreement between the observed and simulated reservoir volume was most likely the result of a local precipitation event that went undetected by the Micronet or Mesonet network as a result of its coarse resolution. For the groundwater domain, daily NSE values were equal to 0.50, 0.63, and 0.69 for the Eakly, Core2, and Alfalfa monitoring wells, respectively, while the PBIAS ranged between −0.1% and 0.2%. Overall, the model was able to simulate the general trends in groundwater levels especially the recovery in 2015 after the shift in the climatic regime from dry to wet period. However, short-term fluctuations (interannual) were not properly represented by the model. This was explained as the smoothing effect of recharge patterns occurring across the watershed that resulted from the spatial discretization of clustered, spatially disconnected HRUs derived from ArcSWAT [36]. HRUs are created based on unique combinations of soil, land use, and topography that can result in spatially disconnected areas clustered to a single SWAT process representation. Consequently, in a cluster of HRUs one or more of these unique combinations may represent contiguous or discontiguous locations within the same sub-basin with dissimilar hydrological responses treated as equal by the SWAT model.
The calibration of the reservoir storage showed a good agreement between the measured and simulated volume (NSE = 0.84; PBIAS = −0.27%) except for the underestimation of a peak in May 2014. This disagreement between the observed and simulated reservoir volume was most likely the result

Evaluation of Scenarios
Historic records of streamflow and groundwater levels followed the bimodal regime of the precipitation with the highest values observed in May-June and September-October. For these periods, the mean monthly streamflow ranged from 0 to 3.0, 0 to 1.1, and 0 to 0.9 m 3 s −1 at the Cobb, Lake, and Willow Creek monitoring gauges, respectively. In 2015, the FCREW experienced a significant shift in its climate regime moving from a dry period lasting approximately six years (2009-2014) to a wet period that started with an extreme rainfall event at the second trimester of 2015. Reservoir levels and groundwater elevations were significantly affected during the dry period, presenting a drop of nearly 29% in the reservoir conservation pool and 14%, 33%, and 41% in the water table depth at the Eakly, Core2, and Alfalfa monitoring wells, respectively. However, a recovery of 42%-70% in the groundwater levels was observed in the following three years after the shift (2015-2017).
The three land use scenarios were expected to result in the recovery of groundwater levels in the FCREW as a result of the reductions in well extractions. These scenarios were also expected to reduce the contribution of irrigation to runoff, potentially generating a decrease in the streamflow and reservoir storage. The most significant changes in the streamflow means are shown in Tables 1 and 2 and Figure 6, and the most likely events for reservoir storage and groundwater levels are shown in Figures 7 and 8. For the latter, the darker colored areas correspond to the most likely events to occur in the next decades. The dotted lines depict the primary or clearly defined continuous pathways (likely to most likely events) and the secondary or likely to less likely alternative pathways.

Baseline
The statistical tests performed on the means of seasonal streamflow under current land use and varying climate (baseline; Sc0) showed different hydrologic responses in the eastern, central, and western parts of the FCREW. The analysis of mean annual streamflow under projected climate at Willow Creek (eastern part), near the reservoir (Figure 1), indicated an upward trend with a maximum increase of 44% to 52% from the 2018-2042 mean in the following two 25-year periods, respectively ( Table 1). These changes in the mean were observed across all seasons but were on average seven percent greater in winter. In Cobb Creek (western part), streamflow did not change significantly for the 2043-2067 period but were expected in the winter and summer of the 2067-2092 period. Lake Creek manifested the least change compared to the other two sub-watersheds (Table 1). Overall, if the land management remained the same, climate will negatively impact the streamflow in the eastern part of the FCREW, where the surface domain of the hydrologic system seemed to be more sensitive to the increasing temperature trends predominant at the last three decades of the century (+0.75 • C decade -1 ). Despite the fact that there were significant changes in the mean streamflow at the Willow Creek (eastern part), its area just represents 9% of the FCREW, while the Cobb Creek in the west corresponds to almost 43%. It is also important to note that the key component in this observed increment in the mean streamflow in Willow Creek was the irrigation, which was triggered more often to overcome the increasing evapotranspiration demand. Projected climate exacerbated the existent dry conditions of the FCREW in the west in response to long-term warming trends rather than decreasing precipitation patterns, which according to the CMIP5 RCP8.5 projections were not expected. Similar findings regarding projected climate have been discussed in other studies at the local and state scales [43,44].

Interscenario Analysis
The statistical test between the baseline and the three land management scenarios after their implementation in 2043 suggested that there were significant differences in the seasonal means of the baseline and those of the different scenarios. The differences were further explored by comparing the changes in the cumulative density function between the baseline and the scenarios focusing in the low, middle, and high streamflow values defined by the 25th and 75th percentiles.
Between the baseline and Scenario 1 (-50% irri. depth), the mean annual streamflow increased up to 18%, 9%, and 23% at the Cobb, Lake, and Willow Creek sub-watersheds, respectively (Table 2). This outcome contradicted the initial hypothesis on the streamflow reduction due to a potential decrease in the contribution of irrigation to runoff. The 50% reduction in the irrigation depth considered in Scenario 1, rather than reducing the amount of extracted water for irrigation in a fixed period, reduced the number of days between applications. This means that the SWAT auto-irrigation schemed applied less water but was triggered more often to sustain plant growth. This led to a major extraction of water from the aquifer system in the long-term and a more frequent rewetting of the soil profile, increasing the contribution to the runoff and baseflow [47].
Seasonal cumulative density functions of streamflow for the 2043-2067 and 2068-2092 periods under the different scenarios are shown in Figure 6. Lake Creek's response was similar to Cobb Creek and was not shown in Figure 6. It can be seen that Scenario 1 affected streamflow values predominantly above the 25th percentile except for the 2068-2092 period at the Willow Creek sub-watershed and winter 2043-2067 at the Cobb Creek, where streamflow distribution was shifted to the right by approximately the same amount. Based on the results, a higher irrigation frequency can increase the probability of occurrence of higher-magnitude streamflow. Scenarios 2 and 3, representing the transition of 50% and 75% of agricultural irrigated lands into non-irrigated rangelands, respectively, generated a negative impact on the streamflow. Both scenarios signaled a reduction in the mean streamflow over the entire hydrologic year. However, this reduction, understood as a relative change from the baseline means, were statistically significant for all the seasonal values −12% to −3% at Willow Creek and only for Scenario 3 in winter and spring at the Cobb Creek (−10% to −8%) ( Table 2). The cumulative density functions ( Figure 6) confirmed that streamflow means did not differ significantly from the baseline in the Lake Creek and that its distribution was not expected to be affected as well. For the Willow Creek, on the other hand, the decreasing patterns on streamflow were slightly more significant for all values above the 25th percentile. A similar behavior can be observed for Cobb Creek over 2068-2092 ( Figure 6).
increase the probability of occurrence of higher-magnitude streamflow. Scenarios 2 and 3, representing the transition of 50% and 75% of agricultural irrigated lands into non-irrigated rangelands, respectively, generated a negative impact on the streamflow. Both scenarios signaled a reduction in the mean streamflow over the entire hydrologic year. However, this reduction, understood as a relative change from the baseline means, were statistically significant for all the seasonal values −12% to −3% at Willow Creek and only for Scenario 3 in winter and spring at the Cobb Creek (−10% to −8%) ( Table 2). The cumulative density functions ( Figure 6) confirmed that streamflow means did not differ significantly from the baseline in the Lake Creek and that its distribution was not expected to be affected as well. For the Willow Creek, on the other hand, the decreasing patterns on streamflow were slightly more significant for all values above the 25th percentile. A similar behavior can be observed for Cobb Creek over 2068-2092 ( Figure 6). Generally, results suggested that streamflow will be affected by changes in the irrigation practices. The 50% reduction in the irrigation depth (Scenario 1) will induce a higher application frequency and, hence, favor the contribution of runoff and baseflow to streams, which signaled positive relative differences from the baseline up to 24% at the Cobb and Willow Creek subwatersheds. This increment was observed to be more significant in streamflow values greater than the 25th percentile of the baseline distribution. Therefore, higher-magnitude events can become more frequent at the FCREW requiring a major reservoir release during the flooding season. The reduction in the irrigated area under Scenarios 2 and 3 indicated a decrease in streamflow of up to 12% across Generally, results suggested that streamflow will be affected by changes in the irrigation practices. The 50% reduction in the irrigation depth (Scenario 1) will induce a higher application frequency and, hence, favor the contribution of runoff and baseflow to streams, which signaled positive relative differences from the baseline up to 24% at the Cobb and Willow Creek sub-watersheds. This increment was observed to be more significant in streamflow values greater than the 25th percentile of the baseline distribution. Therefore, higher-magnitude events can become more frequent at the FCREW requiring a major reservoir release during the flooding season. The reduction in the irrigated area under Scenarios 2 and 3 indicated a decrease in streamflow of up to 12% across the watershed. Although, it may not directly affect agricultural production, as it predominantly relies on the groundwater resources, reduction in streamflow can affect the Fort Cobb reservoir operations. In fact, Garbrecht and Schneider [48] reported that the difference in the reservoir inflow during dry and wet periods in FCREW has been recorded to be of up by 100%.

Baseline
The likelihood for change in the groundwater levels followed similar patterns to those observed for the streamflow under the integrated climate and land use scenarios. Generally, the most likely path of GW levels in Core2, located in the east of FCREW, under all the scenarios evaluated in this study did not show significant differences from the baseline. It described an upward trend right after 2028 with a high probability (P > 0.47), setting the water levels close to saturation at the upper boundary of the GW domain in 2039, where it stayed until the end of the century. The recovery response in Core2 followed the increase in the mean annual streamflow (up to 52%) observed at the Willow Creek sub-watershed, confirming the relevant role of the surface water-groundwater interactions in the budget of the RSA [49,50].
The baseline for the Eakly monitoring well in the west, similar to the Cobb Creek streamflow, also located in the west, described a marked decreasing trend (most likely; 0.2 < P < 0.53) projecting a depletion in the groundwater levels equivalent to an increase in the water table depth of 99%-120% for the 2018-2097 period (Figure 7a; Sc0). This depletion was more pronounced during the first three decades of the projection, with the GW level changing up to -64% in 2048 compared to that in 2018 (-2.07 m decade -1 ), while for the remaining five decades the most likely path showed a depletion of just 21%-34% (-0.68 m decade -1 ). This behavior was expected at this location because of the additional pressure in the aquifer system to cover the higher evapotranspiration demand driven by the projected increasing temperatures (+6.3 • C at the end of the century). Moreover, the western part, where this well is located, was the driest, with some years recording 10%-22% less of precipitation than the central and eastern counterparts. Although the 99%-120% (7-10 m) projected increase in the water table depth may have seemed low compared to the mean aquifer thickness (84 m), it actually represented the combined effect of almost three of the hypothetical sustained ten-year droughts estimated by Ellis [51] for FCREW. Similar depletion patterns have been observed in the High Plains Aquifer and are directly associated with an overtopping of the aquifer storage beyond the rate of available recharge [52].
The most likely patterns of groundwater levels under the baseline in the Alfalfa monitoring well were more stable, fluctuating between -7% and 41% from the 2018 level through 2098 but with an upward trend (+3.8 m decade -1 ) during the 2026-2036 period and a depleting pattern from 2046 to 2065 (-1.7 m decade -1 ; Figure 7b; Sc0). The closeness of the Fort Cobb Reservoir, whose target elevation pool was set at 409 mamsl, to Alfalfa might have influenced the stability of the observed level. Additionally, notice that this monitoring well was located at the central part of the watershed as the Lake Creek gauging station, which did not signal significant changes in its mean under the climate projections.

Interscenario Analysis
For Scenario 1, the 50% reduction in the irrigation depth produced a recovery in the groundwater levels at the Eakly and Alfalfa monitoring wells approximately 5-7 years after the scenarios were applied (2043) and with probabilities ranging from 0.15 to 0.56 (likely to most likely). This recovery represented an increase in the water table of 100% in the following 27 and 6 years, respectively (Figure 7a,b; Sc1). This outcome might appear contradictory because, despite that the irrigation depth was reduced by 50%, more water was extracted in the long-term as the irrigation frequency increased to cover additional evapotranspiration demands driven by the increasing temperatures. However, the reduction in the irrigation depth lessened the extraction rates, putting them in equilibrium with the water releasing capacity of the RSA and hence relieving the excessive pressure over the system observed in the baseline. Moreover, the water movement in the soil profile became more dynamic, favoring the recharge into the aquifer (water recycling) [7,47,53]. The quick response observed at the Alfalfa monitoring well (17.3 vs 5 m decade -1 in Eakly) is explained by the denser extraction network in the area, which in case of a generalized reduction in the irrigation depth translated into a fast pressure release in the aquifer system. The GW levels at Eakly also depicted a secondary path (less likely to likely; 0 < P < 0.18) after 2052, projecting an increasing depletion until the end of the century, which did not significantly differ from the observed in the baseline's most likely path. Note that the ensemble describing the recovery pathway as most likely indicated that this kind of measure paved the way to a possible mitigation of the effects of climatic and anthropogenic stressors on the groundwater resources. The most likely patterns of groundwater levels under the baseline in the Alfalfa monitoring well were more stable, fluctuating between -7% and 41% from the 2018 level through 2098 but with an upward trend (+3.8 m decade -1 ) during the 2026-2036 period and a depleting pattern from 2046 to 2065 (-1.7 m decade -1 ; Figure 7b; Sc0). The closeness of the Fort Cobb Reservoir, whose target elevation pool was set at 409 mamsl, to Alfalfa might have influenced the stability of the observed Figure 7. Groundwater levels (m) at the (a) Eakly and (b) Alfalfa monitoring wells for the baseline (Sc0) and three integrated climate and land use scenarios (Sc1-Sc3). Primary paths (white dotted lines) represent clearly defined continuous pathways that incorporate likely to most likely events, while secondary paths (black dotted lines) describe likely to less likely alternative pathways.
As for Scenarios 2 and 3 (50% and 75% AGRR to RGNE), their most likely paths (0.15 < P < 0.44) showed a recovery in the groundwater levels right after 2058, 15 years after the scenarios where implemented. These results were more significant for Scenario 3 where 75% of the agricultural lands were turned into rangelands, which translated to reduction in irrigated area by 75%. This scenario led groundwater levels close to surface ground level at the Eakly and Alfalfa monitoring wells in the course of 40 and 24 years (3.4 and 2.2 m decade -1 ), respectively (Figure 7a,b; Sc3). Scenario 2, being a reduction in the irrigated area of 50%, showed a recovery in the groundwater levels (depth to water table) of just 22% at Eakly, while at Alfalfa a similar recovery as in Scenario 3 was observed. Notice that in the western part of the watershed, to achieve the same recovery pattern in groundwater levels observed by reducing the 50% of the irrigation depth (Scenario 1), 75% of the cropland needed to be turned into range lands (Scenario 3). The groundwater levels at Alfalfa also described a secondary path as the one observed under the baseline (decreasing trend after 2066) but with a probability of less than 0.15.
Overall, the scenarios considered in this study showed a probable recovery of the depleting groundwater levels in FCREW, especially in the western part. The most likely paths of all the scenarios considered in this study were able to depict an increase in the water table of 100%. However, despite the positive impacts of Scenarios 2 and 3, its implementation may be less effective than Scenario 1, which signaled recovery responses 1-3 decades earlier. Moreover, a land use change of 50%-75% of the watershed area may be less feasible than promoting irrigation application efficiency because of its economic consequences. In fact, the Water for 2060 initiative prioritized education and incentive programs for water-efficient crop irrigation equipment conversion and practices in order to limit the Oklahoma water use by 2060 to that recorded in 2012 [54].

Impacts on Reservoir Storage
Baseline The baseline scenario indicated that the 409 m level defined as the target pool elevation by USACE can be maintained with a high probability (0.26 < P < 0.51) through 2062 (Figure 8; Sc0), with maximum and minimum values for the most probable events within the same range observed for the 2011-2017 period (62.5-105 MM m 3 ). After 2062, the uncertainty was significantly higher, as shown by the widening band, where two possible trajectories (0.1 < P < 0.28) manifested. The first path (more likely) described an abrupt drop of 98% in the pool storage over the course of 15 y, setting it close to the dead pool stage (1.55 MM m 3 ), while the second path (slightly less likely) maintained the level around the target storage (91.1 MM m 3 ) until 2090 but then abruptly dropped towards extremely low levels ( Figure 8; Sc0). This mainly was due to the decrease in the mean annual streamflow at Cobb Creek (main tributary of the FCR) after 2068 as response to the projected increasing temperatures, which also led to higher evaporation rates. The results for the streamflow, reservoir storage, and groundwater levels showed how significantly dry conditions in the eastern part of the FCREW are expected to be exacerbated in the last 3 decades of the projection period.
Interscenario Analysis Scenario 1, driven by changes in the irrigation depth, presented a similar behavior to the baseline (two likely paths after 2062) but favoring the stationarity of the storage within the historic records (62.5-105 MM m 3 ) after 2062 (Figure 8; Sc1). As expected, the statistical tests confirmed that the difference between the means for the most likely paths with respect to the baseline for the 2043-2097 period was significant for Scenario 1, depicting an increase of 57%. The last two scenarios, driven by a reduction in the irrigated area resulting from the conversion of croplands into non-irrigated grasslands, also exhibited a most likely decreasing path in the storage after 2062 (0.1< P < 0.47) (Figure 8; Sc2-3). This reduction in the irrigated area favored a recovery in the groundwater levels but affected the water balance in the surface, lessening the water recycling effect produced by irrigation.
Overall, the FCR was expected to reduce its release after the mid-2060s to compensate for the decline in its inflow coming mainly from the Cobb Creek sub-watershed and for the losses from increasing evaporation and water demand from Anadarko and Chickaska. around the target storage (91.1 MM m 3 ) until 2090 but then abruptly dropped towards extremely low levels (Figure 8; Sc0). This mainly was due to the decrease in the mean annual streamflow at Cobb Creek (main tributary of the FCR) after 2068 as response to the projected increasing temperatures, which also led to higher evaporation rates. The results for the streamflow, reservoir storage, and groundwater levels showed how significantly dry conditions in the eastern part of the FCREW are expected to be exacerbated in the last 3 decades of the projection period. Figure 8. Reservoir storage (m 3 × 10 6 ) for the baseline (Sc0) and three integrated climate and land use scenarios (Sc1-Sc3). Primary paths (white dotted lines) represent clearly defined continuous pathways that incorporate likely to most likely events, while secondary paths (black dotted lines) describe likely to less likely alternative pathways.

Interscenario Analysis
Scenario 1, driven by changes in the irrigation depth, presented a similar behavior to the baseline (two likely paths after 2062) but favoring the stationarity of the storage within the historic records (62.5-105 MM m 3 ) after 2062 (Figures 8; Sc1). As expected, the statistical tests confirmed that the difference between the means for the most likely paths with respect to the baseline for the 2043-2097 period was significant for Scenario 1, depicting an increase of 57%. The last two scenarios, driven by a reduction in the irrigated area resulting from the conversion of croplands into non-irrigated grasslands, also exhibited a most likely decreasing path in the storage after 2062 (0.1< P < 0.47) ( Figure  8; Sc2-3). This reduction in the irrigated area favored a recovery in the groundwater levels but affected the water balance in the surface, lessening the water recycling effect produced by irrigation.
Overall, the FCR was expected to reduce its release after the mid-2060s to compensate for the decline in its inflow coming mainly from the Cobb Creek sub-watershed and for the losses from increasing evaporation and water demand from Anadarko and Chickaska. Figure 8. Reservoir storage (m 3 × 10 6 ) for the baseline (Sc0) and three integrated climate and land use scenarios (Sc1-Sc3). Primary paths (white dotted lines) represent clearly defined continuous pathways that incorporate likely to most likely events, while secondary paths (black dotted lines) describe likely to less likely alternative pathways.

Conclusions
We developed a semicoupled SWAT-MODFLOW model (SWATmf) to assess the future patterns of streamflow, groundwater levels, and reservoir storage from three management practices under the most aggressive climatic scenario (RCP8.5) at the FCREW. Results showed that all the scenarios, aimed at decreasing groundwater extractions or implementing conservation measures, were expected to generate minor impacts on the streamflow. However, the projected climate change drove significantly negative impacts in the streamflow at the western part of the watershed (up to -18% in winter). Thus, the Fort Cobb Reservoir was expected to reduce its release after the mid-2060s to maintain its current target level and compensate for the losses from increasing evaporation and water demand from Anadarko and Chickaska. This could affect all the activities downstream of the dam that rely on the reservoir release.
The groundwater levels were expected to decrease in the western part of FCREW resulting in a probable decrease in the water table by 99%-120% at the end of the century if the current land use were maintained. The mitigating scenarios, determined by changes in the irrigation depth (Scenario 1) and the total irrigated area (Scenarios 1 and 2) signaled a recovery response 7-10 years after changes were applied. The scenarios focusing on reducing groundwater extractions were found to elicit faster hydrologic systemic responds than the two that implemented conservation measures. Furthermore, the introduction of new technologies does not contravene the tradition and current economic activities of the local communities. The reduction of the irrigated area, although efficient, would imply cessation of agricultural activities and conversion of the landscape to a different land use.
The probabilistic approach employed in this study proved to be reliable in identifying future most likely patterns of water resources at the watershed scale from model outcomes due to a diverse set of climate projections. Particularly, this methodology improved the assessment of the environmental risk by providing most likely pathways rather than limiting the analysis to the averaged climate trend as is usually done. It is, however, important to note that, despite the fact the modeling framework used in this study described the spatiotemporal fluctuations of streamflow and groundwater levels well, the model underestimated the baseflow contribution to streamflow. This may be due to the lack of feedback from the MODFLOW domain to the SWAT shallow aquifer. Although the results of this study may be considered for land management planning in the FCREW, the baseline model was validated in a relatively short period (10 years) subjected to the available environmental data.
The framework that we developed in this study facilitated the evaluation of the possible impacts of different land use management under different climate projections and allowed us to estimate the most probable state of our water resources several decades after the implementation of changes. This will enable stakeholders to formulate timely adaptation and mitigating strategies to adopt to natural and anthropogenic changes.