Assessing the response of groundwater quantity and travel time distribution to 1 . 5 , 2 and 3 degrees global warming in a mesoscale central German basin

Groundwater is the biggest single source of high-quality fresh water worldwide, which is also continuously threatened by the changing climate. This paper is designed to investigate the response of regional groundwater system to the climate change under three global warming levels (1.5, 2, and 3 ◦C) in a central German basin (Nägelstedt). This investigation is conducted by deploying an integrated modeling workflow that consists of a mesoscale Hydrologic Model (mHM) and a fullydistributed groundwater model OpenGeoSys (OGS). mHM is forced by five general circulation models under three represen5 tative concentration pathways. The diffuse recharges estimated by mHM are used as outer forcings of the OGS groundwater model to compute changes in groundwater levels and travel time distributions. Simulation results indicate that under future climate scenarios, groundwater recharges and levels are expected to increase slightly. Meanwhile, the mean travel time is expected to decrease compared to the historical average. However, the ensemble simulations do not all agree on the sign of relative change. The ensemble simulations do not show a systematic relationship between the predicted change and the warming level, 10 but they indicate an increased variability in predicted changes with the enhanced warming level from 1.5 to 3 ◦C. This study indicates that a higher warming level may introduce more uncertain and extreme events for the studied regional groundwater system.


Introduction
The availability, sustainability, and quality of water resources are threatened by many sources, among which the changing climate plays a critical part (Stocker, 2014).A significant sign of climate change is the global warming, which has been evidenced by the analysis of long-term air temperature records.Not only the earth surface temperature shows a constant warming trend, the sea surface temperature has also increased (Stocker, 2014).There has been adequate proofs that the massive greenhouse gas emissions since the eighteenth century accelerate the global warming process (Stocker, 2014).Consequently, it is urgently needed to estimate the change of meteorological variables (e.g., precipitation and temperature) in the future global warming scenarios.General circulation models (GCMs) combined with different emission scenarios or representative concentration pathways (RCPs) have been widely employed for climate impact study (Collins et al., 2013;Thober et al., 2018;Marx et al., 2018).
Climate change may significantly alter the pattern of terrestrial hydrological processes, influence the spatial and temporal behavior of shallow water storages, and manipulate the degree and frequency of extreme events such as floods and droughts (Van Roosmalen et al., 2009;Sridhar et al., 2017;Thober et al., 2018;Marx et al., 2018).Hydrological processes and states (e.g., evapotranspiration, soil moisture, and potential recharge) are tightly coupled with current climate and meteorological variables (e.g., precipitation, humidity, atmosphere temperature).The impact of climate change on the terrestrial water cycle is, unfortunately, uncertain.Climate model projections show a good consistency in future global averaged trends, but may disagree on the magnitude of regional-scale variables, particularly when precipitation projection is involved in (Meehl et al., 2007).Many past studies devote to estimate the control and uncertainty of climate change on hydrological states and fluxes (Hunt et al., 2013;Samaniego et al., 2018;Renée Brooks et al., 2010;Hattermann et al., 2017;Goderniaux et al., 2009).Among them, some studies indicate that the frequency and intensity of extreme events (e.g., soil moisture drought, heat wave) may be exacerbated owing to anthropogenic warming (Samaniego et al., 2018;Kang and Eltahir, 2018;Marx et al., 2018).The global water scarcity is likely to be exacerbated due to the potential decline in fresh water resources under the 2 • C global warming scenario (Schewe et al., 2014).
As the single biggest source of world's fresh water supply, groundwater plays a critical role in the sustainability of terrestrial ecosystem and the environmental consequences of climate variability.Globally, groundwater makes up 35% of the total freshwater withdrawals, constituting approximately 36%, 27% and 42% of water consumption for households, manufacturing, and agriculture, respectively (Döll et al., 2012).Although the general knowledge of scale-dependent hydraulic properties of the subsurface hydrologic systems is still quite limited, they prove to be increasingly influenced by anthropogenic factors (Küsel et al., 2016).Worldwide groundwater system can be affected by climate variability directly by change in recharge or indirectly by change in groundwater abstraction (Taylor et al., 2012).Furthermore, these effects may be adjusted through anthropogenic activities such as land use change.Many recent studies devote to evaluate the impact of climate change on groundwater availability (Engdahl and Maxwell, 2015;Goderniaux et al., 2015;Jackson et al., 2011;Taylor et al., 2012;Van Roosmalen et al., 2009;Stisen et al., 2011;Woldeamlak et al., 2007;Maxwell and Kollet, 2008;Havril et al., 2017).These past studies often use coupled climate-land-surface-subsurface models to investigate the potential response of groundwater storages to the outer forcings under different climate scenarios.Compared with the land surface processes, the groundwater reservoir is less vulnerable to extreme events (Maxwell and Kollet, 2008).The slow response of groundwater to meteorological variability can be explained by the highly dynamic surface water/groundwater interaction, the existence of variably thick unsaturated zone, and the big volume of groundwater storages.Quantification of uncertainty in future water resource projections and travel times (decades to centuries) of regional groundwater system is critically important for regional water sustainability.
Due to the diverse patterns of the terrestrial water cycle in regions under different climate conditions, climate change will have diverse impacts on the groundwater recharge change.Sandström (1995), for instance, found that in Tanzania, a 15% decline in precipitation, without any change in air temperature, will result in a 40-50% decline of groundwater recharge, indicating a potential amplified change of recharge compared to that of precipitation.While some studies found a increasing trend of recharge (Brouyère et al., 2004;Van Roosmalen et al., 2009), others indicate that climate change is likely to lead to decreased recharge rates (Pulido-Velazquez et al., 2015;Woldeamlak et al., 2007;Havril et al., 2017).The changes of recharge, regardless of sign of change, will significantly influence the groundwater levels, and may lead to ecological problems such as the vanishing of wetlands (Havril et al., 2017).Modification of groundwater recharge will control the flow paths and travel times of pollutants, which is critical to the sustainability of regional groundwater system.Moreover, the modification of groundwater recharge can change the age distribution for water in both the vadose zone and the saturated zone, as well as significantly change the composite age distribution (Engdahl and Maxwell, 2015).
Groundwater travel time distribution (TTD) is a robust description of the storage and transport dynamics within aquifers under various external forcings.It has many implications for hydrogeological and environmental studies.For instance, significant time-lags of the streamflow response to external forcings have been observed by multiple studies (Howden et al., 2010;Stewart et al., 2012).Besides, the legacy pollutants in groundwater reservoirs can have a great impact on the total pollutant loads for agricultural catchments (Wang et al., 2016;Van Meter et al., 2017).Groundwater TTD, as a lumped description of the heterogeneous aquifers, sheds light on the assessment of groundwater responses to non-point source contamination subjected to a changing climate and land use (Böhlke, 2002;Engdahl and Maxwell, 2015).

