Integrated assessment of the impact of climate and land use changes on groundwater quantity and quality in the Mancha Oriental system ( Spain )

Climate and land use change (global change) impacts on groundwater systems cannot be studied in isolation. Land use and land cover (LULC) changes have a great impact on the water cycle and contaminant production and transport. Groundwater flow and storage are changing in response not only to climatic changes but also to human impacts on land uses and demands, which will alter the hydrologic cycle and subsequently impact the quantity and quality of regional water systems. Predicting groundwater recharge and discharge conditions under future climate and land use changes is essential for integrated water management and adaptation. In the Mancha Oriental system (MOS), one of the largest groundwater bodies in Spain, the transformation from dry to irrigated lands during the last decades has led to a significant drop of the groundwater table, with the consequent effect on stream–aquifer interaction in the connected Jucar River. Understanding the spatial and temporal distribution of water quantity and water quality is essential for a proper management of the system. On the one hand, streamflow depletion is compromising the dependent ecosystems and the supply to the downstream demands, provoking a complex management issue. On the other hand, the intense use of fertilizer in agriculture is leading to locally high groundwater nitrate concentrations. In this paper we analyze the potential impacts of climate and land use change in the system by using an integrated modeling framework that consists in sequentially coupling a watershed agriculturally based hydrological model (Soil and Water Assessment Tool, SWAT) with a groundwater flow model developed in MODFLOW, and with a nitrate mass-transport model in MT3DMS. SWAT model outputs (mainly groundwater recharge and pumping, considering new irrigation needs under changing evapotranspiration (ET) and precipitation) are used as MODFLOW inputs to simulate changes in groundwater flow and storage and impacts on stream–aquifer interaction. SWAT and MODFLOW outputs (nitrate loads from SWAT, groundwater velocity field from MODFLOW) are used as MT3DMS inputs for assessing the fate and transport of nitrate leached from the topsoil. Three climate change scenarios have been considered, corresponding to three different general circulation models (GCMs) for emission scenario A1B that covers the control period, and short-, mediumand long-term future periods. A multi-temporal analysis of LULC change was carried out, helped by the study of historical trends (from remote-sensing images) and key driving forces to explain LULC transitions. Markov chains and European scenarios and projections were used to quantify trends in the future. The cellular automata technique was applied for stochastic modeling future LULC maps. Simulated values of river discharge, crop yields, groundwater levels and nitrate concentrations fit well to the observed ones. The results show the response of groundwater quantity and quality (nitrate polPublished by Copernicus Publications on behalf of the European Geosciences Union. 1678 M. Pulido-Velazquez et al.: Global change impact on Mancha Oriental groundwater lution) to climate and land use changes, with decreasing groundwater recharge and an increase in nitrate concentrations. The sequential modeling chain has been proven to be a valuable assessment tool for supporting the development of sustainable management strategies.

