How Do Gaining and Losing Streams React to the Combined Effects of Climate Change and Pumping in the Gharehsoo River Basin, Iran?

Projections of potential impacts of climate change and groundwater abstraction on gaining and losing streams, particularly in ephemeral river basins exhibiting sporadic and intricate flux exchanges, have remained largely unexplored. To fill this gap, we propose a promising modeling scheme based on the new fully integrated hydrological model SWAT‐MODFLOW‐NWT, calibrated and validated for 1978–2012, to quantify the intertwined surface‐groundwater interactions under a conjuncture of three climatic emission scenarios (RCP 2.6, 4.5 and 8.5) and two groundwater pumping variants: “pumping” (extending current groundwater utilization into the future) and “nonpumping” (assuming a complete cease of pumping in the future). By forcing the integrated model with future downscaled climate predictors of CanESM2 under the aforementioned RCPs for three time slices up to year 2100, projections of various water resources components for the Gharehsoo River Basin (GRB), in northwestern Iran were made. Results demonstrate that because of a general decrease of future precipitation, though with ups and downs across the total projection period, most of the surface and ‐subsurface budget quantities and fluxes are substantially affected. In particular, future groundwater discharge (baseflow) to the gaining streams will be more influenced by the “pumping” variant (increasing and decreasing for “nonpumping” and “pumping”, respectively) than the concentrated groundwater recharge from the losing streams (decreasing and increasing for “nonpumping” and “pumping”, respectively). Future water yield and groundwater storage will also diminish and, surprisingly, this cannot be alleviated by future “nonpumping”, indicating the groundwater overutilization is the compelling reason for the future water scarcity in the GRB, rather than climate change alone.


Introduction
According to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change (IPCC), "climate change over the 21st century is projected to reduce renewable surface water and groundwater resources in most dry subtropical regions (robust evidence, high agreement), intensifying competition for water among sectors (limited evidence, medium agreement)" (IPCC, 2013). In fact, water resources have been affected by global change all over the world in a wide verity of unpredictable ways (Green et al., 2011;Mani et al., 2016). Potential impacts of climate change on surface hydrological processes have been discussed in further details, but as emphasized by Holman, 2006, Green et al., 2007, Green et al., 2011, Bovolo et al., 2009, Nkhonjera & Dinka, 2017, Nistor & Mindrescu, 2019, not many studies have drawn attention to impacts posed on the subsurface hydrological processes.
The fluctuations of the flow in gaining and losing streams in a river network connected hydraulically to the aquifers are highly associated with the recharge variabilities (Taie Semiromi & Koch, 2019). In addition, since climate change produces major effects on the recharge rate, providing a reliable estimation of the later necessitates simulation of the atmospheric, surface, and subsurface hydrologic processes using a fully integrated mechanistic hydrological model (P. Goderniaux et al., 2009). Recharge amounts can be estimated using a wide variety of techniques ranging from simple to complex approaches. These are conventionally classified as water budget-based methods (e.g., Tilahun & Merkel, 2009), groundwater modeling techniques (e.g., J. W. Chen et al., 2016), surface-water data-based techniques (e.g., Meyboom, 1961), physical methods: unsaturated zone (e.g., Min et al., 2015), physical methods: saturated zone (e.g., Varni et al., 2013), chemical tracer approaches (e.g., Wu et al., 2016), and heat tracer methods (e.g., Taniguchi & Sharma, 1993). Notwithstanding their differences, none of the aforementioned approaches can address flux and feedback exchanges taking place between the surface and groundwater domains. Due to the fact that the meteoric infiltration (recharge) is reliant upon precipitation and evapotranspiration over the land surface, evapotranspiration within the vadose and saturated zone, and ultimately river-aquifer interactions, this feedback becomes an inherent component of the water cycle. Thus, to accurately quantify these processes, applying a fully integrated mechanistic hydrological model is inevitable. This is owing to the fact that the involved processes, including the flux feedback, can be simulated both simultaneously and separately. Taking into consideration either one part of the hydrological system, namely, surface water or groundwater, will yield unrealistic, inaccurate, and potentially unusable results in the context of climate change impact assessments (P. Goderniaux et al., 2009).
Despite the big advantages of using fully integrated mechanistic models for the assessment of climate change impacts on the subsurface hydrological components, most studies have enjoyed employing simplified schemes thus far for such a purpose. For instance, Z. H. Chen et al. (2002) investigated climate change impacts on a Canadian aquifer using a simple linear transfer function, which connected the groundwater head fluctuations to the variabilities of precipitation and temperature. Obviously, these types of simplifications between recharge and groundwater head introduce considerable uncertainty. To reduce the latter, created as a result of an independently assessment of surface and groundwater domains under a changing climate, fully coupled/integrated physics-based and spatially distributed hydrological models have recently come to the fore (Taviani & Henriksen, 2015). Thus far, a number of fully coupled surface-subsurface hydrological models have been developed such as HydroGeoSphere (Brunner & Simmons, 2012), CATHY (Camporese et al., 2010), ParFlow (Maxwell et al., 2007), InHM (VanderKwaak & Loague, 2001), OGS (Kolditz et al., 2012), and SWAT-MODFLOW-NWT (Bailey et al., 2016). In comparison with other integrated physic-based hydrological models, SWAT-MODFLOW-NWT takes advantages of two models, namely, SWAT and MODFLOW-NWT, where either is the most commonly used model for surface and subsurface 10.1029/2019WR025388