Study Area
As a sub-basin of the Thuringian basin, the Nägelstedt basin is located at central Germany and it has an area of about 850 km 2 (Figure 1).It is a headwater catchment of the Unstrut river.The Unstrut river is a typical, meandering lowland river with only stands for warm temperate, fully humid, and warm summer climate (Kottek et al., 2006).It shows a leeward decreasing trend 10 of areal precipitation and an rising mean air temperature from the eastern Hainich ridge to the Unstrut Valley (Kohlhepp et al., 2017).In the larger Thuringian basin, groundwater has been intensively extracted for domestic, industrial and agricultural uses.
About 70% of the fresh water requirement for Thuringia is satisfied by groundwater (Wechsung et al., 2008).
The extremely fertile soils in the meadows (wet black soil and loess) make the Thuringian basin one of the best agricultural basins in Germany.Approximately 88% of the total land use of this region is regarded as arable land (Wechsung et al., 2008).
At the same time, the proportion of woodland and grassland has fallen sharply, leading to an extreme reduction in biodiversity in these areas (Wechsung et al., 2008).The nitrogen inputs from agriculture fluctuate with time and positions from 5 kg/ha to 31 kg/ha (Wechsung et al., 2008).
The stratigraphy in this area is characterized by a succession of carbonate-siliciclastic alternations.The main aquifer system consists of several sedimentary rocks, including the Middle Keuper (km), the Lower Keuper (ku), the Upper Muschelkalk (mo), the Middle Muschelkalk (mm), and the Lower Muschelkalk (mu) (Seidel, 2003).The Middle Keuper consists of a marly series with gypsum and dolomite, whereas the Lower Keuper is consitituted of grey clays, schieferletten, and dolomitic limestone.The Upper Muschelkalk (Hauptmuschelkalk) is mainly made up of shelly limestone, marl and dolostone.The Middle Muschelkalk consists mainly of evaporites (gypsum, anhydrite and halite), meanwhile the Lower Muschelkalk consists of limestone and marls (Seidel, 2003;McCann, 2008;Jochen et al., 2014).Karstification occurs in the Muschelkalk formation, but has proved to be limited or concentrated in specific zones in this area (Kohlhepp et al., 2017).
The Nägelstedt basin is chosen as the study area for the following reasons: (1) It is a typical agricultural basin where potential non-point source contamination may threaten the sustainability and resilience of groundwater, and (2) the critical zone (CZ) in Nägelstedt basin has been comprehensively investigated using infrastructure platform from the Collaborative Research Center AquaDiva (Küsel et al., 2016;Kohlhepp et al., 2017).