Abstract.Climate and land use change (global change) impacts on groundwater systems cannot be studied in isolation.Land use and land cover (LULC) changes have a great impact on the water cycle and contaminant production and transport.Groundwater flow and storage are changing in response not only to climatic changes but also to human impacts on land uses and demands, which will alter the hydrologic cycle and subsequently impact the quantity and quality of regional water systems.Predicting groundwater recharge and discharge conditions under future climate and land use changes is essential for integrated water management and adaptation.In the Mancha Oriental system (MOS), one of the largest groundwater bodies in Spain, the transformation from dry to irrigated lands during the last decades has led to a significant drop of the groundwater table, with the consequent effect on stream-aquifer interaction in the connected Jucar River.Understanding the spatial and temporal distribution of water quantity and water quality is essential for a proper management of the system.On the one hand, streamflow depletion is compromising the dependent ecosystems and the supply to the downstream demands, provoking a complex management issue.On the other hand, the intense use of fertilizer in agriculture is leading to locally high groundwater nitrate concentrations.In this paper we analyze the potential impacts of climate and land use change in the system by using an integrated modeling framework that consists in sequentially coupling a watershed agriculturally based hydrological model (Soil and Water Assessment Tool, SWAT) with a groundwater flow model developed in MODFLOW, and with a nitrate mass-transport model in MT3DMS.SWAT model outputs (mainly groundwater recharge and pumping, considering new irrigation needs under changing evapotranspiration (ET) and precipitation) are used as MODFLOW inputs to simulate changes in groundwater flow and storage and impacts on stream-aquifer interaction.SWAT and MODFLOW outputs (nitrate loads from SWAT, groundwater velocity field from MODFLOW) are used as MT3DMS inputs for assessing the fate and transport of nitrate leached from the topsoil.Three climate change scenarios have been considered, corresponding to three different general circulation models (GCMs) for emission scenario A1B that covers the control period, and short-, medium-and long-term future periods.A multi-temporal analysis of LULC change was carried out, helped by the study of historical trends (from remote-sensing images) and key driving forces to explain LULC transitions.Markov chains and European scenarios and projections were used to quantify trends in the future.The cellular automata technique was applied for stochastic modeling future LULC maps.Simulated values of river discharge, crop yields, groundwater levels and nitrate concentrations fit well to the observed ones.The results show the response of groundwater quantity and quality (nitrate pol-

Introduction
Future climate and land use changes might substantially modify the hydrological processes and consequently, having an impact on groundwater quantity and quality.While climate change will alter hydrological conditions due to changes in the major climate variables (air temperature, precipitation and evapotranspiration), groundwater resources will be impacted by climate change not only through their interaction with surface water bodies (e.g., lakes and rivers), but also, indirectly, through the recharge process (Jyrkama and Sykes, 2007).Therefore, quantifying climate change impacts on groundwater requires assessment of future evolution of the hydrological variables that interact with the groundwater system (river stages, groundwater recharge, pumping, pollutant leaching, etc.).However, a proper assessment of those variables is difficult, as they depend on multiple physical factors subject to high temporal and spatial variability.
Recent climate and land use change studies in the Mediterranean region pointed at temperature and evapotranspiration increase, especially during summer (Chaouche et al., 2010;Candela et al., 2012;Molina-Navarro et al., 2014;Ertürk et al., 2014).Although the existence of a decreasing trend in the yearly precipitation has not been unanimously reported in all the studies, they agree in pointing at a change in the in-year precipitation pattern, with higher precipitation during fall and winter and lower rainfall during summer (Chaouche, 2010;Molina-Navarro et al., 2014).Those modified patterns, combined with land use changes, will decrease the amount of water recharge to the aquifers (Candela et al., 2009;Ertürk et al., 2014).This will provoke lower water tables and decreased groundwater discharge, which may consequently reduce stream baseflow and impact water supply and groundwater dependent ecosystems (Kløve et al., 2014).These effects will be higher in heavily committed groundwater bodies, and will be intensified if water abstraction is increased to meet a growing demand for water.Predicting recharge and discharge conditions under future climatic and land use changes is essential for integrated water management and adaptation.
Numerical simulation models representing the spatial heterogeneity and temporal variability, if properly designed and calibrated, provide the most adequate way to estimate the impacts of climate and land use changes on groundwater systems.Numerical hydrological models have been used to estimate global change impacts on the surface processes of a watershed (e.g., Caballero et al., 2007;Mango et al., 2011;Shrestha et al., 2013;Ma et al., 2014).Those models are also useful tools to deal with the uncertainties on the impacts of climate change associated with the hydrological processes (Kingston and Taylor, 2010;Xu et al., 2011).In order to assess the impacts of future conditions (climate, land use, water demands, adaptation, etc.) on groundwater systems, some form of coupling between hydrological and hydrogeological processes must be used (Holman et al., 2012;Pulido-Velazquez et al., 2015).However, integrated models' development continues to be a challenge, as it requires assumptions that hinder detailed assessments about certain variables.The sequential coupling of numerical models can be used as an alternative to adequately assess climate and global change impacts.Those approaches generally employ hydrological models capable of representing the land phase of the hydrological cycle (infiltration, recharge, river flow, surface nitrate flow, nitrate leaching, etc.), while groundwater quantity and quality impacts are assessed through groundwater flow and transport models (e.g., Candela et al., 2009;Narula and Gosain, 2013).Their linkage is made through an output-input scheme, in which some of the hydrological models' outputs (groundwater recharge, nitrate leaching, etc.) are used as inputs in the groundwater models.Therefore, sequential coupling is able to keep the necessary detail in the key processes in the system, while providing a comprehensive description of the way those processes are interconnected.
The main goal of this study is the integrated assessment of climate and land use change impacts in the Mancha Oriental groundwater system (Spain), considering groundwater quantity and quality (nitrate pollution).The approach consists of the sequential coupling of the hydrological model Soil and Water Assessment Tool (SWAT), the groundwater flow model MODFLOW, and the groundwater transport model MT3DMS.Section 2 presents the case study, describes the methodology and provides information about the calibration and validation processes.Section 3 presents the results for the combinations of climate and land use change scenarios, key findings and novelty features.Section 4 incorporates the main conclusions for the case study and for the methodology 2 Materials and methods

Case study: the Mancha Oriental system in Spain
The Mancha Oriental system (MOS) is located in an area of semiarid climate in the southwestern part of the Jucar River basin, mainly within the Albacete province, Spain (Fig. 1).The study area covers about 8400 km 2 consisting of plains surrounded by mountain ranges that delimitate the borders between the Jucar, the Tajo (Tagus), the Guadiana and the Segura river basins.The northern area is dominated by three main rivers: the perennial Jucar River and its seasonal tributaries, Valdemembra and Ledaña.The southern portion is a former endorheic plain which was artificially connected to the Jucar River in the nineteenth century by the Maria Cristina Channel.The main land use of the area is agriculture, being especially characteristic the circular-shape groundwater-irrigated crops, devoted mainly to corn, wheat and barley.
The Mancha Oriental aquifer is made of three main layers over a Triassic (Keuper) impervious base.The bottom layer is a confined Jurassic limestone aquifer with a thickness of 250 m, which emerges, getting unconfined, in certain areas in the south and west of the study zone.The mid-layer is formed with Cretaceous limestone mainly confined but emerging in the north and east.The upper layer, which holds the majority of the pumping wells, is a Miocene limestone unconfined aquifer with some clay intercalations, up to 150 m thick.A detailed geological description can be found in Sanz (2005) and Sanz et al. (2009Sanz et al. ( , 2011)).
In the last 25 years, an important transformation from dry to irrigated lands has taken place in the MOS, with the development of an intensive irrigated agriculture that represents one of the main factors in the current economic development of the region.More than 80 000 ha. of land equipped with modern technologies are currently irrigated, with most irrigation coming from groundwater.The main crops are wheat, corn, barley and alfalfa, with a significant share of the crop production still dependent on subsidies from the EU Common Agricultural Policy (CAP), and with some growing areas of vegetables and vineyard.The aquifer has been subject to an intensive groundwater overexploitation since the 1980s, which has resulted in a continued drop of groundwater levels, especially in the southern area where irrigated crops concentrate.Groundwater levels have declined up to 60 m with respect to the levels in the 1980s, with up to 3000 Mm 3 of groundwater storage depletion from 1980 to 2005 (Sanz et al., 2011).In 2009 the groundwater body was classified as in bad quantitative status by the Jucar River Basin Agency (Confederacion Hidrografica del Jucar, CHJ) due to the unbalance between the available groundwater resource (as defined by the EU Water Framework Directive) and abstractions.The water table decline has caused the extinction of wetlands and lagoons, such as the Acequion or Salobral lagoons.The stream-aquifer interaction with the Jucar River has been also substantially affected: formerly, the aquifer discharged into the Jucar River flow, while today the river recharges the MOS (Sanz et al., 2011).This has led to a significant depletion of streamflow in the Jucar River with important environmental consequences (such as the drying of a significant reach of the Jucar River in the summers of 1994 and 1995), provoking conflicts with downstream uses.The Mancha Oriental groundwater body is in bad quantitative status, and it also fails in reaching the good chemical status because of the increasing groundwater nitrate pollution due to the intensive use of fertilizers in agriculture (CHJ, 2009a).Nitrate concentrations in the aquifer system have locally reached values up to 125 mg L −1 (Moratalla et al., 2009), far away of the standard of 50 mg L −1 .Different strategies for controlling groundwater nitrate pollution in the area have been proposed, including fertilizer quotas and fertilizer taxes (Peña-Haro et al., 2010, 2014).
In order to deal with these issues, the CHJ is proposing, in the framework of the new Jucar River Basin Management Plan (PHJ), actions for a more sustainable management of the Mancha Oriental aquifer, including demand reduction (mainly by irrigation efficiency improvement in some cases), some substitution of water sources from groundwater to Jucar water diversions, and an intense research on the aquifer and water uses behavior (CHJ, 2009b).The control of groundwater abstractions and water use by remote-sensing and personal inspections (Castaño et al., 2010), combined with collective actions through groundwater user associations (Lopez-Gunn, 2003), are helping to stabilize groundwater abstractions.In this context, climate change is likely to exacerbate the groundwater management issues, as decreasing groundwater recharge is likely to increase the pressure on the quantity and quality status of the aquifer.

Climate change scenarios
The climate change scenarios rely on the Intergovernmental Panel on Climate Change (IPCC) Special Report on Emissions Scenarios (SRES) A1B emission scenario.Climate data were obtained from three different general circulation model (GCM) drivers: CNRM (National Centre of Meteorological Research), ECHAM5-r3 (European Centre for Medium-Term Weather Forecast) and HADCM3-Q0 (Hadley Centre).These scenarios have been downscaled using the SMHIRCA 3.0 RCM (Regional Climate Model) of the Swedish Meteorological and Hydrological Institute (SMHI).This study is one of the 16 case studies in the EU GENESIS Project, which deals with climate and land use impacts on groundwater and dependent ecosystems.For all case studies, the same dynamic downscaling method was applied by SMHI, which provided the meteorological forcing time series for the climate change scenarios (Kjellstrom et al., 2011;Nikulin et al., 2011).Daily time series of the relevant meteorological variables were provided for the 1961-2100 period, corresponding to a control period , and the short-term (2010-2040), medium-term (2040-2070), and long-term (2070-2100) scenarios.
A comparison was made for the control period  between the climate model simulations and the historical time series, in order to check if they fit the observed patterns in the main statistics of temperature and precipitation.Historical data sets were obtained from the Spanish Meteorological Agency (AEMET).Figure 2 shows the comparison of monthly means and standard deviations for maximum temperature (at Albacete-Los Llanos station) and total precipitation.The CNRM scenario is the one that better resembles the mean maximum historical temperature, with larger differences in the other two.Regarding monthly standard deviations of maximum temperature, the HADCM3 scenario shows the largest discrepancies with respect to the historical data.Figure 2 shows some differences on the mean precipitation monthly pattern, especially in April and November.Some deviations are also found in the comparison of the standard deviations.There is not a single model that clearly outperforms all in terms of agreement with the observation.In any case, we did not use this comparison for selecting climate projections, since the performance of the model in the past control period does not guarantee improvements in prediction accuracy under a non-stationary climate (Reifen and Toumi, 2009;Teutschbein and Seibert, 2012).Instead, we decide to use a set of models that provides a range of plausible future projections.This will allow for assessing the sensitivity of the model results to the uncertainty on the climate projections.
The temperature and precipitation series for the three climate scenarios show a steady increase in temperature, with a decreasing trend in precipitation, more accused in the longterm (Fig. 3).The other meteorological variables analyzed (relative humidity, solar radiation and wind speed) did not show any clear trend.Finally, CO 2 concentrations were derived from the Integrated Science Assessment Model (ISAM) and BERN carbon cycle models, both used in the IPCC 2007 Fourth Assessment Report's climate projections, for the SRES A1B emission scenario (Nakicenovich and Swart, 2000).A concentration value of 421.67 ppmv was considered for the short-term period, while concentrations of 531.67 and 665.50 ppmv were considered for the medium-and longterm, respectively.

Land use change projections
Four land use change (LUC) scenarios have been considered in this case study: -NLUC: no land use change or baseline scenario.This scenario has been included to be compared against the other ones, allowing the identification of synergic effects of climate and land use changes.It has been defined by all the three time periods (from 2010 to 2100), and its land use pattern corresponds to the current one.
-ST: short-term scenario based on a multi-temporal analysis of historical LUC changes and their key drivers, future EU scenarios and a combination of LUC allocation techniques.This scenario has been defined for the 2010-2040 period.
-MLTII: medium-and long-term scenario with increased irrigation.It considers the LUC trend observed in the last 20 years, characterized by a change from nonirrigated to irrigated crops of about 10 % of the agricultural area with respect to the ST scenario.We can assume that this is partially supported by subsidies coming from the EU CAP.This scenario has been defined for 2040-2070 and 2070-2100.
-MLTDI: medium-and long-term scenario with decreased irrigation.It considers the policies recently set up by the Farmers' User Association of Mancha Oriental (JCRMO), and potential energy and water price increases, the latter coming from the application of the pricing policies required by the EU Water Framework Directive (WFD).These driving forces would cause a decrease in the irrigated area of about 20 % with respect to the ST scenario, which would return to be operated without irrigation.This scenario has been defined for 2040-2070 and 2070-2100.For the short-term LUC scenario, a multi-temporal analysis of land use change (Oñate-Valdivieso and Bosque-Sendra, 2010) was carried out, based on the study of historical trends in crop patters (derived from remote-sensing images annually processed for the area; Calera et al., 1999Calera et al., , 2005Calera et al., , 2012) ) and key driving forces for explaining LUC transitions.A spatially explicit model was use for that purpose, the Land Change Modeler Module (Eastman, 2006).This model provides a set of tools to perform historical analysis, trends and future scenario projections based on geographic information system (GIS) techniques.Data sets include series of remote-sensing images, driving forces and scenarios from regional projects.LUC images from CORINE Land Cover Project were used as baseline (1990 and2006;Fig. 4).
CORINE Land Cover (CLC) images have 100 m resolution and provide a LUC classification scheme widely used in all of Europe (Feranec et al., 2010).Based on these images, historical analysis was performed to identify, evaluate and select the most significant Land use and land cover (LULC) transitions.To develop LUC future scenarios, all transitions need to be modeled according to the most likely driving forces.Driving forces were chosen from the available literature and statistical data with spatial representativeness (INE, 2011) from a wide range of biophysical, social, economic and political factors.A set of 20 driving forces were finally selected and spatially represented in raster format (Table 1), using a Cramer's V test to select relevance of each driving force in every transition with a threshold value of 0.15 (Eastman, 2006).
Markov chains and European scenarios and projections (Eururalis 2.0 and Image 2.2) have been used to define trends in future land use in Europe.The EU project Eururalis 2.0 (Klijn et al., 2005) was chosen because its focus on rural areas.Eururalis developed Europe's future land use scenarios for 2010, 2020 and 2030 (Westhoek et al., 2006;  Rienks, 2007), all available online to assist decision makers on the most likely scenarios of the Common Agricultural Policy (CAP).The predicted LUC scenarios were obtained from a combination of suitability maps and transition probability matrices extracted from all data sets described.A multi-criteria evaluation (MCE) method using Artificial Neural Networks (ANN) was applied to obtain suitability maps (Oñate-Valdivieso and Bosque Sendra, 2010) where change could happen for each transition selected.The multilayer perceptron (MLP) with three layers, backpropagation of error algorithm and the sigmoid transformation function were used to train every ANN.The transition probability matrices (TPM), based on the Markov's chain theory, express the probability of one LUC to change to another in a determined period of time.TPM were obtained from historical data and from the Eururalis land use scenarios.Applying the Markov process theory to the historical data (2000)(2001)(2002)(2003)(2004)(2005)(2006), transition matrices for the 2000-2010 period were derived.For the 2000-2020 period, transition matrices were obtained from the tendencies showed in the Eururalis scenario.Finally, the cellular automata algorithm fed with the suitability maps of every transition and the transition probability matrices produced the final suitability maps and the land use scenarios.A detailed explanation of the method can be found in Henriquez-Dole (2012).

Modeling framework
In order to assess the impacts of the predicted future conditions (climate, land use, water demands, adaptation, etc.) on groundwater systems, a link has to be established between those forcings and the groundwater system responses.This requires the practical integration of operational models that not only represent all of the relevant processes in the hydrologic system in a physically meaningful way but also are simple enough to allow large-scale basin-wide applications (Sophocleous and Perkins, 2000).Different processes of the hydrologic cycle need to be modeled.The characterization of the land phase of the hydrological cycle is essential for assessing the impacts of climate and land use changes on the temporal and spatial distribution of groundwater recharge and contaminant loadings.Moreover, in basins in which irrigated agriculture is the dominant land use, as in our case study, we also need to involve calculations of plant growth and crop yields, crop evapotranspiration (ET), irrigation applications and groundwater pumping changes.Another additional requirement in this case was to include simulation of nitrate leaching from the crop fertilization, in order to assess impacts of future scenarios on groundwater nitrate pollution.The tool we selected for this purpose was the SWAT, one of the most widely used watershed models worldwide, applied extensively to a broad range of scales and environmental conditions (Gassman et al., 2007(Gassman et al., , 2014)).
While the quasi-distributed SWAT model is capable of properly simulating the spatiotemporal distribution of groundwater recharge rates (at the spatial resolution given by their hydrologic response units, HRUs), its groundwater module is lumped.Therefore, distributed groundwater parameters (such as hydraulic conductivities and storage coefficients) cannot be represented, and the approach is very limited for expressing the spatial distribution of groundwater levels and groundwater flow dynamics.Different SWAT-MODFLOW integration and coupling approaches have been proposed in the literature to deal with this issue (Sophocleous and Perkins, 2000;Kim et al., 2008).
The proposed modeling framework sequentially couples the SWAT watershed model with the fully distributed groundwater model MODFLOW (McDonald and Harbaugh, 1988), and finally the multi-species transport model MT3DMS (Zheng and Wang, 1999) for simulating the fate of the nitrate leached into the aquifer system.In this approach, SWAT model outputs are used as MODFLOW inputs, and SWAT and MODFLOW outputs are used as MT3DMS inputs (Fig. 5).This coupled framework was used to assess the climate and land use change impacts on hydrology, groundwater flow and nitrate leaching and transport in the MOS.

Watershed modeling
SWAT is a conceptual basin-scale continuous-time quasidistributed watershed simulation model for predicting the impacts of management on hydrology, sediments, agricultural production and chemical yields (Arnold et al., 1998;Neitsch et al., 2005;Gassman et al., 2007).Calculations are done on a daily basis, offering results at both daily and monthly timescales.
The SWAT interface tool in ArcGIS (ArcSWAT), is used to develop the model input data sets.The SWAT model requires a wide range of data depending on the modeled processes.It firstly divides the study area into different sub-basins regarding the river network.The sub-basins are further discretized into HRUs consisting of homogenous soil, slope and land use combinations (equal combinations associated with different sub-basins remain as separated HRUs).At the HRU scale, SWAT simulates the processes specified by the user, being able to perform calculations related to hydrological processes, sediment transport, nutrient cycle, soil temperature, crop growth and pesticides management (Arnold et al., 1998).

DEM
Soil data

SWAT model
Hydro-geological data Geological data

Groundwater abstractions
Groundwater recharge

Nitrate leaching Groundwater flows
Groundwater velocities

Groundwater nitrate concentrations
Figure 5. Modeling framework adopted.
Water balance is the driving force behind all the processes in SWAT because it impacts plant growth and the movement of sediments, nutrients, pesticides and pathogens (Arnold et al., 2012).SWAT calculates the hydrology at each HRU by means of the water balance equation, which includes daily precipitation, evapotranspiration, percolation, runoff and return flow components.Every process included in the model can be solved using different methodologies.In our case, the surface runoff derived from daily rainfall is estimated with a modification of Soil Conservation Service (SCS) curve number method (included in the model).The percolation through each soil layer is predicted using storage routing techniques combined with crack-flow model (Arnold and Williams, 1995).The evapotranspiration was estimated using Hargreaves formula.Finally, the flow routing in the river channels is computed using the variable storage coefficient method (Williams, 1969).
SWAT uses a single plant growth model to simulate all types of land cover and differentiates between annual and perennial plants.The plant growth model is used to assess removal of water and nutrients from the root zone, transpiration, and biomass/yield production (Arnold et al., 2012).Planting, harvesting, tillage passes, nutrient and pesticide applications can be simulated for each cropping system with specific dates or with a heat unit scheduling approach.The irrigation applications can be simulated for specific dates or with an autoirrigation routine, which triggers irrigation events according to a water stress threshold.
The nitrogen (N) processes and soil pools simulated by SWAT are described in Neitsch et al. (2005).SWAT monitors five different pools of N in the soil.Two of them are inorganic forms of N: NO − 3 and NH + 4 .The other three are organic forms of N: fresh organic N associated with crop residue and microbial biomass, and the stable organic pool associated with the soil humus.SWAT is capable of simulate N fixation by legumes, fertilizer inputs and nitrogen in the rainfall as well.

M. Pulido-Velazquez et al.: Global change impact on Mancha Oriental groundwater
The SWAT model for the MOS has been fed with the following inputs, divided into four categories: the digital elevation model (DEM), climate data, soil data and land use data.The DEM used in SWAT has a 670 m × 670 m cell size, being used to delimitate 35 sub-basins (using the ArcSWAT watershed delineator tool) and to derive the slope map.
Daily scale historical climate data records have been obtained from the Climatic Research Unit (CRU) website (available at http://badc.nerc.ac.uk/data/cru/),AEMET raster database (available at http://escenarios.aemet.es/),and data records obtained in several stations within the case study area from the AEMET and the Integrated Service of Irrigation Advising (SIAR).Air relative humidity, wind speed and solar radiation were only available for the full period in the Albacete-Los Llanos weather station.In order to complete the data gaps, a spatial correlation analysis has been carried out, using the available data in all the weather stations.Correlation coefficients were obtained with respect to the Albacete-Los Llanos for estimating the unrecorded values as a function of the observed ones through regression analysis.Humidity correlation coefficients ranged between 0.82 and 0.97, wind speed between 0.59 and 0.96 and solar radiation between 0.90 and 0.99.
Soil type data were obtained from the FAO's Digital Soil Map of the World (there is no Spanish database for all Spanish territory).The original resolution of the FAO's world digital soil map is 1 : 5 million, which presents soil parameters in grid cells of 5 latitude/longitude resolution.In the case study area Glegyc Cambisol (Bg), Chromic Luvisol (Lc) and Calcic Cambisol (Bk) soil types have been found with different textures, representing a total of seven soil categories (1 for Bg, 1 for Lc and 5 for Bk).Land use data have been obtained from the CLC project for year 2000.CLC data were obtained at 1 : 100 000 scale, with a 100 m × 100 m resolution mesh and a minimum land use area of 25 ha.The CLCobtained data have been compared with land images (Calera et al., 2005(Calera et al., , 2012;;Henriquez-Dole, 2012), in order to check its availability.The agriculture land was divided into a crop mosaic consisting of corn (40.4 %), wheat (23.6 %), barley (18.0 %) and alfalfa (9.0 %), onion (4.5 %), and sugar beet (4.5 %).
The HRUs were obtained as combination of 12 land use categories (Fig. 4), 7 soil categories and 1 slope category.In a preliminary SWAT HRU definition, an excessive amount of HRUs were found.In order to reduce it, a filter was applied by not considering, within one specific sub-basin, land use and soil types whose surface area is below 10 % of the total sub-basin area excluding irrigated land.After these operations, a final number of 445 HRUs were obtained.Then, the crop management practices (fertilizer application, seed time, irrigation and harvest timetable; the water source, hydric stress threshold, etc.) were introduced in the SWAT model using the ArcSWAT interface.The crop management parameters were based on the standard practices of the farmers in the watershed.

Groundwater flow model
MODFLOW is a fully distributed model that solves the three-dimensional groundwater flow equation using finitedifference (FD) approximations.The model calculates the hydraulic head at each cell of the FD grid (for which the aquifer properties are assumed to be uniform), and from there, the flow between cells, stream-aquifer or lakegroundwater interaction, flows through drains.For that purpose, it requires geological and hydrogeological inputs such as top and bottom layer elevations, hydraulic parameters at the grid (hydraulic conductivity, storage coefficient), boundary and initial conditions and external stresses.
The MOS groundwater model consists of seven hydrogeological units (HU), three of them are considered as aquifers (HU2, HU3 and HU7) and the other as aquitards (Sanz et al., 2011(Sanz et al., , 2009)).The hydrogeological unit 7 (HU7) is present throughout the MOS and is composed of limestone and fractured dolostone.The HU3 is only present in the northeast part of the study area and is composed of fractured limestone and dolostone.The upper aquifer, the HU2, which is located in the central part is composed of an alternate sequence of marl-lime and marl.Six hydrogeological domains can be identified in the MOS: northern (ND), central (CD), El Salobral-Los Llanos (SLD), Moro-Nevazos (MND), Pozo Cañada (PCD) and Montearagón-Carcelén (MCD) domains.According to Sanz et al. (2011), there is hydraulic connection among the ND, CD and SLD domains, but not among MND, PCD and MCD or among these three and ND, CD and SLD.The Jucar River is the most important surface body and it is hydraulic connected to the aquifer, mainly to the HU2.
The groundwater flow was simulated using MODFLOW (McDonald and Harbaugh, 1988).The model was discretized into 114 columns, 129 rows and 6 layers, with a cell size of 1 km 2 .The main aquifers are represented by layers 2, 4 and 6 (HU7/HU6/HU5, HU3 and HU2), while the other layers are semipermeable units (HU4, lower HU1 and upper HU1).The total number of pumping wells considered was 1776, most of them located in layers 2 and 6.To the west of the studied area, the western Mancha System is assumed to contribute with groundwater discharge into the MOS.This was simulated by a constant-head boundary condition.The initial heads, hydraulic conductivity and storage coefficients were obtained from Sanz (2005), although the latter values were further modified during the calibration process.The hydraulic conductivity varies between 0.05 and 500 m day −1 and the storage coefficient, between 1 × 10 −4 and 1 × 10 −5 .Recharge and pumping values were obtained from the SWAT model outputs, assuming that farmers pump the optimal amount of water required by their crops.

Groundwater nitrate transport model
MT3DMS (Zheng and Wang, 1999) is a 3-D groundwater solute transport model that solves the groundwater transport equation using a finite-difference approximation, discretizing the spatial domain into cells in which equal characteristics and solute concentrations are assumed.Nitrate transport parameters for the Mancha Oriental aquifer were estimated using the values reported in the literature.The model simulates advection and diffusion in nitrate transport in the saturated zone.
For the MOS model, initial concentration values were interpolated from data reported in Moratalla et al. (2009).Nitrate leaching loads entering the aquifer were obtained from SWAT outputs.Only agricultural sources were considered.Since only few scattered measured values of nitrate concentrations were available for calibration, the calibration mainly consisted in matching the maximum nitrate concentrations simulated with the ones reported in Moratalla et al. (2009).

Models' calibration
In the proposed modeling framework, given the sequential use of models, the SWAT, MODFLOW and MT3DMS models were calibrated independently.

SWAT model
The SWAT calibration procedure followed four steps, consisting of river flow, groundwater recharge, crop yield and nitrate leaching calibrations.A preliminary sensitivity analysis was run using the SWAT-CUP software (Abbaspour, 2012), which integrates various calibration/uncertainty analysis procedures for SWAT.From this analysis we obtained that the most sensitive parameters of the model were the SCS curve number (CN2), the groundwater discharge coefficient (α) into surface water, the travel time coefficient from soil to shallow aquifer, the coefficient of water finally lost as evaporation and the two primary evapotranspiration parameters.The fact that the CN2 parameter was the most sensitive highlights the importance of considering land use changes in the analysis (as they affect the CN2).
Two river gauging stations were selected to calibrate the SWAT model's surface hydrology.The first one, Los Frailes (08036 station), is located at the very center of the case study zone, close to the confluence of the Jucar and Valdemembra rivers.The second one, Alcala del Jucar station (08144), is placed downstream of the confluence between the Jucar and Ledaña, and therefore downstream of the reach of stream-aquifer connection within the Mancha Oriental aquifer (Fig. 1).Daily flow data from 1994 to 2004 were used to calibrate the SWAT model, with data from 1991 to 1993 as the warm-up period.Figure 6 depicts the comparison between the simulated and observed monthly values.Both plots show a proper performance of the SWAT model, which ade-quately reproduces the Jucar River discharges.For the 08036 station, we obtained a monthly correlation coefficient of 0.92 and a Nash-Sutcliffe efficiency (NSE) coefficient equal to 0.84, while for station 08144 those values were 0.83 and 0.68, respectively.The lower performance shown by the calibration coefficients in Alcala station seems to be related to the higher drainage area, increasing the complexity for the calibration.The validation of SWAT for streamflow was done for the 2005-2006 period with acceptable results: daily correlation coefficients equal to 0.80 in Los Frailes gauging station, and to 0.52 in Alcala del Jucar, while daily Nash-Sutcliffe indexes were equal to 0.68 in Los Frailes and 0.33 in Alcala del Jucar.Monthly coefficients were not obtained since there are only 24 monthly data records available.
Groundwater recharge values obtained from SWAT have been compared with the ones reported in previous studies in the same case study.Sanz et al. (2011) estimated a mean recharge of 320 Mm 3 year −1 , while the Jucar River Basin Management Agency (CHJ, 2013) assessed a recharge between 238 and 334 Mm 3 year −1 for the 1994-2004 period.The SWAT model yields a mean annual groundwater recharge of 310 Mm 3 for the same period.
The average net value obtained with SWAT for streamaquifer interaction was also compared to those reported by the Geological Survey of Spain (IGME) for the area.The values previously reported by the IGME expressed a range from 40 to 60 Mm 3 year −1 (IGME-DGA, 2010), where SWAT reported an average value of 45 Mm 3 year −1 .
Unlike other previous studies using SWAT (e.g., Narula and Gossain, 2013), SWAT calibration has been also based on the crop simulation component.The simulated irrigation volumes and crop yields have been compared with the ones reported from crop surveys and experimental data in the zone.These data have been obtained from the deliverables of the Agronomic Technical Institute of the Province (Instituto Técnico Agronómico Provincial, ITAP) (López Urrea et al., 2003).As shown in Table 2, the SWAT model performance is consistent with the expected values.
Regarding nitrate leaching, there is no possible comparison with historical data, as there are no historical records.Given the usual absent of data, the calibration of nitrate loads is often based on observed nitrate concentrations in the existing piezometers and the crop growth component.Nitrate leaching calculations use SWAT leaching equations, information on fertilizer types and application, soil characteristic (including initial nitrate concentration) and other factors such as the percolation factor affecting nitrate transport.
In order to check the SWAT model performance, the simulated nitrate leaching values have been compared with those provided in previous studies about nitrate leaching for different crops (Martin-Benlloch, 2012) using the Environmental Policy Integrated Climate Model (EPIC) modeling tool.As Table 2 shows, SWAT outcomes are in agreement with those reported in the detailed study using EPIC.

MODFLOW model
MODFLOW calibration has been carried for the same period as the SWAT model (1994-2004 period), using 24 piezometers.Figure 7 depicts the calibration results obtained in 4 piezometers located in different zones of the Mancha Oriental aquifer.The seasonal fluctuations are mainly due to pumping during the irrigation season.The R 2 coefficients for the goodness of fit for heads are 0.69, 0.62, 0.27 and 0.72.The model performance closely resembles the historical records at the observation wells and, therefore, the MOD-FLOW model has been adequately calibrated.

MT3DMS model
The parameter adjustment for the MT3DMS model was hindered by the lack of data, especially by the lack of a continuous time series with at least one record per month.The initial concentration values were interpolated from data reported in the literature (Moratalla et al., 2009).The nitrate load entering the aquifer was obtained from the nitrate leaching variable in SWAT.Since only few scattered observed values of nitrate concentrations are available for the calibration, the process mainly consisted in matching the maximum nitrate concentrations simulated and their spatial distribution pattern with the ones reported in the literature (Moratalla et al., 2009).Figure 8 shows the results of this match, consisting of a graph in which the maximum nitrate concentrations obtained by MT3DMS are showed, comparing them to the literature-reported maximum value in the observation points marked (plotted besides each location).As can be seen, MT3DMS results fit the observed values for the majority of the control points.

Scenario results
Each scenario run refers to a single combination of climate change, period and land use scenarios, up to a total number of 12 (Table 3).

Impacts on groundwater recharge
With changing precipitation and temperature patterns and values ( climate change scenarios agree in a reduction of the mean recharge over the twenty-first century, especially in the longterm (Fig. 9), and a reduction with respect to the historical recharge in the last decades, estimated as 310 Mm 3 year −1 .The values of the short-term mean recharge are similar for the CNRM and HADCM3 scenarios, despite the fact that the precipitation is higher in CNRM than in HADCM3 (Fig. 3).
In this case, the difference in precipitation is not fully transferred to recharge.However, in the ECHAM5 scenario a lower precipitation results in a lower groundwater recharge.
A remarkable recharge reduction can be noticed in the longterm for the ECHAM5 and CNRM scenarios, associated with the steep decrease in precipitations (Fig. 3), whereas the HADCM3 scenarios show no decrease in either average recharge or precipitation.Furthermore, it is observed that the simulations with land use MLTII show the highest groundwater recharge, independently of the climate change scenario.This pattern is associated with the increase of irrigation area in this scenario (the irrigated area reaches a maximum of 1000 km 2 in MLTII).Assuming optimal pumping for plant growing as SWAT does, greater irrigated areas lead to higher percolation rates.Irrigation prevents the soil water deficiency in semiarid regions during the dry periods, enabling deep percolation in the early rainfall and, when applied in excess, irrigation water may also provide an additional percolation flux (Seiller and Gatt, 2007), although in this case the irrigation efficiencies are high (the most common irrigation system is center pivots).
Climate change seems to be the main driver of groundwater recharge change, given that the differences between climate scenarios with the same land uses (up to around 30 % in the long term) are larger than the differences when considering different land use scenarios (up to around 5 %).

Impacts on groundwater quantity
Figure 10 shows the evolution of the groundwater levels associated with the aquifer layers 2 and 6, which hold the   The results show cycles of decreasing and increasing groundwater levels, but without a generalized descending trend (Fig. 10).In layer 6, the MLTDI scenarios present higher levels than those corresponding to MLTII, with the land use change as the major driving force (scenarios with the same land use scenarios but with different climate change scenarios are very similar).On the other hand, scenarios with different land use scenarios but the same climate change scenario show higher differences.For layer 2, the groundwater levels showed for the MLTDI scenarios are higher than for the MLTII, but the oscillations of levels between scenarios are more significant and, during some periods, the levels in the MLTDI scenarios are located below the MLTII ones.To sum up, the effect of land use change on groundwater  levels seems to be higher than the effect of climate change since the LUC modify the amount of pumping, while causing just slight changes in the recharge.Table 5 shows the mean annual pumping for the different scenarios and periods, which can be compared with the historical pumping of 406 Mm 3 year −1 (Sanz et al., 2011).
The lower amount of pumping in the MLTDI scenarios leads to higher groundwater levels.However, no general decline trends can be noticed in the whole set of scenarios.The cycles noticed in the time series are mainly associated with climate variability.

Impacts on groundwater quality
Groundwater quality impacts of climate and land use changes have been analyzed in terms of groundwater nitrate concentrations using the MT3DMS model with the SWAT model outputs (recharge, pumping, nitrate leaching).With regard to nitrate leaching, all the future scenarios show an increase with respect to the average value reported by SWAT for the calibration period, which was equal to 56 kg NO 3 ha −1 .This situation is consistent with recent studies in which an increase of nitrate leaching in groundwater was predicted for the twenty-first century (Stuart et al., 2011).Regarding Fig. 11, the scenarios show higher differences among scenarios when advancing through the twenty-first century.While similar leaching values were found for the short term (between 75 and 85 kg NO 3 year −1 ), the results for the medium term and long term show a broader range, from 65 to 105 kg NO 3 year −1 .The scenarios associated with the MLTII land use scenario (the largest irrigated area) are generally the ones with higher nitrate leaching, while the MLTDI scenarios offer the lower bound of leaching values (Fig. 11).
Nitrate concentrations in the Mancha Oriental aquifer increase in nearly all the observation points (Fig. 12).Although with a more or less constant nitrate load (nitrate leaching) depending on the scenario, the groundwater nitrate concentrations increase over time, with the highest concentrations for the MLTII scenario.Given the absence of groundwater nitrate attenuation processes, the maintenance of the nitrate leaching load causes an increase in nitrate concentrations over time.This shows that, in order to reduce groundwater concentrations, a significant reduction in fertilizer applications and better management practices would be necessary.All climate change and land use change scenarios show the same trend at a particular point.Regarding climate change, the CNRM scenarios appear as the ones with higher concentrations in most of the case study area.Land use change scenarios do not show in general a significant effect on groundwater nitrate concentration.

Discussion and conclusions
Global (climate and land use) change impacts on water quantity and quality in the Mancha Oriental aquifer system (MOS) have been assessed through the sequential coupling of three models: a watershed hydrological model (SWAT), a groundwater flow model (MODFLOW) and a groundwater transport model (MT3DMS).The models have been successfully calibrated and validated using a large array of sources: gauged streamflow time series, previous assessments of groundwater recharge, crop yields and applied irrigated water (essential for ensuring adequate modeling of crop responses in an essentially agricultural basin), and records of groundwater heads and nitrate concentrations in the Mancha Oriental system (MOS).Spatial and temporal variability of water quantity and quality has been assessed for combinations of three climate change scenarios (corresponding to different GCMs) and three land use change scenarios for nearly all the twenty-first century.In this way, the study has shown that the integrated sequential use of the three models offers a valuable tool for assessing the impacts of land use and climate change pressures.
The modeling results for all scenario runs point out the same future trend: a certain reduction of groundwater recharge, which will be further exacerbated in the long term (except for the HADCM3 scenarios).The combined analysis with the land use change (LUC) scenarios implicitly introduces variations in the amount of groundwater pumping as well.In the end, groundwater levels will be affected by the combined effect of recharge and pumping changes.In this sense, the results show cycles of decreasing and increasing groundwater levels, but without a generalized descending trend.Although climate change seems to be the main driver of the changes in groundwater recharge (the differences among climate scenarios with the same land uses are larger than those introduced by the land use scenarios), groundwater levels seem to be more affected by the LUC because of the effect of changes in pumping.The results also show an increase in groundwater nitrate concentrations driven by the maintenance of nitrate leaching loads.For that, all the scenarios analyzed with MT3DMS showed increasing concentrations.The increased irrigation scenarios were the ones presenting the highest nitrate concentrations.
Although further studies at the basin level would be needed, it can be concluded that global change could become a great threat and challenge to the Mancha Oriental system and, given the importance of stream-aquifer interaction, to the whole Jucar River basin.This will increase the current challenges for aquifer management, already officially classified as in poor status, for reaching the required quantitative and chemical standards of the EU WFD and guaranteeing a sustainable exploitation.
The new PHJ already considers a program of measures designed to prevent further groundwater level depletion, enhance groundwater-head recovery, and revert the increasing nitrate concentrations.
Finally, some reflections in light of the analysis conducted for the case study are -Global change must be taken into account in the planning and assessment of management policies for a sustainable development of the Mancha Oriental system, as the expected results without considering global change impacts are likely to be too optimistic.The currently intended program of measures regarding the MOS might have to be redefined under the expected climate and LUC conditions, searching for robust solutions to adapt to global change conditions.The integrated modeling framework developed in this work could be useful for the assessment of the potential effect of adaptation measures to cope with future global change.
-One essential issue when dealing with assessing future global scenarios is uncertainty.Uncertainty arises from the climate and land use projections, as well as from the impact modeling exercise.Uncertainty on climate and land use projections was partially addressed by considering a set of scenarios, showing the dispersion of results in terms of groundwater recharge, and groundwater quantity and quality.Uncertainty on the impact modeling exercise is quite complex to address, since there are many factors involved (data uncertainty, uncertainty on the model conceptualization, on the model calibration, on the coupling of models, etc.).In our case, we just assessed the robustness of the SWAT model results to the calibrated parameters with a sensitivity analysis.Further research would be need on addressing uncertainty in such a complex framework in a holistic and integrated way.
-Climate and land use changes impact not only groundwater levels and stream-aquifer interaction but also groundwater quality.Increasing groundwater nitrate concentrations can be anticipated, due to the continuous intense use of fertilizers in agriculture over time.In this context, economic instruments might have an essential role in enhancing a sustainable management of diffuse pollution for the future (Peña-Haro et al., 2014).However, water quality projections of this study are subject to a high uncertainty, since no formal calibration or validation regarding nitrate leaching was possible.Further improvements of the modeling chain for nitrate leaching, and collection of additional data and field experiments, would be required in order to obtain more accurate assessments.Although some studies have shown that the leaching rate may increase under future climate scenarios, the implications of climate and land use change on the rate of nitrate leaching are not yet fully understood, given the complex mix of factors that simultaneously affects the leaching process and the usually lack of enough site-specific data (Stuart et al., 2011).
-Participatory processes engaging the relevant stakeholders are essential in the successful definition and implementation of sustainable adaptation measures for groundwater management.In this sense, techniques such as the multi-attribute value theory already applied to the case study (Apperl et al., 2015) can be useful for ranking measures based on the stakeholder preferences and values and anticipating potential conflicts among competing uses.

Figure 2 .
Figure 2. Monthly mean and standard deviation comparison on maximum temperatures at Albacete-Los Llanos weather station and total precipitation on the MOS.

Figure 3 .
Figure 3. Temperature and precipitation 10-year moving average values for the climate change scenarios.

Figure 4 .
Figure 4. CORINE Land Cover images for years 1990 and 2006 in the Mancha Oriental aquifer.

Figure 9 .
Figure 9. Groundwater recharge results obtained from SWAT for the short, medium, and long term.

Figure 10 .
Figure 10.Groundwater level simulation over the twenty-first century of layers 2 and 6 of the Mancha Oriental aquifer.

Figure 11 .
Figure 11.Nitrate leaching results obtained from SWAT for the short, medium and long term.

Figure 12 .
Figure 12.Nitrate concentration simulation in the Mancha Oriental aquifer obtained from MT3DMS.

Table 1 .
Selected driving forces for the Mancha Oriental land use change.

Table 2 .
Crop irrigation, crop and nitrate leaching management calibration.

Table 4 .
Monthly average temperature and precipitation over the short, medium and long term.