Water Resources Research
TAIE SEMIROMI AND KOCH hydrological modeling, respectively. SWAT-MODFLOW-NWT has been developed in a manner to produce a complementary ability for having the advantages of either of the two models, while overcoming the potential disadvantages when they are used individually (Taie Semiromi & Koch, 2019). Conventionally, surface and subsurface hydrological models have mostly been coupled with the aim of quantifying the interactions between a river network and its adjoining groundwater system (Bejranonda et al., 2007;Kalbus et al., 2006).
Employing physically based hydrological modeling approaches taking into consideration both surface and subsurface processes in order to project the climate change impacts on surface-groundwater interactions has not been well documented as of yet. To that respect, a few studies have investigated the interconnection between a groundwater system and a river network under a changing climate (Goderniaux et al., 2011(Goderniaux et al., , 2015P. Goderniaux et al., 2009;Saha et al., 2017;Scibek et al., 2007). Moreover, a few works have explored the response of groundwater-dependent ecosystems, including wetlands, lakes, and riparian areas to climate change scenarios (Dwire et al., 2017;Havril et al., 2017;Klove et al., 2014;Taviani & Henriksen, 2015;van Engelenburg et al., 2018).
Although in previous studies climate change impacts on one of the compartments of surfacewater or groundwater has been projected, the climate change impacts upon the surface-groundwater interactions, as reflected in gaining (effluent) and losing (influent) streams, groundwater storage, and the ultimate product of the water cycle in a basin, namely, the so-called water yield, have not been assessed till date. In this regard, Scibek et al. (2007) studied groundwater heads and streamflow hydrographs under climate change and groundwater pumping scenarios. However, they did not assess the impacts of climate change and pumping scenarios on the gaining and losing streams, where the flux exchanges between groundwater domain and river network (Alvarez-Villa et al., 2014;Baalousha, 2012). P. Goderniaux et al. (2009), Goderniaux et al. (2011), and Goderniaux et al. (2015 projected climate change impacts on surface and groundwater components using the HydroGeoSphere (HGS) model, a finite-element coupled surface-subsurface flow model ( Brunner & Simmons, 2012). However, the central focus of these studies was to provide groundwater head and streamflow fluctuations under the Special Report on Emission Scenarios (SRES) A2 emission (medium-high) scenario. Thus, a spatiotemporal analysis of the climate change scenario on the surface-groundwater interactions represented by the gaining/losing streams in relation to their adjoining groundwater domain was not conducted. Notwithstanding the previous studies, Saha et al. (2017) explored the groundwater discharge/baseflow generated from gaining reaches under climate change scenarios A2 and B1 of SRES using a physically based and distributed model, that is, Gridded Surface Subsurface Hydrologic Analysis (GSSHA). Although the baseflow, supplied by the gaining streams, was projected for 2020-2040 for different temporal scales, the spatial variability of this water component across the Kiskatinaw River Watershed in British Columbia, Canada, was not evaluated. In addition, the potential impacts of climate change on the concentrated groundwater recharge, produced by the losing streams/reaches, towards the shallow groundwater were not taken into consideration.
Given the mentioned research gaps in the assessment of climate change impacts on surface-groundwater interactions, in the current study, we introduce a new hybrid modeling scheme for a comprehensive assessment of climate change effects on different tangible components of the hydrological cycle, in particular, the surface-groundwater flux exchange. To that end, we will benefit from a recently developed and well-tested fully coupled hydrological model, namely, SWAT-MODFLOW-NWT (Bailey et al., 2016;Taie Semiromi & Koch, 2019), in which the land surface processes, in-stream, vadose zone, and subsurface hydrological processes are efficiently and consistently coupled. In the light of using this coupled model, highly resolved spatiotemporal outputs of baseflow/groundwater discharge resulting from the gaining streams, as the conventional output of most hydrological models, in response to various climate change scenarios, can be projected. Additionally, the same projection can be made for the losing streams where much of the groundwater recharge takes place. In addition to the surface-groundwater exchange processes, this model enables one to assess a wide range of other water resources components of a catchment such as groundwater storage, total streamflow, and water yield in response to different stressors such as climate change effects.
There is an urgent necessity for climate impact and hydrological studies in the Gharehsoo River Basin (GRB), located in northwestern Iran where an intermittent/seasonal drainage network dictates a complex surface-groundwater interaction mechanism. First, inasmuch as the study area is located in a mountainous region, it is expected to be more sensitive to global climate change (Huber et

Water Resources Research
is the case for many mountainous areas of the world (e.g., the Andes and the Himalayas) (Hu et al., 2013). Second, as reported by IPCC (2013), in many mid-latitude and subtropical dry regions such as the GRB, the mean precipitation will most likely decrease in the future. This will have enormous adverse impacts on the water availability in the GRB. For these reasons, projections of spatiotemporal impacts of climate on a wide range of water resources components, particularly those supplied by gaining and losing streams, under different emission scenarios are very important for proposing effective adaptation and mitigation strategies in this region.
Third, the fragile balance between water supply and demand has recently been lost in the GRB, and supposedly, the ongoing climate change impacts will further worsen the situation. To that respect, Ardabil Regional Water Authority (2013) has expressed a need for climate change impacts assessment studies on the various water resources components of this basin. This is especially required in the GRB where the low flows, as a highly transient flow, are dependent on the sporadic and intricate interactions between the shallow groundwater and losing/gaining streams. Indeed, these low flows can supply a major fraction of agricultural, industrial, and domestic consumptions for this basin, especially, during dry periods, that is, summer.
Thus, the objective of this study is to advance our understanding of gaining and losing streams under climate change impacts by employing a fully integrated hydrological model, SWAT-MODFLOW-NWT, as the latter can resolve several simplifications considered in previously mentioned studies. Above all, in this study, a recently developed state-of-the-art hybrid model employing artificial neural networks (ANNs), discrete wavelet transform (DWT), and quantile mapping (QM) is used for the downscaling of daily precipitation (Taie Semiromi & Koch, 2017). Moreover, this study makes use of the Second Generation Canadian Earth System Model (CanESM2) developed by the Canadian Centre for Climate Modeling and Analysis (CCCma). Climate projections are available under three representative concentration pathways (RCPs), namely, RCP 2.6, RCP 4.5, and RCP 8.5. Moreover, to distinguish man-made effects from climate change impacts on surface-groundwater interactions, the projections will be made under two further variants, namely, (1) "nonpumping" variant (undisturbed system: assuming a complete cease of pumping in the future) and (2) "pumping" variant (disturbed system: extending the current groundwater utilization into the future).

Gharehsoo River Basin
The Gharehsoo River Basin (GRB), situated in the northwest of Iran (Figure 1a), covers an area of 4,193 km 2 . It is located between latitudes 37°46′ to 38°36′N and longitudes 47°46′ to 48°42′E (Javan et al., 2015). Recognized as a semiarid region, dry and wet seasons of the GRB span from May to October and November to April, respectively. Considering the historical daily records obtained from synoptic meteorology stations of the Iran Meteorological Organization and climatology stations of the Iran Ministry of Energy in the study area, the hydroclimate variables, including daily minimum and maximum temperature and precipitation of 12 weather stations (Figure 1b), distributed over the GRB for the observed period 1978-2005 were collected. The statistics and geographical characteristics of the selected weather stations are listed in Table 1 and 2. On average, the annual precipitation across the GRB is around 305.90 mm (Table 1); however, their average varies spatially from 267.11 mm for the plain area/aquifer of the basin to 357.36 mm in the mountainous zone (Figure 1b). The precipitation falls mainly in the cold months, that is, November to April. According to the Köppen climate classification, the prevailing climate in this region is "cold and semiarid" (BSk).

Ardabil Aquifer
The Ardabil plain/aquifer (38°30′ to 38°27′N and 47°55′ to 48°20′E), located within the GRB, covers an area of 1,073 km 2 (Figures 1a-1d). The aquifer is encircled by two high mountains, Alborz and Sabalan. On average, the areal annual precipitation over the aquifer is about 270 mm, and the prevailing climate identifies May and August as the wettest and driest months of the year, respectively (Ardabil Regional Water Authority, 2013; Nourani et al., 2015). With an average temperature of 9°C, the Ardabil aquifer is known for being one of the coldest regions in Iran, as characterized by an average of 130 freezing days in a year 10.1029/2019WR025388

Water Resources Research
TAIE SEMIROMI AND KOCH (Vousoughi et al., 2013) (Table 2). Demographically, according to the national census conducted in 2011, the population of the plain is around 564,000 that inhabits two big cities, including the provincial city Ardabil, and 88 villages (Statistical Center of Iran, 2011). Farming is the major professional occupation in this plain. As a result of massive expansion of Ardabil city, developing irrigated agriculture, and industrial developments, an increasing pressure on the aquifer has been exerted since the 1980s. The annual water consumption of the drinking, industrial, and agricultural sectors amounts to 26, 4, and 177 million cubic meters per year (MCM/yr), respectively; out of which 89% is obtained from groundwater (Ardabil Regional Water Authority (2013). The groundwater supply emanates from 2,622 active pumping wells, 36 qanats, and 77 springs operating across the Ardabil aquifer (K. Kord & Asghari Moghaddam, 2014).

Hydrogeology and Boundary Condition
Quaternary alluvial deposits resulting from the erosion of the Alborz and Sabalan mountains form the Ardabil plain ( Figure 2a) (K. Kord & Asghari Moghaddam, 2014). In the light of geophysical surveys, pumping test analysis, and drilling logs undertaken for the Ardabil aquifer which can provide useful information on the structure of aquifer systems (Jellalia et al., 2015), the Ardabil aquifer is found to be composed of gravel, sand, and a small fraction of clay. The aquifer can be divided into two layers, including a major uppermost unconfined aquifer with a thickness varying from 20 to 40 m, and a tiny confined layer within its central part. However, since these two layers are laterally hydraulically connected, isolating them from one another cannot be a plausible option, especially, when a surface-groundwater interaction investigation is aimed. Most of the pumping wells have been drilled through the unconfined aquifer, with only a small fraction penetrating into the lower layer which lies between a depth of 110 and 220 m (M. Kord et al., 2013). The thickness of the alluvium gradually decreases from east to west. The greatest alluvial thickness of the aquifer,  reaching to 220 m and representing a thick saturated zone, composed mainly of coarse grain sediments, is located in the east and southwest parts of the plain/aquifer. The groundwater heads measured by the piezometric stations, used for monitoring the qualitative and quantitative characteristics of the groundwater in the Ardabil aquifer, are employed to draw the groundwater contour map, as well as to delineate the flow direction ( Figure 2b). The groundwater heads vary from 1,529 to 1,280 MSL (meters above mean sea level) and follow the prevailing topography ( Figure 1c); therefore, the cardinal direction of the groundwater flow is towards the northwest of the plain (Figure 2b). To assign the hydraulic aquifer properties including specific yield (Sy) and transmissivity (T), pumping induced drawdowns at 40 observation wells/piezometers, distributed evenly across the aquifer, were used. Following Todd and Mays (2005), pumping wells having no observation wells were employed for the estimation of the transmissivity (T). The aquifer's T was found to vary widely from 50 to 2,200 m 2 /day, and its Sy from 0.021 to 0.14 (M. Kord et al., 2013).
The boundary conditions of the numerical groundwater flow model (MODFLOW-NWT, Niswonger et al., 2011) were assigned by respecting the geological formations encompassing the aquifer as well as the observed groundwater contour pattern (see Figure 2). To that end, flow/head boundary conditions were specified along most sections of the Ardabil aquifer boundary by means of MODFLOW's FHB-package (Flow Head Boundary), thereby attributing in a general way no-flow, inflow, and outflow boundaries, as illustrated in Figure 2b. Groundwater inflow and outflow into the model domain were calculated using Darcy's law and specified along the northwest boundary of the aquifer where not only all surface water is drained out of the main outlet of the basin (Samian station) ( Figure 1b), but where also groundwater outflow occurs. Some sections of the aquifer boundary, that are in contact with impervious formations, were specified as no-flow, although according to the groundwater contours (see Figure 2b), some inflow may enter the aquifer. Due to the fact that some installed piezometers are quite close to the aquifer boundary, time-variant specified heads (MODFLOW's CHD-package) and recognized as Dirichlet conditions were assigned at places where historical observed groundwater heads were available, for example, P18 ( Figure 2b).

Constructing, Calibrating, and Validating the SWAT-MODFLOW-NWT Coupled Model
The coupled hydrological model employed for this study is the so-called SWAT-MODFLOW-NWT. The model is based on integrating MODFLOW-NWT (Newton-Raphson formulation of the Modular Three-Dimensional Finite-Difference Groundwater Flow model of the eminent MODFLOW model of McDonald et al. (1988)) (Niswonger et al., 2011), with SWAT (Soil and Water Assessment Tool; J.G Arnold et al., 1998), as developed by Bailey et al. (2016). For the purpose of coupling, the SWAT and MODFLOW-NWT models were independently constructed, calibrated, and validated against discharge measured at five hydrometric stations ( Figure 1b) and groundwater heads recorded at 68 observation wells (Figure 2b), respectively. Subsequently, a recalibration was made for the hydraulic conductance of the stream network, where the SWAT and MODFLOW-NWT models' flow domains essentially meet. The finally coupled SWAT-MODFLOW-NWT model was constructed, calibrated, and validated for the 1978-2012 time period wherefore the first 3 years were considered as a spin-up period to allow the water cycle to reach an equilibrium and thus to ensure an initial state (Daggupati et al., 2015). Therefore, the complete study period, that is, 1978-2012, was split into three sub-periods, 1978-1980, 1988-2012, and 1981-1987 for spin-up, calibration, and validation, respectively. The magnitude and direction of the flow between a gaining/losing stream and its adjacent groundwater domain are determined according to the MODFLOW's River package. To that respect, SWAT calculates the stream depth based on Manning's equation, assuming a trapezoidal shape for the reaches/streams, and the SWAT-estimated volume of flow (Neitsch et al., 2011). Then, by adding the computed stream depth to the stream bed elevation extracted from the digital elevation model (DEM), the effective stream stage is calculated to be used in the SWAT-MODFLOW-NWT model (Bailey et al., 2016).
To parameterize the river-aquifer interaction of the coupled model and thus providing a high resolution of the spatial variability of this interaction, the river network was incised in the MODFLOW-grid cells using GIS operations. As a result, 2,581 reaches were created. The magnitude and direction of flow, characterized in gaining and losing reaches, were determined using an adaption of Darcy's law (equation 1).

Water Resources Research
TAIE SEMIROMI AND KOCH where QRIV n (L 3 T −1 ) quntifies the volumetric flow between a reach and its adjacent aquifer (McDonald, Harbaugh,, & Geological Survey (U.S.), 1988) and CRIV n (L 2 T −1 ) is the hydraulic conductance of the reach's bed material as calculated by equation 2.
where K n is the hydraulic conductivity of the reach's bed material with thikness M n , L n is the length of the reach through the grid cell (i, j, k), and W n represents the reach's width.
In integrated hydrological modeling, the river-aquifer interaction is normally parameterized using the river conductance CRIV n , as characterized via a head-dependent flux (Cauchy-type) boundary condition (Cousquer et al., 2017;Furman, 2008;Pascal Goderniaux et al., 2009). As a lump parameter, it cannot be measured at the field scale due to its combination of different parameters/processes. Accordingly, it is calibrated using a history matching approach (Ebel et al., 2009;McDonald et al., 1988;Rushton, 2007), wherefore the simulated variables such as groundwater heads and streamflow discharge are assessed against their observed values. Nevertheless, the calibrated values should be further check out by the soft knowledge of the modeler (Anderson et al., 2015).
In the present study, the hydraulic conductance of the reaches CRIV n was calibrated using a trial-and-error procedure by minimizing the discrepancy between observed and SWAT-MODFLOW-NWT-simulated groundwater/streamflow discharge. Ultimately, the water yield as the final product of the river-aquifer flux exchange, surface runoff, lateral flow, and tile flow was calculated as: where SURQ denotes the surface runoff/direct runoff, LATQ the interflow/lateral flow, TILEQ the tile flow, GWQ the groundwater discharge/baseflow generated from the gaining reaches, and SWGW indicates the concentrated recharge towards the adjacent aquifer produced in the losing streams (Park & Bailey, 2017).
It should be noted that in the application area, that is, the GRB, there is no tile drainage system. Therefore, the term TILEQ in equation 3 was "zero" in this study.
The details of the individual calibrations and validations of each model (SWAT, MODFLOW-NWT) and of the coupled one are described by Taie Semiromi and Koch (2019). In principle, this integrated hydrological model affords a dynamic variable exchange between SWAT and MODFLOW-NWT. To that respect, first, the deep percolation/recharge estimated by SWAT is passed to MODFLOW-NWT; second, the river heads calculated by SWAT are fed to the River package (RIV) of MODFLOW-NWT; and third, the computed groundwater discharge/baseflow of the gaining reaches is delivered to SWAT and, similarly, the concentrated groundwater recharge of the losing reaches computed by the River package (RIV) of MODFLOW-NWT is subtracted from the total streamflow (direct surface runoff + baseflow + lateral flow). A flowchart of the exchanged variables of the proposed coupled modeling as described is sketched in Figure 3a; with the details of the various methodological steps provided by Bailey et al. (2016) and Taie Semiromi and Koch (2019).

Downscaling of Hydroclimate Variables
In the present study, daily minimum and maximum temperature at five weather stations distributed across the GRB (Figure 1b) were downscaled using the well-known Statistical DownScaling Model (SDSM) (Wilby et al., 2002) which has found innumerable applications in climate change studies of the last two decades (see Wilby & Dawson, 2013, for a review). The daily historical data records available at these stations representing the reference period fall within the 1978-2005 time period. Climate projections were available for RCP (Representative Concentration Pathways) 2.6, RCP 4.5, and RCP 8.5 up to year 2100. These RCPs represent the range of GHG emissions reported in the wider literature well; they include a stringent mitigation scenario (RCP 2.6), two intermediate scenarios (RCP 4.5 and RCP 6.0), and one scenario with very high GHG emissions (RCP 8.5). After screening the predictors, based on the strongest correlations between various climate predictors (e.g., near-surface temperature, and mean sea level pressure) and observed weather variables (min and max temperatures), depending on the weather station between two and three predictors The daily precipitation was downscaled for 12 weather stations ( Figure 1b) using a state-of-the-art hybrid statistical downscaling model that takes advantages of artificial neural networks (ANNs), discrete wavelet transform (DWT), and quantile mapping (QM) (Taie Semiromi & Koch, 2017). Using a stepwise regression procedure statistical relationships between the National Centers for Environmental Prediction (NCEP)-predictors and the daily precipitation at each station for the reference period 1978-2005 were established and thus the most suitable predictor(s) were screened, similar to above. Daily predictors of CanESM2 under RCPs 2.6, 4.5, and 8.5 were then employed to force the developed hybrid model for the future projection period 2006-2100. The coupled wavelet and neural network (DWT-ANN) model employs various subseries components of the DWT-multiresolution analysis of Level 4 of the selected climate predictors at a weather station. These decomposed predictors were used as inputs to train (calibrate) the model with respect to the observed precipitation (predictand) as output. Afterwards, the DWT-ANN-simulated daily precipitation was bias-corrected using the QM algorithm.

Projection of the Effects of Climate Change Scenarios and Pumping Variants on the Various Water Resources Compartments
Once the daily maximum and minimum temperature and precipitation data were successfully downscaled, the coupled hydrological model, that is, SWAT-MODFLOW-NWT, was forced by the downscaled hydroclimate variables under six scenarios/variants (three × two) including three RCPs 2.6, 4.5, and 8.5, with each of them applied to the coupled model under two groundwater utilization varinats, namely, "pumping" and "nonpumping". By doing so, one can differentiate the net impacts of climate change from that of Water Resources Research groundwater overutilization. Given the fact that the Ardabil aquifer has been strictly regulated to not be further exhausted by new pumping wells (Ardabil Regional Water Authority, 2013), the number of pumping wells and therefore the amount of groundwater extracted for the future projection period for the "pumping" variant were assumed to be as much as the last year of the hydrologic model simulation (2012), namely 2,622 wells, extracting a total of 160 MCM per year (K. Kord & Asghari Moghaddam, 2014).This pumping rate corresponds to an equivalent water height across the total area of the Ardabil aquifer of about 150 mm/yr, which is nearly 15 times higher than the total diffuse surface-to-groundwater recharge of 10 mm/yr during the historic period, demonstrating clearly the extreme overutilization of this aquifer. For the "nonpumping" variant on the other hand, all these wells were shut off in the coupled hydrological model during the future simulation periods, in order to be able to separate out the hydrological effects due to groundwater pumping from those due to climate change. A similar approach was taken by Scibek et al. (2007)  Finally, the spatiotemporal variability of the various water resources components of the GRB was extracted for three future time slices: 2013-2040 (near future), 2041-2070 (middle future), and 2071-2099 (distant future) and under the six aforementioned climate/pumping condition combinations. These water resources components include meteoric infiltration (diffuse recharge), baseflow emanating from the gaining streams, concentrated groundwater recharge contributed by the losing streams, baseflow index (amount of baseflow divided by amount of total streamflow), groundwater storage, and water yield as calculated by equation 3. It should be noted that since the coupled hydrological model produces the output variables at a daily temporal scale, they were then aggregated and averaged to produce monthly outputs.

Downscaling of Min and Max Temperature
The results of the SDSM-downscaled min temperature show (see Figure 4) that, in comparison with that of the reference period , a notable increase in min temperature, particularly, during winter (e.g., January and February), will occur under all three emission scenarios (RCP 2.6, 4.5, and 8.5), wherefore (and not unexpectedly) this increase is more pronounced for RCP 8.5 and the distant future (2071-2100). Nevertheless, a minor decrease in min temperature is to be expected in September for all stations except for Sareyn, for which the coldest days records on a national scale have been reported several times (Ardabil Regional Water Authority, 2013 Similarly, as illustrated by Figure 5, the downscaled future max temperature will be considerably raised at all stations and months, except for April. Similar to the downscaled min temperature, the max temperature will be dramatically increased when going from RCP 2.6 to 8.

Water Resources Research
and 2-4°C in autumn. Moreover, they found that the minimum and maximum temperature will occur in months January and July, respectively. Overall, their results are fairly consistent with the projections made in the present study for the distant future (2071-2100) under RCP 4.5.

Downscaling of Precipitation
The daily precipitation projected by the DWT-ANN-QM model for the three future time slices , and 2071-2100 under the three emission scenarios RCPs 2.6, 4.5, and 8.5, are illustrated in Figure 6. These projections manifest that, compared with the reference period , the future

10.1029/2019WR025388
Water Resources Research precipitation will decrease under all three scenarios, particularly for the wet season, that is, from January to May for the climate stations Shams Abad, Namin, while less reductions at station Ardebil under all three RCPs, and even a precipitation increase is to be expected during the dry period, that is, June to September at some stations, such as Hir and Shams Abad. This projected rise in precipitation in the future wet seasons over the GRB has also been corroborated by Javan et al. (2015) under the SRES scenario B2 as their findings show that the precipitation will be increased for the future period (2071-2100) for the months January, February, March, September, and December. Moreover, our results indicate further that the reduction in the monthly precipitation under the extreme RCP 8.5 is surprisingly less than that projected under RCP 4.5. This peculiar behavior can be associated with a higher perceptible water supply of the atmosphere resulting from a higher temperature increase under RCP 8.5 than RCP 4.5 (Huang et al., 2014). Overall, the precipitation will decrease in the near future (2013-2040), then it will rise in the middle future (2041-2070), and again decrease in the distant future (2071-2100). In Table 3, the estimated average of the annual precipitation for the reference period  and three future periods (2006-2040, 2041-2070, and 2071-2100) under RCPs 2.6, 4.5, and 8.5 is listed. The numbers in the table corroborate Figure 6; that is, for most of the weather stations, the average annual precipitation will decrease for the three future time periods, relative to the historic one. As indicated by a " * " in A climate study conducted by Terink et al. (2013), using nine GCMs and SRES A1B, 2013 demonstrated that the annual precipitation will decline by 15-20% for the majority of countries across the Middle East and Northern Africa for the two future periods 2020-2030 and 2040-2050, which appears to confirm the reduction witnessed in the present study, namely, a 12% decline in the precipitation over the GRB, on average. Nonetheless, in some studies, by contrast, a rise in the future precipitation has been reported. For instance, Modarres et al. (2018) projected the variability of the maximum daily rainfall across the North and Northwest of Iran for 2020-2049 using six GCMs under the two scenarios A2 and B1. In particular, their results suggest the projected precipitation to vary from −16 to 7% under scenario A2, compared with the reference period . Likewise, depending on the GCM used, the future maximum rainfall may change from −16 to 10% under scenario B1, relative to the reference period.

Projection of Diffuse Recharge Under the Various Climate Change Scenarios
The SWAT-MODFLOW-NWT-estimated monthly recharge (deep percolation) over the GRB for the reference simulation period , the projected near future (2013-2040), middle future (2040-2070), and distant future (2071-2099) periods under RCPs 2.6, 4.5, and 8.5 are represented in Figure 7. One can notice that the recharge will be increased from months of January to April (mostly winter season) for all three future time periods under RCPs 2.6 and 4.5, while such an increase is only expected for January, February, and April for the middle and distant futures under RCP 8.5. However, in the rest of the months, the recharge will be nearly the same as that of the reference simulation period (e.g., July to October, except for the distant future period 2071-2099) or will decline (e.g., March and May under RCP 8.5). In addition to these future changes of the magnitude of the recharge, Figure 7 shows further that its timing will most likely also be changed, due to an anticipated winter warming (Woldeamlak et al., 2007). Thus, one can notice from Fgure 7 that, whereas the maximum recharge occurs in March during the reference simulation period, it will be retarded and moved to April, and this holds for all three future time periods and the three RCPs.
It is interesting to note that this one-month retardation of the future maximum recharge occurs despite the fact the projected precipitation (see Figure 6), will decrease in April, irrespective of the future time slices and emission scenarios considered. However, as a result of the notable projected increase in the min and max temperature for the future winters (see Figures 4 and 5), and due to the fact that a major fraction of the precipitation is received and stored as snow, particularly in the mountainous region of the basin (Ardabil Regional Water Authority, 2013), an earlier hillslope snowmelt runoff, and correspondingly, earlier recharge may be expected (Hayashi & Farrow, 2014). In this study, in contrast, the maximum recharge is projected to be shifted 1 month later, namely, April. This happens because a major portion of the precipitation, as the main driver of the hydrological cycle (e.g., J. G. Arnold et al., 1993), is evapotranspirated during the winter seasons, owing to the considerable projected rise in the min and max temperature. However, as Figures 4

Water Resources Research
and 5 indicate, the temperature will not be raised in the future months of April, so that more water for the diffuse recharge will be available during that month.

Projection of Surface Runoff Under the Various Climate Change Scenarios and Pumping Variants
The SWAT-MODFLOW-NWT-predicted climate change impacts on direct runoff/surface runoff shown in Figure 8 demonstrate that for all three future periods and all three emission scenarios RCP 2.6, 4.5, and 8.5, this hydrological quantity will considerably decline (Figure 8), compared to that of the reference simulation period  and this hold for both ("pumping" and "nonpumping") variants. As illustrated in Figure 8, this reduction in the surface runoff is more (expectedly) pronounced when going from RCP 2.6 to 8.5 as well as from the near future (2013-2040) to the distant future (2071-2099). However, the pumping variants exert no effect on the magnitude and timing of the surface runoff. Interestingly, similar to the projected recharge above (see Figure 7), the maximum projected surface runoff will shift from March to April, for the same reasons as discussed earlier. Indeed, surface runoff is partitioned from the same driver (effective precipitation) of the hydrological cycle, with the latter predicted to be majorly reduced during the future winter seasons owing to a projected large rise in the min and max temperature during that period. And as a nonsignificant change in the temperature projected for April, less evapotranspiration is expected, so that, in turn, a larger fraction of the precipitationwhich was shown to decrease under all emission scenarios and all time slices in April (see Figure 6)may partition off into the runoff during that month. This peculiar behavior may be ascribed to some counterbalance of the nonsignificant changes in the min and max temperature and the projected reductions in the precipitation. In addition, given the fact that, on average, 24%, 25%, 37%, and 67% of the total precipitation is received as snow at Ardebil, Samian, Nir, and Lay stations, respectively (Ardabil Regional Water Authority, 2013), and owing to an increase in the precipitation

Water Resources Research
projected for the winter for most time slices and emission scenarios, the stored snow will start melting in warmer months (e.g., April). Accordingly, this can further increase the occurrence of the maximum surface runoff, witnessed for April. These results are to be contrasted with the statements of the review of Huber et al. (2005) and of the Canadian watershed study of Scibek et al. (2007) who found that an earlier snowmelt and, correspondingly, an advanced hydrograph peak generally arises in mountainous regions in the wake of global climate change. Unlike the findings of this study, Scibek et al. (2007) also discovered that a noticeable increase in winter discharge will occur in future climate scenarios, most likely due to an increase of rain and snowmelt volumes during the winter season under a warmer climate.

Projection of the Gaining Streams' Reactions to the Various Climatic Scenarios and Pumping Variants
Due to the fact that the coupled hydrological model SWAT-MODFLOW-NWT has greatly enhanced the estimation of the baseflow/groundwater discharge and since the baseflow generally sustains the river flow, particularly during dry periods in semiarid regions like the GRB, an understanding of the impacts of climate change on this component of the total streamflow is of paramount importance. As baseflow is essentially related to groundwater flow, which, in turn, is strongly affected by external stresses, like pumping, the net impact of the latter on baseflow needs to be investigated and separated from that of climate change. Thus, in the following, the baseflow is analyzed for the two variants of SWAT-MODFLOW-NWT model, that is, "pumping" and "nonpumping." The results, shown in Figures 9a-9c, indicate that for the "nonpumping" model variant, for all future time slices and all of emission scenarios the baseflow will be substantially increased, irrespective of the reduction and increase obtained earlier for the precipitation and the min and max temperature, respectively. However, the same is not true under the "pumping" model variant, since in that case the baseflow will markedly decrease under all emission scenarios and time slices, except for the "near future" when such a baseflow decline cannot be detected. Such a great sensitivity of the baseflow to (1) Since the modeled area/Ardabil aquifer is a plain (see Figure 1c) where the groundwater level is close to the land surface, even a slight change of pumping (and less so of climatic) stress can exert a big influence on the groundwater heads (Melo & Wendland, 2017); (2) the groundwater draft from the Ardabil aquifer has dramatically expanded from 35 MCM in 1978 to 160 MCM in 2012, (corresponding to a water height of 150 mm/yr, that is, about 15 times the present-day total diffuse annual groundwater recharge of 10 mm/yr), thus reversing the current groundwater utilization into the "nonpumping" variant (describing an undisturbed system in which the utilization is satisfied by natural springs only) which will lead to a rise of the groundwater heads and, as a result, to an increased discharge contribution from the groundwater towards the gaining streams, albeit a reduction (by-12% on average) in the precipitation and an increase in the min and max temperature.
In Figure 10, the monthly time series of the historic and projected baseflow under the three emission and the two pumping variants are provided, in addition. First of all, one can notice from the three panels of the different RCPs (2.6, 4.5, 8.5) that the differences in the baseflow time series are minor, which indicates that the impacts of the emission scenarios on the hydrological quantity are more or less also negligible. However, there are systematic differences between the two pumping variants. For the "nonpumping" variant, one can notice a discontinuous jump of the baseflow in year 2013, showing the immediate reaction of the system, once the pumping stress is shut-off. From that time onwards, there is monotonous increase of the baseflow, up to about year 2030, and a steady decline thereafter until the end of the future projection period (year 2100), demonstrating the "kick in" of long term climate change. This general course of the baseflow holds also for the "pumping" variantassuming the current groundwater utilization volume of 160 MCM up to year 2100 -, though its runs with a more or less parallel decrement from its "nonpumping" counterpart. However, under this "pumping" variant, climate change appears to have a measurable effect on the baseflow towards the gainng streams only from about year 2045 onward. All this is a consequence of the fact that the groundwater heads will fall on account of the balance lost between the total groundwater recharge and discharge. Accordingly, owing to Darcy's law implemented in the MODFLOW's RIV package (equation 1), the concentrated groundwater discharge contribution towards the gaining streams will be reduced in the wake of a decline of diffuse groundwater recharge, ensued, in turn, from the reduction of the projected have broadly linked climate change-induced precipitation changes to spatiotemporal variabilities of groundwater recharge. However, whereas there is an unequivocal rise in temperature due to the global warming, the projected recharge could increase or reduce, depending on how precipitation reacts to climate change effects (Jyrkama & Sykes, 2007).
To characterize in more detail how the baseflow, relative to the surface runoff and lateral flow, will react to the climate change impacts, we calculated the so-called baseflow index BI = baseflow/(surface runoff + lateral flow + baseflow) for the two pumping variants, using these SWAT-MODFLOW-NWT-simulated variables for the historic and projected periods. The results shown in Figure 11 indicate that for both pumping variants, under all three emission scenarios, the future baseflow index will be considerably increased for all months, except for April. Moreover, similar to the previously projected timings of the recharge and surface runoff (1 month delay in the maximum surface runoff and recharge, relative to the reference simulation period) and, unlike the findings of Scibek et al. (2007) who revealed a 1 month earlier peak discharge, the timing of the maximum baseflow index is also projected here to be shifted by one month, Figure 10. SWAT-MODFLOW-NWT-simulated historic and projected monthly time series of groundwater discharge (baseflow) to the river network (gaining streams) under the three emission scenarios RCP 2.6 (a), 4.5 (b), and 8.5 (c), and and the two pumping variants ("pumping" and "nonpumping"). Note that the noticeable breaks in 1988 of the historic period result from a a three-year spin up period (an adjustment process for the model to reach an "optimal" state, where internal stores such as groundwater storage move from the estimated initial condition to an "optimal" state (Kim et al., 2018)) considered before running the calibration and validation periods of the coupled model.

Water Resources Research
namely from August in the historical period to September in the future projection period (see Figure 4). This shift may be attributed to the projected maximum monthly max temperature witnessed in future August (see Figure 5), but reduced min temperature in future September (see Figure 4). This spectacular rise of the projected baseflow index, in comparison with that of the reference simulation period, can be associated with the noticeable decline of the projected for the direct surface runoff (see Figure 8) which makes up the substantial fraction of the total streamflow in comparison to the baseflow and lateral flow. As such, due to the fact that the surface runoff appears in the denominator of the BI-definition index, a smaller surface runoff contribution, projected under the different climate change scenarios here, will deliver a bigger baseflow index. Furthermore, since the baseflow appears in both numerator and denominator of the BI formula and, owing to the fact, that the surface runoff and lateral flow are basically not influenced by either of the two pumping variants (see Figure 7), the baseflow index, in turn, will not be affected neither by the latter, and this is exactly what is seen in Figure 11.

Projection of the Losing Streams' Reactions to the Climatic Scenarios and Pumping Variants
Similar to the gaining streams in the previous section, the reactions of the losing streams to the different climate change scenarios and pumping variants are shown in Figure 12 by the time series of the concentrated groundwater recharge, produced under influent conditions from the losing streams. These graphs indicate that there is a noticeable difference of this quantity between the historic and the projected period, particularly under the emission scenarios RCP 4.5 and 8.5. However, compared with the baseflow of the gaining streams of the previous section, the concentrated recharge of the losing streams here is (1) more than an order of magnitude smaller and (2) less affected by the two pumping variants, that is, "nonpumping" and "pumping. This can be attributed to the fact that the depth to the water table is very shallow in most parts of the Ardabil aquifer, thus allowing for more groundwater contribution towards the river network (effluent conditions), so that influent conditions from losing streams rarely occur (Taie Semiromi & Koch, 2019). Such influent conditions take place mainly across the southern part of the aquifer at locations where a considerable depth to water table has been established due to groundwater overutilization (Figure 2b) (Taie Semiromi & Koch, 2019). Moreover, as discussed earlier, the surface runoff will decrease under the different climate scenarios considered (see Figure 8), which indicate less water will be available to be recharged over

10.1029/2019WR025388
Water Resources Research the banks and bed of the losing streams to the groundwater. Surprisingly, in contrast to the projected rise of the baseflow under the "nonpumping" variant in the previous section (see Figure 10), Figure 12 shows that the concentrated groundwater recharge from the losing streams will slightly drop in the future under this "nonpumping" variant. This can be ascribed to an overall future rise of the groundwater table, expected in this case. Under this circumstance, as the groundwater heads rise, according to the river-groundwater interaction formulation of Darcy's law (equation 1), the direction of the exchange flow will be reversed, so that some losing streams become gaining ones. In addition to these temporal assessments of the river-groundwater interactions above, the spatial variability of the losing and gaining streams across the study region will be analyzed in the subsequent section.

Projection of the Spatial Variability of the Gaining and Losing Streams Under the Different Climatic Scenarios and Pumping Variants
One of the distinct advantages of the current SWAT-MODFLOW-NWT coupled hydrological model over previously developed separate surface-and groundwater flow models is that the former is capable to produce high-resolution outputs of the surface-groundwater flux exchanges. According to the flow directions between the 2,581 reaches/streams and the adjacent aquifer, obtained by overlaying the MODFLOW-grids of the groundwater domain and the river network (Taie Semiromi & Koch, 2019) (Figure 1d), we identified 922 and 1,659 gaining and losing streams, respectively. However, the results of the SWAT-MODFLOW-NWT-spatially variable projections shown in Figures 13-15 for the three climate scenarios, RCP 2,6, 4.5 and 8.5, respectively, illustrate that the number of the gaining and losing streams will vary widely in the future, depending on the climatic scenarios and the pumping variants considered. Thus, for RCP 2.6 ( Figure 13), one can notice that under the "pumping" variant, the maximum daily

10.1029/2019WR025388
Water Resources Research baseflow yielded by the gaining streams (effluent conditions) will be increased from 622.71 m 3 /day for the historic period (red bars) to 709.29 m 3 /day for the time slice 2041-2070 (the period for which the precipitation is expected to rise) (green bars). Meanwhile, for the "pumping" variant, although there is still an increase of the future maximum daily baseflow, the latter is only about half of that for the "nonpumping" variant (from 622.71 to 673.98), which indicates that pumping reduces clearly the discharge of the groundwater aquifer to the gaining streams. The barplots of Figure 16 condense the numbers of the spatial variability maps above by quantifying the total groundwater discharge/recharge to/from the river network for the three emission scenarios and the two pumping variants. For example, the panel of Figure 16b shows that, in spite of the aforementioned increases of the future maximum daily baseflow, the total amount of daily groundwater discharge (baseflow) across the study region will slightly and continuously drop for the future time slices for the "pumping" variant, irrespective of the RCP. Nonetheless, the same does not hold for the "nonpumping" variant ( Figure 16a) where the projected total average daily baseflow will grow from its historic amount of 63416 m 3 /day to values higher than 70000 m 3 /day depending on the time slice, irrespective of the RCP. This increase witnessed for the future maximum daily baseflow (groundwater discharge) under the projected climate change in the study region under the "nonpumping" variant scenarios is most likely a result of the more frequent extreme precipitation events that are to be expected in the wake of the global warming predicted over most mid-latitude land masses and across wet tropical regions (a global-scale intensification of heavy precipitation) (IPCC, 2013). This will lead, consequently, to more diffuse surface-to-groundwater recharge, giving rise, in turn, to a shallower groundwater table, thus intensifying the groundwater discharge contribution (baseflow) towards the river network (gaining streams). This is particularly true for the Ardabil aquifer, where the depth to the water table is already very shallow (Ardabil Regional Water Authority, 2013), so that the groundwater residence time of discharge to the streams will be quite short (Modica & Buxton, 1998).

Figure 16.
Historic and future (for the three time slices as indicated) projected daily total groundwater discharge (top panels) and recharge (bottom panels) to/from the river network under the three emission scenarios (RCPs) for "nonpumping" (left panels) and "pumping" (right panels) variants.

Water Resources Research
Similar to the gaining streams above, the spatial distribution of the concentrated maximum groundwater recharge (influent condition) of the losing streams taking place over the bed and bank of the streams is illustrated in Figures 13-15 by the blue and brown bars for the historic and future period, respectively. One can see from the size of these bars, that neither the climatic scenario nor the pumping variant exerts a notable influence, unlike of what has been observed for the baseflow of the gaining streams. On the other hand, as shown by the two lower panels of Figure 16, the total daily volume of the concentrated groundwater recharge from the losing streams over the whole aquifer will, however, significantly differ for the three emission scenarios and the two pumping variants. Thus, whereas for the "nonpumping" variant, the future total stream-groundwater recharge is more or less the same as that of the historic period,though somewhat diminished for the most extreme RPC 8.5 climate scenario -, the opposite is true for the "pumping" variant, where this quantity increases by up to~20%, namely for the distant future time slice (2071-2100). Nevertheless, in spite of a much higher number of the losing streams, that is, 1,659 rather than 922 for the gaining streams above, the concentrated groundwater recharge from the losing streams is about three orders of magnitude smaller than the groundwater discharge to the gaining streams (baseflow). This is due to the fact that both the hydraulic conductance and the head differences between the losing streams and the adjacent groundwater aquifer are generally rather small (Irvine et al., 2015), namely in the eastern, flatter topography of the GRB. Thus, the maximum daily concentrated groundwater recharge from the losing streams of about 0.4 m 3 /day, projected under RCP 2.6 and "pumping" for the middle future (2041-2070), takes place mainly across the western part of the aquifer (see Figure 13d).
According to the precipitation rise projected under RCP 2.6 for the time slice 2041-2070 (see Figure 6), a growth in the hillslope surface runoff and, thus, more water to the river network is to be expected. But, on the other hand, the impacts of the higher precipitation on the runoff generation for the scenario and time slice considered may be masked by the unequivocal temperature rise projected for the GRB (Figures 4 and 5). This may be the reason that no noticeable changes can be spotted for the concentrated groundwater recharge from the losing streams, with the projected maximum concentrated groundwater recharge differing only by 0.35 to 0.40 m 3 /day (see Figures 15a and 13d). Under the "pumping" variant, as an increased groundwater head drawdown is to be expected (Figure 17), the head differences between the losing river's stage and the groundwater table will be increased, so that, the concentrated groundwater recharge will grow. For that reason, the latter will occur under the "pumping" variant. Likewise, the projected daily average of the concentrated groundwater recharge volume over the whole aquifer, in relation to that of the historic period, is also anticipated to increase under the "pumping" variant, regardless of the climatic scenarios considered ( Figure 16d). As expected, the concentrated groundwater recharge from the losing streams will fall noticeably under the "nonpumping" variant, and its minimum will take place under RCP 8.5 (Figure 16c). On the other hand, it may be possible that some of the losing streams will convert into gaining streams under the "nonpumping" variant, when a rise in the groundwater heads may occur again. In conclusion of this section, it has been shown clearly that the groundwater contribution, that is, the baseflow towards the gaining streams is to be recognized as the dominant surface-groundwater interaction process for the GRB also in the future, extending the corresponding findings of Taie Semiromi and Koch (2019) for the historic period.

Projection of Groundwater Storage Under the Various Climatic Scenarios and Pumping Variants
In Figure 17 the climate change impacts on the temporal evolution of the groundwater storage of the aquifer are shown for the different RCPs and the two pumping variants. These time series indicate that under all of three climatic scenarios and for both the "nonpumping" and "pumping" variants, the groundwater storage will monotonously decrease in the future though with a gentler slope than what is observed for the historic period. Surprisingly, even for the "nonpumping" variant, one can see an ongoing reduction of the groundwater storage in the future, and this occurs, despite the fact that groundwater depletion that has been occurring on account of a dramatic overutilization in the recent past (Ardabil Regional Water Authority, 2013; Taie Semiromi & Koch, 2019), should have been reversed to a large extent under this "nonpumping" variant. Thus, by comparing the curves for the two pumping variants in Figure 17, one can conclude that the unsustainable groundwater utilization in the GRB is the determining factor for the observed decrease of the groundwater storage there for the "nonpumping" variant, exacerbated, moreover, by climate change impacts that add even more pressure on the system. Therefore, for this "nonpumping" variant, one can

10.1029/2019WR025388
Water Resources Research clearly notice that even by assuming an undisturbed system ensured by eliminating the abstraction performed by the 2,622 farming pumping wells extracting around 150 MCM water per year, the initial condition of the first year of the historic simulation period, that is, 1980, cannot be recovered even at the end of the projection period in year 2099. In accordance with the findings of this study, Kirby et al. (2016) affirmed that the future planned irrigation developments in five main regions of Bangladesh can have a larger impact on the water balance components than the projected climate change itself. This is also the case in the present study here, as one cannot distinguish any difference in the groundwater storages projected under the three climatic emission scenarios (RCP 2.6, 4.5 and 8.5) depicted in the three panels of Figure 17. This insensitivity may be explained by the only minor contribution of the diffuse recharge (meteoric infiltration), amounting to 17.8 MCM (equal to only 11% of the total water budget) (Ardabil Regional Water Authority, 2013; Taie Semiromi & Koch, 2019). In addition, in the wake of the earlier projected precipitation decreases (see Figure 6) for the three RCPs considered, the future diffuse surface-to-groundwater recharge will not only become smaller, but also induce more or less equal groundwater storages for these three climate emission scenarios.

Projection of Water Yield Under the Various Climatic Scenarios and Pumping Variants
Ultimately, we calculated the future total monthly water yield (defined in equation 3) of the GRB as the final product of all surface-and subsurface hydrological contributions, including their interactions, extending the work of Taie Semiromi and Koch (2019) who estimated the water yield in the GRB only for the historic reference period. In this way, the impacts of the different climate change scenarios and of the pumping variants on the future water yield can be assessed. Figure 18 illustrates that under all of three RCPs and for both

Water Resources Research
pumping ("nonpumping"and "pumping") variants the future monthly water yield is anticipated to fall beneath that of the historic period except for month of April in the time slice 2041-2070 under the (rather benevolent) RCP 2.6, for which an increase in the precipitation has been found earlier for half of the stations in the GRB (see Table 3). Similar to the future projected surface runoff and recharge (see Figures  8 and 9, respectively), the timing of the maximum water yield will, in relation to the historic period, also be shifted 1 month later, owing to the fact that the major fraction of the water yield emanates from the surface runoff which experiences the same fate. It can be further noticed from the various panels of Figure 18, and even more clearly from the time series of the monthly water yield depicted in Figure 19 for the different RCPs / pumping variants that, as expected, the magnitude of the average monthly water yield will generally decline when going from the rather benevolent RCP 2.6 to the most extreme RCP 8.5 climate scenario. Strangely, no future increase of the average water yield under the "nonpumping" variant can be detected, although some sporadic and irregular bursts of the water yield can take place ( Figure 19). To explain this kind of peculiar behavior, one should be reminded that, even though an appreciable rise of the baseflow, as one of the contributors to the water yield, has been found earlier under the "nonpumping" variant (see Figures 9 and 10), baseflow generally constitutes only an insignificant fraction of the total streamflow in the GRB. This is because the flow regime of the GRB represents an intermittent/ephemeral system, particularly over the last decade (2000-2012) (Taie Semiromi & Koch, 2019), when baseflow contributes to streamflow only in a sporadic and irregular way. In other words, GRB's river system characterized originally by a perennial system has been converting into an intermittent system due to a depletion of the groundwater aquifer with a subsequent drop of the groundwater table, that is, decreasing head differences with the adjacent gaining streams. Furthermore, and more importantly, as surface runoff (but not lateral flow), as the main constituent of the water yield (see equation 3) was shown earlier (see Figure 8) to be not affected by any of the two pumping variants, it is of no surprise that the future water yield will not be impacted neither. Nonetheless, in conclusion of this section, as the surface runoff will generally diminish in the wake of the climate change projected for the various RCP emission scenarios (see Figure 8), less water will be yielded in the future than in the historic period as well. And, counterintuitively, this future adverse situation cannot be remedied neither

Water Resources Research
by shutting off the groundwater pumping wells, as this was shown to have no effect on the future projected (decreasing) surface runoff (see Figure 8) which, in turn, will lead to the diminished future water yield seen in Figures 18 and 19.

Summary and Conclusions
Most water resources studies carried out to date have been directed towards predicting the potential impacts of future climate change on surface water hydrology, while less attention has been paid to such impacts on groundwater hydrology. And even less as of yet are researches geared towards a better understanding of the net effects of the climate change on the fate of gaining and losing streams in a watershed/ aquifer, as this requires a coupled, more complex approach to the intertwined surface-and subsurface hydrological compartments. In the light of this scarcity of appropriate studies regarding this topic, the present study aims, therefore, to shed light on the different processes acting in both the surface-and subsurface water compartment, including, in particular the complex hydrological interactions of the surface and groundwater domain, as they are mainly reflected by the processes taking place in the losing and gaining streams. This task is achieved through the application of the fully coupled surface-subsurface hydrological model SWAT-MODFLOW-NWT to the Gharehsoo River Basin (GRB), with the enclosed Ardabil aquifer, located in northwestern Iran. This coupled hydrological model was forced with the downscaled climate predictors (scenarios) under three representative concentration pathways (RCPs), namely, RCP 2.6, 4.5, and 8.5 for up to year 2100, dividing the whole simulation period in a reference historic period  and three future time slices 2013-2040, 2041-2070, and 2071-2100. Furthermore, to distinguish the influence of the prevalent groundwater overdraft in the GRB on the gaining and losing streams from that of climate change, we also run the climate-forced coupled hydrological model under two future pumping variants, that is, "nonpumping" (assuming a shut-off of all groundwater Figure 19. Similar to Figure 17, but for the monthly total water yield of the GRB.

10.1029/2019WR025388
Water Resources Research extraction at the end of the historic period) and "pumping" (assuming an ongoing of the average groundwater extraction of the recent past also in the future).
First of all, the results of temperature downscaling indicate that a significant rise in monthly min temperature, particularly during winter (e.g., January and February) under all scenarios, wherefore this increase will become more tangible for RCP 8.5 and the distant future (2071-2100) time slice. In contrast, for the future months of September a minor decrease in min temperature is observed for all stations, except Sareyn.
Similarly, the prediction of monthly maximum temperature demonstrates that a considerable future increase for all months, except April at all weather stations. Similar to the downscaled min temperature, the max temperature will also soar dramatically when going from RCP 2.6 to 8. The projections of the precipitation using a recently developed hybrid model DWT-ANN-QM of the authors (Taie Semiromi & Koch, 2017) point out that the daily precipitation under RCP 2.6, 4.5, and 8.5 will, relative to the reference period , slightly decrease for the 2006-2040 time slice, slightly increase for 2041-2070, and decrease again for 2071-2100, resulting overall in an average decrease of about 12% across the GRB by year 2100. Interestingly, this percentage reduction of the precipitation will be more under RCP 4.5 than under RCP 8.5, owing most likely to the fact that the higher temperatures anticipated for the later than the former will bring about higher evapotranspiration, that is, more water vapor into the atmosphere and, thus, will be higher propensity for precipitation under RCP 8.5. On the other hand, for all RCPs, the monthly precipitation will increase in future months of the wet period December-February, meanwhile for all the other months a decrease, except under RCP 8.5 and somehow 4.5, is expected. Using these downscaled climate predictors in the SWAT-MODFLOW-NWT coupled hydrological model, the various water budgets components of the surface-subsurface compartment of the hydrological cycle (surface runoff, diffuse surface-to-groundwater recharge, baseflow, etc.) were projected into the future.
For example, regarding the diffuse groundwater recharge, due to the projected rise in the future precipitation during the winter, an increase of this quantity, relative to the historic period, will also occur across all future periods considered, 2013-2040, 2040-2070, and 2071-2099 and this holds for all three RCPs (2.6, 4.5 and 8.5).
In the other dry season months (e.g., July to October), the diffuse groundwater recharge will be the same as that of the historic period (except for distant future period 2071-2099), because the soil water content will be very low during that season, so that even in case of a storm event, most of the infiltrated water will be trapped due to a very high matric potential of the soil; hence, inhibiting diffuse groundwater recharge. However, the latter is expected to decrease in some future months, such as March under RCP 8.5 and May under all three RCPs. In addition to the change of the magnitude of the recharge, its peak will also be shifted from March to April owing to the anticipated future winter warming. The SWAT-MODFLOW-NWT coupled surface-groundwater model analysis of the GRB of Taie Semiromi and Koch (2019) indicates that the major source of the groundwater recharge for the Ardabil aquifer in the recent past is the concentrated recharge generated from the losing streams, corroborating previous qualitative estimations of Ardabil Regional Water Authority (2013). In contrast, the rate of the present-day meteoric groundwater recharge (diffuse recharge) is with only~12 mm/yr (see Figure 7 and Taie Semiromi and Koch (2019), for details) rather low, due to (1) the prevailing semiarid climate with onlỹ 270/yr mm precipitation which leads to a very high soil matric potential in most of the seasons, so that a major fraction of the total infiltration (gross recharge) is retained in the vadose zone, in turn, leaves less water for the deeper net groundwater recharge; (2) the build-up, impervious surface area of Ardabil city, as the capital of the same province, covers a noticeable fraction of the aquifer's total area (around 70 km 2 ), so that diffuse surface-to-groundwater recharge is almost negligible. To quantify the range of uncertainty of the estimated groundwater recharge rate for this aquifer, Healy (2010) suggested that this hydrological variable should be evaluated using different approaches. Although, on average, 89% of the total water consumption in the GRB is supplied from groundwater resources (K. Kord & Asghari Moghaddam, 2014), knowledge of the Ardabil aquifer's hydrogeological characteristics, including recharge rate, suffers from a 10.1029/2019WR025388 Water Resources Research lack of detailed exploration. Thus, our study here can be considered a supporting research task in that regard that can enhance the understanding of the various hydrological processes controlling the interaction of the surface-and subsurface water compartments in the GRB.
The SWAT-MODFLOW-NWT-predictions of the future climate change impacts on the surface-groundwater interactions demonstrates that the baseflow generated from the groundwater discharge to the gaining streams will be more influenced than the concentrated groundwater recharge from the losing streams. However, more importantly, the future fate of the baseflow depends tremendously on the pumping variant used in the model. Thus, whereas for the "nonpumping" variant (shutting off all wells after year 2013) the baseflow will significantly increase in the futurethough tapered noticeably by adverse effects of climate change (less precipitation) in the distant future time period 2071-2099-, the opposite is true for the "pumping" variant (which assumes that the groundwater utilization for the whole projection period (2013-2099) will remain the same as that of the last year of the historic period, 2012), that is, a notable ongoing decrease of the baseflow is obtained under all RCPs, again more strongly in the distant future period. These projected changes of the baseflow affect also the future baseflow index (BI), though, because of its definition (ratio of baseflow to total streamflow), somewhat differently as baseflow itself. Thus, BI rises spectacularly in the future, irrespective of the pumping variant employed. This is due to the dominant effect of the surface runoff (in the denominator of the BI-equation) which will strongly decrease in the future due to climate change (decreasing precipitation), and this is, of course, not affected by the pumping variant. In contrast to the high sensitivity of the future baseflow to the two pumping variants, this does not hold for the concentrated groundwater recharge produced by the losing (influent) streams. First of all, the total amount of concentrated recharge in the GRB is up to nearly two orders of magnitude smaller than the previously discussed concentrated discharge (baseflow). This is due to the fact that since most of the area of the modeled Ardabil aquifer is rather flat plain, where the groundwater table is near the land surface, positive head difference towards the river gage heights are produced, which leads to the sizable generation of baseflow. Furthermore, inasmuch as the future surface runoff was shown to be considerably abated under the different climate forcing scenarios, less streamflow will also be available within the river network for the concentrated groundwater recharge from possibly losing streams. It is interesting to note that, in contrast to the strong future projected augmentation of the baseflow under the "nonpumping" variant, the future concentrated groundwater recharge from the losing streams will be slightly reduced, because the future groundwater table is expected to fall less under this "nonpumping" variant and due to an additional small future increase of the diffuse surface-to-groundwater recharge (see Figure 7), which could cause a reversal of the head gradient between the river and the unconfined upper aquifer layer, so that, eventually, a losing stream (influent conditions) may switch to a gaining stream (effluent condition).
In addition to the temporal assessment of the climate change and pumping impacts on the losing and gaining streams, we investigated also their spatial variability. On this subject, the results illustrate that the maximum daily baseflow originating from all gaining streams will rise under all of the three climatic emission scenarios (RCPs) as well as for the two pumping variants, for instance, from 622.71 m 3 /day, for the historic period to 709.29 m 3 /day (in the upstream, hilly sections of the GBR) for time slice 2041-2070 (during which an increase of the precipitation is projected) under RCP 2.6 and "nonpumping". In spite of these sporadic and local increases, the future total daily baseflow will slightly decrease for the "pumping" variant, and this holds for all three RCPs. On the contrary, for the "nonpumping" variant the daily baseflow of 63415 m 3 /day obtained for the historic period will soar to values higher than 70000 m 3 /day, depending on the time slice, but irrespective of the RCP emission scenario.
Regarding the spatial variability of the concentrated groundwater recharge from the losing streams, similarly to its future temporal reactions to the climate change scenarios and the pumping variants mentioned earlier, only a low sensitivity of this parameter can be seen in the corresponding maps, unlike for the baseflow. Nevertheless, the future total daily volume of the concentrated groundwater recharge for the whole aquifer will considerably change under the various combinations of emission scenarios and pumping variants. Despite the fact that the number of losing streams, 1,659, is much higher than that of the gaining streams, 922, the concentrated groundwater recharge from the former is up to two orders of magnitude smaller than the groundwater discharge to the latter. The maximum of the "little" groundwater recharge from the losing streams of only 0.4 m 3 /day will occur over the western part of the Ardabil aquifer during the middle 10.1029/2019WR025388 Water Resources Research future (2041-2070), under RCP 2.6 and for the "pumping" variant. And dissimilar to the baseflow, the future total daily concentrated groundwater recharge volume over the aquifer is anticipated to rise under the "pumping" variant, irrespective of the RCP emission scenarios. For the "nonpumping" variant, on the contrary, the future total concentrated groundwater recharge from the losing streams will drop considerably, due to a lesser falling of the groundwater heads and the additional impact of a slightly increased future diffuse surface-to-groundwater recharge (see Figure 7), that forces some of the losing streams to converted to gaining streams.
Regarding the groundwater storage projected under the climatic emission scenarios and the two pumping variants, the results point out that the present-day unsustainable groundwater abstraction (continued into the future in the "pumping" variant) is the real agent driving for an ongoing drawdown, that is, groundwater storage decrease, in the future, even if the 2622 farming pumping wells that extract around 150 MCM water per year are completely shut off after the end of the historic period (2012). And not only that, the initial groundwater storage that existed in the initial years (initial condition of the simulation) of the historic period (e.g., 1980) cannot even be recovered, as the small absolute amount of the future diffuse groundwater recharge, albeit slightly increasing from historic 10 mm/yr to future 11 mm/yr, for RCP 2.6 and 4.5 (though not for RCP 8.5, because of less precipitation and higher evapotranspiration) (see Figure 7), is not able to replenish the depleted aquifer. In conclusion, the delicate balance between the total groundwater recharge and discharge for the Ardabil aquifer in the predevelopment has been irrevocably lost over the last two decades and will continue to do so even more in the future, accentuated by the adverse impacts of climate change (reduced precipitation in the near future, together with a small degree of diffuse surface-to-groundwater recharge in the GRB).
Eventually, the effects of the different climate emission scenarios and pumping variants on the total water yield of the GRB are determined. The results show that for both pumping variants, this quantity will drop for all future months, except for April in time slice 2041-2070, under the (benevolent) RCP 2.6, the time period for which an increase in the precipitation is projected. Surprisingly, "nonpumping" does not noticeably improve this adverse future water yield situation.This is, because the major fraction of the water yield is contributed by the surface runoff in the GRB and, for this reason, the same findings for the latter are obtained for the water yield as well, that is, a shift of its apex from historic-time March to April. Overall, less water will be available in the GRB in the future than in the historic period, wherefore, expectedly, this adverse effect will be the strongest for the most extreme RCP 8.5 climate emission scenario.
Although the results of the present study indicate that by combining a complex semidistributed surface water model (SWAT) with a 3D groundwater model (MODFLOW-NWT) significant enhancements in the understanding of the interaction of surface-subsurface hydrological processes can be achieved, there are several remaining, general and particular issues yet, that should be addressed in further studies. First of all, one concern is the general uncertainty attributed to the parameterization, parameter estimation, and to the observations required for setting up, driving and evaluating a hydrological model need (Wagener & Gupta, 2005;Y. Liu & Gupta, 2007). This is particularly true for integrated hydrological models such as SWAT-MODFLOW-NWT used here which requires a large number of input parameters related to the land surface, vadose zone, and groundwater domain hydrological processes, many of which could only be prescribed here on an ad-hoc, experience-based manner (Bailey et al., 2016;Taie Semiromi & Koch, 2019). Thus, in spite of our awareness of the ensuing uncertainties in the simulated surface-groundwater interactions, and so impacting the overall water balance in the study basin (GRB) region, such a general uncertainty analysis has not been performed in the present study. Suffice to say here that the hydraulic conductance of the stream-aquifer interface (river banks and bottom beds) appears to have a pivotal role in that regard.
Secondly, since future climate change scenarios are also faced with large uncertainties, the effects of latter on the surface-groundwater interactions should also be investigated. Thirdly, future land-use and land-cover change (LULCC), not considered in the present study, although certainly being subject to change by 2100, will have additional implications for the water resources budget of the GRB. The most recent SWAT-model version allows to implement updated LULCC during the ongoing simulation (Moriasi et al., 2019) and future investigations should take this into consideration. With these reservations, the water resources results obtained for the GBR in the present study should be seen more as trends rather than exact predictions. Notwithstanding these caveats, as the climate change impacts on the water components were projected under both current (historic) groundwater utilization ("pumping") and "nonpumping" variants, the here proposed fully coupled surface-subsurface hydrological model, that is, SWAT-MODFLOW-NWT, can make it possible to estimate the maximally sustainable groundwater utilization amount. Consequently, the adverse impacts of future climate change on the water resources can be alleviated. The influence of most likely occurring future increments of the groundwater demand was not considered in the current study, because Ardabil Water Authority has prohibited further groundwater utilization from the Ardabil aquifer for the purpose of recovering the ongoing sharp groundwater head drawdown. Ergo, model adjustments to that regard are recommended in future studies, since there is no doubt that the groundwater demand will rise in the future, not only due to a shortage of surface water (triggered by a decrease in precipitation and increase in temperature), but also due to population growth as well as expanded irrigated agriculture.
Although the findings of the present study may not be comparable to other hydrogeomorphic regions in the world because temperature and, particularly, precipitation vary spatially to a large extent over the globe, our proposed hybrid modeling scheme can successfully be applied to any other study region for future water resources management and planning in the wake of different anthropocentric and hydroclimate stressors, provided that GCMs and/or RCM-climate predictors under different GHG emission scenarios are available which are, however, new sources of uncertainty.