Methodology and Materials
To investigate the impact of different climate change scenarios, we modified the modeling framework originally developed by Thober et al. (2018) and Marx et al. (2018) from EDgE and HOKLIM projects by means of coupling it to a three-dimensional subsurface model.Specifically, we use temperature and precipitation derived from five GCMs under three different RCPs to force the mesoscale Hydrologic Model (mHM), aiming to derive the land surface fluxes and states under different future warming scenarios.The projected recharges from mHM calculations are fed to the groundwater model OpenGeoSys (OGS) for the assessment of groundwater quantity and travel time distributions.

Climate scenarios
We use five General Circulation Models (GFDL-ESM2M, HadGEM2-ES, IPSL-CM5A-LR, MIROC-ESMCHEM, and NorESM1-M) gathered from the Coupled Model Inter-comparison Projects 5 (CMIP5) to provide the climate variables to the mHM model.Temperature and precipitation are derived from these GCMs under three representative concentration pathways (RCPs; RCP2.6, RCP6.0, and RCP8.5), which are available from ISI-MIP project (Warszawski et al., 2014).RCPs are representations of emission scenarios, with RCP2.6, RCP6.0, and RCP8.0 representing low, medium, and high emission scenarios, respectively.This multimodel ensemble approach enables the consideration of uncertainty in climate modeling.Climate variables from GCMs are further downscaled to a 0.5 • C spatial resolution by means of a trend-preserving bias correction approach  (Hempel et al., 2013).The trend-preserving bias correction approach is capable of representing the long-term mean and extremes of catchment state variables (Hempel et al., 2013).The 0.5 degree data is further interpolated into 5 × 5 km 2 grids by means of external drift kriging (EDK) approach.The EDK approach can incorporate altitude effects at sub-grid scale, and has been successfully used in many past studies (Wood et al., 2011;Zink et al., 2017;Thober et al., 2018).
We use the period 1971-2000 to represent current climate conditions because 1991-2000 is the latest decade that are available in the GCM data.The GCM data from this period serves as a baseline scenario for the future projection of climate change.
A time-sampling approach is applied to estimate the time span for different global warming level of 1.5, 2, and 3 • C (James et al., 2017).The five GCMs have different degrees of climate sensitivity due to the different climate projections, therefore providing different meteorological forcings to the mHM model.Specifically, different time periods of 1.5, 2, and 3 • C global warming are estimated by five GCMs under three RCPs (Table 1).We note that some combinations of GCMs and RCPs cannot be identified for the future climate projection before 2099, resulting in a total of 35 GCM/RCP combinations being used in this study (Table 1).

The mesoscale Hydrologic Model (mHM)
The disaggregated meteorological data are used as meteorological forcings of the mesoscale Hydrological Model (mHM) for a daily simulation.mHM is a spatially explicit distributed hydrologic model that applies grid cells as primary hydrologic units, and accounts for multiple hydrological processes including infiltration, surface runoff, evapotranspiration (ET), soil moisture dynamics, snow accumulation and melting, groundwater recharge, and discharge generation.mHM is forced by hourly or daily meteorological forcings (e.g., precipitation and temperature), and uses accessible physical characteristics including soil textural, vegetation, and geological properties to estimate the spatial variability of parameters by means of its unique Multiscale Parameter Regionalization (MPR) technique (Samaniego et al., 2010;Kumar et al., 2013).The MPR technique is capable of coping with fine-scale features because the effective model parameters are regionalized on the basis of the underlining subgrid-scale information using a consistent upscaling algorithm.The mHM simulations have been successfully established across Europe, and the simulated land surface fluxes have been verified by eddy-covariance stations across Germany (Zink et al., 2017).

OpenGeoSys (OGS)
The porous media simulator OpenGeoSys (OGS) is used to simulate regional groundwater flow and transport processes.OGS has been successfully coupled to mHM through a coupling interface -mHM-OGS (Jing et al., 2018b).The coupling interface interpolates the grid-based recharge produced by mHM into the nodal recharge values spreading over the top surface of OGSmesh.In doing so, mHM and OGS are dynamically coupled as a surface-subsurface model such that the potential recharge produced by mHM can be fed to OGS and serves as outer forcing of the groundwater module (Jing et al., 2018b).Specifically in this study, we feed the projected 5 × 5 km 2 recharge from mHM under future climate scenarios to the coupling interface (mHM-OGS) to force the groundwater model.OGS is based on the finite element method (FEM) and solves the partial differential equations (PDEs) of fluid flow by means of linear/non-linear numerical solver.OGS is capable of simulating single processes including saturated zone flow, unsaturated zone flow, and solute transport, as well as coupled processes including saturated/unsaturated flow, multi-phase flow, and reactive transport.Specifically in this study, OGS is used to compute three dimensional saturated zone flow.
Moreover, a Lagrangian particle tracking method -namely random walk particle tracking (RWPT) -is used to track flow pathways and compute travel time distributions of water parcels (Park et al., 2008a, b;Jing et al., 2018a).The RWPT method assumes that the advection process is deterministic, while the diffusion/dispersion processes are modeled stochastically (Park et al., 2008a).The RWPT method has been widely used to account for reactive transport processes and travel times (Park et al., 2008b;Jing et al., 2018a;Engdahl, 2017).

mHM model setup
Fed by the five GCMs, the down-scaled meteorological dataset with a spatial resolution of 5 × 5 km 2 is used as the outer forcing of mHM.The model is set up across Europe using land use dataset and is forced with spatially distributed meteorologic observations obtained from E-OBS dataset (Haylock et al., 2008).Global parameters of mHM are calibrated against discharge observations using the GRDC database.All ensemble simulations are established with the same morphological, land use, and soil type data in order to keep the relevant parameters consistent throughout this study.Furthermore, the mHM model was validated using observations from many gauging stations across Europe with a period 1966-1995 (Marx et al., 2018).The calibration-constrained parameter set is derived and used for groundwater recharge projection.The projected groundwater recharge, with a spatial resolution of 5 × 5 km 2 , is further downscaled to a 250 × 250 m 2 spatial resolution using the bilinear interpolation in the study site for OGS groundwater model.

OGS model setup
A 25-m Digital Elevation Model (DEM) is used to determine the outer bounds of the catchment and the top surface elevation of three-dimensional model domain.A three dimensional stratigraphic mesh is set up on the basis of above information and bore log data from Thuringian State office for the Environment and Geology (TLUG) (Fischer et al., 2015).The mesh consists of 293,041 structured hexagonal elements with a size of 250 m in the x and y direction as well as with a 10 m resolution in the z direction.The parameter zonation approach is used to represent the heterogeneity of hydraulic properties-hydraulic conductivity in this study.The geological zones within the three-dimensional mesh representing Nägelstedt catchment are displayed in Figure 2.Ten different sediment units are delineated based on the stratigraphy in this area, including Middle Keuper (km), Lower Keuper (ku), Upper Muschelkalk 1 (mo1), Upper Muschelkalk 2 (mo2), Middle Muschelkalk 1 (mm1), Middle Muschelkalk 2 (mm2), Lower Muschelkalk 1 (mu1), Lower Muschelkalk 2 (mu2), alluvium, and the uppermost soil layer (Figure 2).The geological unit "alluvium" represents sandy outwash and gravel near streams, whereas "soil" denotes the uppermost soil layer with a depth of 10 m.The post-calibrated values of the hydraulic conductivity in each geological unit obtained from a previous study are assigned to the corresponding geological layers of the mesh   2).
Given that this study is designed to assess the potential response of regional groundwater system to global warming scenarios, a steady-state groundwater system could be assumed.This assumption is made because the future warming level is a long-term average, and on such a temporal scale (decadal), the short-time fluctuations of climate forcings are essentially damped in the regional groundwater system (Maxwell and Kollet, 2008).
The bottom and outer boundaries of model domain are impermeable, and no-flow boundary conditions are assigned onto these geometries.The spatially distributed recharges estimated by mHM under future climate scenarios are mapped onto each grid nodes :::: node of the mesh surface by the model interface (mHM-OGS).Long-term averaged pumping rates are assigned as Neumann boundaries to each production wells :::: well, wherein the pumping rates are obtained from the literature on the basis of long-term historical data (Wechsung et al., 2008).The total long-term averaged pumping rate over the Nägelstedt catchment is 18870 m 3 /day, and it is set constant for all climate scenarios.A fixed head boundary is assigned to the main perennial streams including one mainstream and three tributaries (Figure 1).For the Lagrangian particle tracking model, about 100000 spatially distributed particle tracers are injected through the top surface of mesh.The spatial distribution of particle tracers is consistent with the spatial distribution of simulated diffuse recharges for each climate scenario.(RMSE) ::::: being :::::::: observed.

Results
This section displays the ensemble of simulated changes in groundwater recharges, levels and travel time distributions.For the sake of clarity, we use plus sign to represent simulated values of increases and minus sign to represent decreases.
Generally, calculations indicate that the groundwater recharge rate is expected to be greater than the 1971-2000 average.The increases in groundwater recharge are below 20% in the majority of GCM/RCP realizations, whereas three GCM/RCP realizations suggest a decrease of groundwater recharge in the study area.The simulation under 3 • C warming scenario represents 5 the greatest SD, i.e., the greatest predictive uncertainty.Note that the predictive uncertainty is mainly introduced by the climate projection using various GCM/RCP combinations, given that mHM is the only hydrologic model used in this study and the parameter values are the same for all simulations.
5.2 Impact ::::::: Climate :::::: impact : on groundwater levels Changes of simulated spatially distributed groundwater levels under future climate scenarios using the minimum, median, and 10 maximum projected recharges are shown in Figure 5. Generally, the areas of topographically-driven flow (e.g., slope) appears to be more sensitive to the changes of recharge compared to the lowland plain.Under 1.5 • C warming scenario, simulated groundwater levels using maximum recharges present an increase ranging from 0 to 10 m compared to those under baseline scenario, whereas those using minimum recharges exhibit a slight decrease.Under 2 • C warming scenario, groundwater levels are expected to increase compared to the base case using median and maximum projected recharges, whereas marginal differences can be found in the simulated levels using minimum projected recharge.Under 3 • C warming scenario, the simulated changes in groundwater levels show the highest variation among the three warming scenarios.Simulations using the maximum recharge suggest a significant increase of groundwater level compared to the 1971-2000 historical average, while simulation using the minimum recharge results in a moderate decrease of groundwater levels (up to a decrease of 5 m at northeastern mountain).The positions of monitoring wells are shown in Figure 1b.
indicate an increased variability in groundwater level change following the increased warming level -which can be evidenced by the increased SD values from 1.5 to 3 • C warming level (Figure 6).
Overall, calculations of spatially distributed groundwater levels help to understand more of the response of groundwater quantity to the projected climate change, but they provide little clue on the change in groundwater transport process.A strong spatial variability in changes of groundwater levels reveals the high sensitivity for climate change in mountainous areas and 5 relatively low sensitivity in lowland plain areas.aquifer system, whereby some geological units present very low hydraulic conductivity values (e.g., mm2 and mu2) and therefore, remarkably slow down the movements of particles in these layers.The mean travel time (MTT), which by definition is the mass-weighted average of travel times for all water parcels within the simulated subsurface system, is a typical metrics for characterizing the timescales of catchment storage.The calculated ensemble averages of MTTs for 1.5, 2, and 3 • C warming scenarios do not exhibit notable differences (79.71, 77.15, and 81.85 years, respectively).projected change in TTDs following the increased warming level -which can be demonstrated by the increased SD values from 1.5 to 3 • C warming (Figure 7).

Discussion and conclusions
We systematically explore the response of a regional groundwater flow system to different global warming scenarios by means of the sequential coupling of a computationally-efficient land surface model (mHM) and a physically-based groundwater model (OGS).The results of ensemble simulations manifest that groundwater recharge is likely to increase moderately for all three warming levels.However, the ensemble simulations do not all agree on the sign of relative change.This finding is consistent with a previous finding that low flows are expected to increase slightly in this region under future climate scenarios , considering that baseflow is the main component of low flow and recharge feeds the baseflow (Marx et al., 2018).Similar increasing trend in climate-induced recharge rates has been suggested by many researchers for other regions such as a northern European catchment (Treidel et al., 2012), high plains of USA (Cornaton, 2012), upper Colorado catchment (Tillman et al., 2016) and Snake River basin (Sridhar et al., 2017).Meanwhile, the seasonal pattern of recharge can be significantly modified by climate change (Chen et al., 2018).
The simulated changes in groundwater levels also manifest similar increases for all three warming scenarios, but show a strong spatial variability depending on the local topography and elevation.These changes can be critical to groundwater/surface water interaction because the increase or decrease in groundwater table would modify the dynamics of groundwater discharge into streams (Havril et al., 2017).Rising :: In ::::: areas :::: with ::::::: complex :::::::::: topography :::: and ::::: dense ::::::: drainage :::::::: network, ::::: rising : groundwater level may activate shallow groundwater flow paths and intensify shallow local flow pathways (Figure 9).This way, the mixing behavior of groundwater storage can also change, because the activation of shallow flow paths will lead to a stronger systematic preference for discharging young water (Kaandorp et al., 2018;Jing et al., 2018a).Those above mentioned processes can be properly reproduced by the coupled model mHM-OGS, given that the coupled model is physically based and is capable to explicitly represent the spatial heterogeneity.Moreover, changes of groundwater levels will impact the land surface processes such as evapotranspiration, soil moisture dynamics and overland flow (Kollet and Maxwell, 2008;Huntington and Niswonger, 2012).
The remarkable influence of climate change on the catchment-scale groundwater travel time distributions is of critical importance to the sustainability of the hydrological system.Simulated mean travel times suggest a moderate decline for three warming scenarios, which is not surprising since the travel time of water parcels is directly controlled by the recharge rate.
With simulated weighted-average MTTs being at a centurial scale, climate-induced variations can significantly affect the longterm sustainability of regional groundwater system.Decline in groundwater MTTs will remarkably shorten the life span of non-point source pollutants (e.g., nutrient and pesticide) in groundwater aquifers and change the spatial and temporal distributions of pollutant concentrations within the aquifer system.Given that the nutrient budget of the connected surface water body is linked with groundwater system, water quality of the surface water body in this region (e.g., the Unstrut river) will response accordingly in the future, although with a long delay.This finding is inline with many recent studies, wherein they highlight the importance of legacy nutrients in catchments as a reason for long-term catchment response (Haygarth et al., 2014;Van Meter et al., 2017).
The steady-state nature of simulations is reasonable for the assessment of long-term climate impact on regional groundwater system because 1) it reduce the computational burden, 2) the temporal fluctuations under the future climate cannot be reasonably projected, and 3) high-frequency fluctuations in external forcings have minor influences on long-term travel time distributions (Engdahl, 2017).However, transient behavior can be very important for many cases where the temporal scale is small and the input forcings are highly dynamic.In recent years, the subject of transient behavior of TTDs has become more and more prevalent in groundwater hydrology (Woldeamlak et al., 2007;Cornaton, 2012;Engdahl, 2017).
We note that the results in this study are only suitable for the Nägelstedt site in central Europe.In other regions of Europe, groundwater recharge change induced by global warming may have distinct behaviors than those shown in this study.For example, some studies indicate an decrease in groundwater quantity in Mediterranean regions due to the decrease in projected precipitation (Pulido-Velazquez et al., 2015;Moutahir et al., 2016).Besides, baseflow is also expected to decrease, leading to an potential increase in drought in Mediterranean regions (Marx et al., 2018;Samaniego et al., 2018).
We only consider the direct impact (i.e., impacts exerted through changed precipitations) of climate change on the regional groundwater system.Interactions between the climate and groundwater are exacerbated by land use change, which is mainly exerted by the intensification of irrigated agriculture.In south Australia and southwest U.S., the transition from natural catchments to rain-fed cropland significantly changes the groundwater storage through the increase in recharge (Taylor et al., 2012).
The indirect influence of global warming to groundwater systems has not been considered in this study.This influence can be a dominant factor threatening the local groundwater system for many regions worldwide (Taylor et al., 2012).Future investigations are needed to incorporate both the direct and indirect impacts of global warming on the sustainability of regional groundwater system.
This increased variability indicates an increased possibility of extreme events in groundwater system following the increase in warming level.Therefore, it is still advisable to reduce :::::: restrain : global warming to 1.5 • C and avoid a global warming of 3 • C.

Figure 1 .
Figure 1.Study area and locations of pumping and monitoring wells within the Nägelstedt basin.Panel a) shows the relative position of Nägelstedt basin in the Thuringian basin, and Panel b) shows the locations of pumping and monitoring wells in Nägelstedt basin.

5
moderate flow velocity under natural conditions.The mountains surrounding the Nägelstedt basin drain almost simultaneously into the Unstrut during heavy precipitation events, which in the past led to regular, prolonged flooding of large parts of the floodplains.The topographic elevations of this catchment range from 164 m at the southeastern lowland to 516 m in the Hainich mountainous region.This region is classified as a Cfb climate region based on the Köppen-Geiger classification, where Cfb

Figure 2 .
Figure 2. Geological zonation and three-dimensional mesh for the aquifer system in Nägelstedt basin (Jing et al., 2018b).Panel (a) underlines the spatial pattern of alluvium and soil layers.Panel (b) further displays the zonation of deep geological units.Full names of legends are listed as follows: km -Middle Keuper, ku -Lower Keuper, mo -Upper Muschelkalk, mm -Middle Muschelkalk, mu -Lower Muschelkalk.

Figure 4 .
Figure 4. Projected changes in groundwater recharge rate under three warming scenarios compared to the baseline scenario 1971-2000.Panel a), b), and c) are the scatter plots showing the individual simulation results, and panel d) is the violin plot showing the uncertainty of ensemble simulations.

Figure 5 .
Figure 5. Contour maps of groundwater levels in Nägelstedt catchment.Panel a) shows the long-term average of groundwater levels in the historical period 1971-2000.Panel b) shows the changes in simulated groundwater levels under 1.5 • C, 2 • C, and 3 • C warming scenarios compared to the baseline scenario 1971-2000 using the maximum, median, and minimum projected recharges.

Figure 6 Figure 6 .
Figure6further shows the changes in groundwater levels at several monitoring wells, of which the locations are displayed in Figure1.In general, changes in groundwater levels are induced by the changes in groundwater recharge such that more groundwater recharge results in higher groundwater levels and vice versa.The uncertainty of groundwater level changes increases following the continuous global warming from 1.5 • C to 3 • C, which can be evidenced by an increasing standard deviation of simulated groundwater levels from 2.20 m for 1.5 • C warming to 4.70 m for 3 • C warming (Figure6).The forecasts of groundwater levels present a wide range of ::::: spread : variation associated with the variability of GCMs.Simulated groundwater levels tend to have the largest increase under three global warming scenarios using meteorological forcings from MIROC-ESM-CHEM model.In contrast, simulated groundwater levels using meteorological forcings from NorESM1-M model show minimal changes compared to the baseline scenario.Although differing in magnitude, the changes of groundwater levels for different wells show a consistent trend (either increasing or decreasing) under the same GCM/RCP realization.The simulations show no systematic relationship between the change in groundwater levels and the change in global warming level, but they do

Figure 7 .
Figure 7. Simulated TTDs in Nägelstedt catchment under 1.5 • C, 2 • C, and 3 • C warming scenarios.Panel a) shows the probability density function (PDF) of TTDs for the ensemble simulations.Panels b), c), and d) show the relative changes of mean travel times (MTTs) under future climate scenarios compared to the base case.

5. 4 Figure 9 .
Figure 9. Conceptual graph showing the influence of rising groundwater level on the regional groundwater flow pattern :: in due to climate change.

Table 1 .
Time periods of 1.5, 2, and 3 • C global warming in five GCMs under three RCPs.