Citation for this paper:
de Graaf, I.E.M., van Beek, R.L.P.H., Gleeson, T., Moosdorf, N., Schmitz, O.,
Sutanudjaja, E.H. & Bierkens, M.F.P. (2017). A global-scale two-layer transient
groundwater model: Development and application to groundwater depletion.
Advances in Water Resources, (102), 53-67.
UVicSPACE: Research & Learning Repository
_____________________________________________________________
Faculty of Engineering
Faculty Publications
_____________________________________________________________
A global-scale two-layer transient groundwater model: Development and application
to groundwater depletion
Inge E.M. de Graaf, Rens L.P.H. van Beek, Tom Gleeson, Nils Moosdorf, Oliver
Schmitz, Edwin H. Sutanudjaja, Marc F.P. Bierkens
April 2017
© 2017 The Authors. Published by Elsevier Ltd. This is an open access article under the CC
BY license. (
http://creativecommons.org/licenses/by/4.0/
)
This article was originally published at:
https://doi.org/10.1016/j.advwatres.2017.01.011
ContentslistsavailableatScienceDirect
Advances
in
Water
Resources
journalhomepage:www.elsevier.com/locate/advwatres
A
global-scale
two-layer
transient
groundwater
model:
Development
and
application
to
groundwater
depletion
Inge
E.M.
de
Graaf
a,b,∗,
Rens
L.P.H.
van
Beek
a,
Tom
Gleeson
c,
Nils
Moosdorf
d,
Oliver
Schmitz
a,
Edwin
H.
Sutanudjaja
a,
Marc
F.P.
Bierkens
a,ea Department of Physical Geography, Faculty of Geosciences, Utrecht University, Netherlands b Department Geology and Geological Engineering, Colorado School of Mines, USA c Civil Engineering, University of Victoria, Canada
d Leibniz-Center for Tropical Marine Research (ZMT), Bremen, Germany e Unit Soil and Groundwater Systems, Deltares, Utrecht, Netherlands
a
r
t
i
c
l
e
i
n
f
o
Article history:
Received 14 November 2016 Revised 27 January 2017 Accepted 27 January 2017 Available online 4 February 2017
a
b
s
t
r
a
c
t
Groundwateristheworld’slargestaccessiblesourceoffreshwatertosatisfyhumanwaterneeds. More-over,groundwaterbuffersvariableprecipitationratesovertime,therebyeffectivelysustainingriverflows intimesofdroughtsandevaporationinareaswithshallowwatertables.Inthisstudy,buildingon previ-ouswork,wesimulategroundwaterheadfluctuationsandgroundwaterstoragechangesinbothconfined andunconfinedaquifersystemsusingaglobal-scalehigh-resolution(5 )groundwatermodelbyderiving newestimatesofthedistributionandthicknessofconfininglayers.Inclusionofconfinedaquifersystems (estimated6–20%ofthetotalaquiferarea)improvesestimatesoftimingand amplitudeof groundwa-terhead fluctuations and changesgroundwater flowpaths and groundwater-surface water interaction rates.Groundwaterflowpathswithinconfininglayersareshorterthanpathsintheunderlyingaquifer, whileflowswithintheconfinedaquifercangetdisconnectedfromthelocaldrainagesystemduetothe lowconductivityoftheconfininglayer.Lateralgroundwaterflowsbetweenbasinsaresignificantinthe model,especiallyforareaswith(partially) confinedaquiferswerelongflowpaths crossingcatchment boundariesaresimulated,therebysupporting waterbudgetsofneighboringcatchmentsoraquifer sys-tems.Thedevelopedtwo-layertransientgroundwatermodelisusedtoidentifyhot-spotsofgroundwater depletion.Globalgroundwaterdepletionisestimatedas7013 km3(137 km3y−1)over1960–2010,which isconsistentwithestimatesofpreviousstudies.
© 2017TheAuthors.PublishedbyElsevierLtd. ThisisanopenaccessarticleundertheCCBYlicense.(http://creativecommons.org/licenses/by/4.0/)
1. Introduction
As theworld’s largestaccessiblesource offreshwater, ground-waterplaysavitalroleinsatisfyingthebasicneedsofhuman so-ciety(Gleesonetal., 2016;UNESCO,2009).Itserves asaprimary source of drinking waterand supplies waterfor agricultural and industrial activities. Duringperiods oflow ornorainfall, ground-waterstorageprovidesanaturalbufferagainstwatershortage, pre-servesevaporationinareaswithshallowwatertables,andsustains baseflows to rivers and wetlands, thereby supporting ecosystem habitats and biodiversity (e.g. Bierkens and vanden Hurk, 2007; de Graafet al., 2013; Wada et al., 2014). Moreover, groundwater oftenflows acrosstopographicalandadministrativeboundariesat
∗ Corresponding author.
E-mail address: idegraaf@mines.edu (I.E.M. de Graaf).
considerable rates, supporting water budgets of receiving catch-ments or aquifers (Schallerand Fan, 2009). However, groundwa-terresourcesareincreasinglythreatenedbyexcessiveabstractions, particularlyinirrigatedareaswhereabstractionratesarehighand groundwater is only slowly renewed (Gleeson et al., 2011). The mostdirecteffectsofgroundwater depletion arefallingwater ta-blesandtheirreversiblelossofgroundwater storage.Asa conse-quence, pumpingcosts increase, and baseflowto rivers and wet-landsisreduced, negativelyaffecting ecosystems,andcompaction oftheemptiedporespacemayleadtolandsubsidence.
Despitetheimportance ofgroundwater andtheexplicit treat-mentofgroundwaterrechargeandpumpinginglobalstudies(Döll et al., 2014a; Wada et al., 2010), most global-scale hydrological modelsdonot includea groundwaterflow component. The obvi-ous reason for this missing component is the lack of consistent global-scale hydrogeological information. Such data is neededto obtaina realistic physical representationofthe groundwater
sys-http://dx.doi.org/10.1016/j.advwatres.2017.01.011
tem, allowing forsimulation of groundwater head dynamics and lateralflows(e.g.Fanetal.,2007;SchallerandFan,2009).Another reasonisthatattheusual timescalesconsidered (fromseasonto years),theamountoflateralflowacrosscellboundariesislikelyto besmallcomparedtotheverticalexchangebetweengroundwater andtheoverlyingsoilthrough recharge,evaporationandcapillary rise(Lametal.,2011).Moreover,manyapplicationsdonotrequire lateralflow. Forinstance,ifthepurposeofthe modelingexercise isto estimate the volume of groundwater depletion using a vol-umebasedapproach(Dölletal.,2014b;Pokhreletal.,2012;Wada etal., 2010), orto relategroundwater storagetolandevaporation (Lam etal., 2011),it generallyissufficient tomodelgroundwater asasinglereservoirwithnoconnectionstoneighboringcells.
However, there are a number of applications that do require lateralflow to be explicitly modeled. First, even though the vol-umesofgroundwatertravelingacross cellsis currentlylimited,it becomesmoreandmoreimportantathigherresolutions(Bierkens etal., 2015; Wood et al., 2012). This is certainly the case when groundwaterabstractionsarelargeandincreasing.Second, partic-ularlybelow megacities in deltas andfor irrigated agriculture in developedcountriessuchastheU.S.,groundwaterabstractionsare oftenfrom(semi-)confinedaquifers.Inconfinedaquifersthe spa-tial influence of abstractions is much larger than in unconfined aquifers,likelyexceeding cellboundaries.Third,thereis consider-ableinteraction betweensurfacewaterandgroundwater, even in confinedaquifers,withthelargerrivers.Notaccountingforthis in-teractionmayoverestimategroundwaterdepletion duetothe ne-glectofincreasedcapture(seee.g.Konikow,2011foracritique on theflux-basedapproach).Thisgroundwater-surfacewater interac-tion is difficult to parameterize in a land-surface modelwithout lateralflow.Fourth,itisimportanttoeventuallysurpassthemere estimationofdepletionvolumesandmovetoestimatinghead de-clinesifoneaimstoestimatebywhattime inthefuture ground-waterbecomesunattainableorwhenfallingheads willresults in riversrunning dry duringlow flow periods. Moreover, a ground-waterflow modelwillbemoresuitable toestimatethe effectsof climatevariabilityandchangeon ecologicallyrelevant groundwa-terdepthsandlowflows(Fan,2015).Finally,eventhoughthe vol-umesofwater crossingcellsmay belimited, flowpaths crossing aquiferboundariesorevenlargercatchment boundaries(deGraaf et al., 2015; Schaller and Fan, 2009) can be long and the asso-ciated groundwater age considerable. Groundwater age dates are alikely additional source ofdata to constrain subsurface proper-ties(aquiferdepth andpermeability)by lackof anydirect obser-vationsofthesequantities(Gleesonetal.,2016).Itisbeneficial to usetheseage dates to their full potential ina groundwater flow model.
In this paper we present the results of a two-layer transient groundwaterflow modelthat has beenbuiltto simulate ground-waterhead dynamicsaffected by changes inclimate andhuman wateruse.Inparticular,itismeanttoestimategroundwater deple-tionvolumesandassociatedheaddeclines,therebyaccountingfor groundwater-surfacewaterinteractionandthepresenceof confin-inglayers. The globalmodel isa first steptoward assessinghow long groundwater reserves will last. It also prepares the ground forfuturehyper-resolutionlandsurfacemodelswheretheeffectof groundwaterconvergenceonevaporationandrunoff cannolonger beparameterizedbuthastobemodeledexplicitly(Bierkensetal., 2015).
Previous work on global-scale groundwater models has been doneby Fanetal.(2013) anddeGraafetal.(2015).Thestudyof Fanetal.(2013)producedafirsthigh-resolutionglobal groundwa-ter table map. However, their method doesnot include hydroge-ologicalinformationonaquifers (e.g.depthsandtransmissivities). Inadditionthehydrologicalconnectionbetweengroundwaterand rivers,whichistheprimarydrainageofgroundwaterinhumid
re-gions, is modeled only implicitly. Moreover, their model requires calibrationtoheadobservationsandonlyreturnsthesteadystate headdistribution.ThestudyofdeGraafetal.(2015)presentsthe firsthigh-resolutionglobal-scalegroundwatermodelincluding hy-drogeologicalinformationandaccountingforgroundwater-surface waterinteractions.Lateralgroundwaterflowsforsingleunconfined aquifersaresimulatedatsteady-state.Arelativesimplemethodfor aquifer parameterization is used based on available global data-sets of lithology(Hartmann andMoosdorf, 2012) and permeabil-ity(Gleesonetal.,2011)suchthatthemethodprovidesresultsfor data-poorenvironments.Also,aquiferthicknessisderivedglobally using terrain attributes. The resultsare promising. This isshown bythehighcorrelationbetweenobservedandsimulatedaveraged groundwaterheads,althoughshallowgroundwaterdepthsinlocal alluvialaquifersinmountainvalleys,smallerthanthegrid resolu-tion(<10 km),arenot captured.Also,thestudy’sgreatest lim-itation is that only thesteady-state head distribution isreported andtemporalchangesingroundwaterheadpatternsduetoclimate orhuman impacts are not considered.Also, only a single uncon-finedaquiferisdescribedandvitalinformationontheaccessibility andqualityofthegroundwaterwaterresourceismissing.This in-formationis neededto simulatetheaquifersensitivity tostorage changes due to climate fluctuations or abstractions. Abstractions arepreferentially positioned inconfinedaquifersystems,asthese systemsarelesssensitiveforclimatefluctuationsandgroundwater quality is oftenbetter than fromunconfined systems (Foster and Chilton,2003).
The first objective of this study is therefore to represent the groundwater systemmorerealistically, by includingaquiferswith aconfininglayer,inordertosimulategroundwaterheaddynamics affected by changes inclimate andhuman wateruseadequately. The groundwatersystemis describedin moredetailby including informationonthepresenceofconfininglayers,leadingtothe cat-egorizationofworld’saquifersinconfinedandunconfinedsystems. The improvedphysicalgroundwaterrepresentation isexpectedto lead toa morerealistic estimate of aquifersensitivityto changes instoragebyclimateorhumanimpacts.
The second objective of this study is to provide estimates of groundwaterdepletion andheaddeclinesworldwide.Previous es-timates of groundwater depletion (e.g. Döll et al., 2012; Wada et al., 2010; Konikow, 2011; Pokhrel et al., 2012) vary by an or-der of magnitude. Forexample theflux-based approach of Wada etal.(2010)estimategroundwaterdepletionat283± 40 km3yr−1
2010astheresidualofgrid-basedgroundwaterrechargeand with-drawal. The approach of Wada etal. (2010) overestimates deple-tionasitdoesnotaccountforincreasedcaptureduetodecreased groundwaterdischarge, norforlong-distancesurfacewater trans-ports. The volume-based approach of Konikow (2011) estimates groundwater depletion at145± 39 km3yr−1 between2001and
2008basedondatafortheUSandotherbigaquifersystemsby ex-trapolationofgroundwaterdepletiongloballyusingtheaverage ra-tioofdepletiontoabstractionsobservedintheUS.Ourstudywill be thefirst tosimulatechanges inglobalgroundwater storagein relationtoclimateandgroundwaterpumpingwhileincluding lat-eralgroundwaterflowandgroundwater-surfacewaterinteraction. Includingthesecomponentswillleadto amorerealisticestimate ofgroundwaterdepletion,thatislowerthanearlierflux-based es-timatesandmoresimilartovolume-basedestimates.
2. Methods 2.1. General
In this study a two-layer groundwater model for the terres-trial part of the world was developed (excluding Greenland and Antarctica). Themodel consistsoftwo parts:(1) thehydrological
Fig. 1. A) Original model structure of the hydrological model PCR-GLOBWB. B) The groundwater store (store 3 in PCR-GLOBWB) is replaced by a groundwater model. The groundwater model is forced with net groundwater recharge minus capillary rise, and surface water levels calculated with PCR-GLOBWB. C) A cross-section illustrating the difference between simulated regional-scale groundwater levels and perched water levels. The latter are simulated as interflow in PCR-GLOBWB.
model PCR-GLOBWB (van Beek etal., 2011) and (2) the ground-water model using MODFLOW (McDonald and Harbaugh, 2000). Bothmodels were run at5 arc-minute resolution(approx.10x10 km atthe equator)over the period1960–2010. The groundwater modelwasforcedwithoutputsfromPCR-GLOBWB,specificallyvia groundwaterrechargeandriverdischarges(Fig.1).Abrief descrip-tionofthemodelsandcouplingisgivenhere,amoredetailed de-scriptionisgivenindeGraafetal.(2015).
2.1.1. Globalhydrologicalmodel
The model PCR-GLOBWB simulates hydrological processes in and between two vertically stacked soil stores (maximum thick-ness0.3and1.2 mrespectively)andoneunderlyinggroundwater store.Themodelwasrunat5 resolutionatadailytimestep.For the climate forcing, monthlydata fromCRU TS2.1 (Mitchell and Jones, 2005) wasdownscaled usingERA-40 (Uppala etal., 2005) andERA-interim(Deeetal.,2011) reanalysistoobtaina daily cli-mateforcing(seedeGraafetal.,2013foramoredetailed descrip-tionofthisforcingdata-set).Foreachtimestepandeverygridcell thewaterbalanceofthesoilcolumniscalculatedbasedonthe cli-matic forcingthat imposes precipitation(rain orsnow depending on temperature),referencepotentialevapotranspiration, and tem-perature.Verticalexchangebetweenthesoilandgroundwater oc-cursthroughpercolationandcapillaryrise.Intheoriginal version ofthePCR- GLOBWBmodel nolateralgroundwaterflow is simu-lated.Instead,groundwaterdynamicsaredescribedbymeansofa linearstorage-outflowrelation(Store3inFig.1A).Specificsurface runoff snow-melt, interflow, andbasefloware accumulated along the drainage network that consistsof laterally connected surface waterelementsrepresentingriverchannels,lakes,orreservoirs.A kinematicwave approximationata sub-dailytimestepisusedto routesurfacewatertothebasinoutlet.
PCR-GLOBWB also simulates human water use by calculating ateach time step sectoralwaterdemand (i.e.theamount of wa-terthat wouldbeabstractedifsufficientwaterwereavailable)for irrigation, industrial,domestic, and livestock), the associated wa-ter withdrawals (partitioned into groundwater and surfacewater based on availability and pumping capacity), consumptive water use,andreturnflows. Irrigationissimulatedby addingadditional watertosoilsasafunctionofsoilwaterdeficitneededtoachieve potentialevaporation.Sectorialwaterdemandswereestimated us-ing country statistics and population numbers downscaled to 5 resolution.Toapproximateeconomicdevelopmentoverthemodel period,dataofGrossDomesticProduct,electricity,andhousehold consumptionwereused.Industrialwaterdemandswerekept con-stantover theyear,whiledomestic demandreflectseasonalityin
relationtoairtemperature.Waterrecyclingratiosforindustryand domesticuseareadoptedfromWadaetal.(2014)andwere calcu-latedper countryonthe basis ofGDPandlevel ofeconomic de-velopment.Irrigation demands were calculated assuming optimal crop growth, using maximum crop transpiration, and accounting forbaresoilevaporationandsoilwateravailability.Lossesduring abstractionand transport are calculated using a country specific irrigationefficiencyfactor(used fromSiebertandDöll,2010). Irri-gationevaporationlossesduringapplicationarenotaccountedfor. Also, domestic irrigationis not includedas itis very small com-paredtocrop irrigation. Totalwaterdemandwasdynamically al-locatedtosurfacewater(rivers,lakes,andreservoirs)and depen-dentonactualavailabilityintheresourceandaccountingfor feed-backssuch asreturn flows of unconsumed abstracted water and groundwater-surfacewater interactions (following de Graafetal., 2013).
2.1.2. Globalgroundwatermodel
Atwo-layerMODFLOWmodel(McDonaldandHarbaugh,2000; Schmitzetal.,2009)replacesthelineargroundwaterstoreof PCR-GLOBWB.Aquiferpropertieswereprescribedandincludeconfined andunconfinedaquifersystems(discussedinthenextparagraph). TheMODFLOWmodelwasforcedwithoutputsfromPCR-GLOBWB, specificallytheflowbetweenlayer2and3(Fig.1)consistingofnet recharge(recharge minus capillaryrise) andriver levels. The net rechargeistheinputfortheMODFLOW’srecharge(RCH)package, while the MODFLOW’s river (RIV) and drain (DRN) packagesare usedtoincorporateinteractionsbetweenthegroundwaterand sur-facewater.As PCR-GLOBWB runson aCartesian (orregular) grid (geographic projection)we have to take account of the fact that cellareasandvolumesoftheMODFLOWgriddonotvaryinspace. Thisis done by adjusting rechargeand the storage coefficient in the RCH andBCF packagesto account forvarying cell areas and volumes associated on a geographic grid (see Sutanudjaja et al., 2011 and de Graaf et al., 2015 for details). Apart from the river waterlevels, the boundary conditions of the globalgroundwater modelwere obtainedasfixed heads set equalto globalsea-level (reference level 0). Initial conditions were obtained by a steady state-statesimulationofhydraulicheadwithaveragerechargeand thenwarmingupthemodelbyrunningtheyear1960back-to-back for20yearstoreachdynamicequilibrium.
Threelevelsofgroundwater-surfacewaterinteractionsare dis-tinguished:(1) large rivers, wider than 25 m, (2) smaller rivers with a width smaller than 25 m, and (3) springs and streams higherupinthevalley.
1.For large rivers, interactions are governed by actual ground-water heads and river levels (HRIV [m]). The latter was cal-culated from river discharges (Qchn [m3s−1]) simulated by
PCR-GLOBWB and using channel properties based (channel widthanddepth)basedongeomorphologicalrelationsto bank-full discharge (Qbkfl [m3s−1]) (Lacey, 1930; Manning, 1891;
Savenije, 2003).The RIV-package calculatesthewater flux be-tween the river and groundwater Qriv [m3d−1]. Water
infil-trates theaquiferifthegroundwaterheadliesbelowtheriver head:Qrivisthenpositive.Waterisdrainedfromtheaquiferif groundwaterheadliesabovetheriverhead:Qrivisthen nega-tive.Qrivforthelargerrivers(Qriv.Big[m3d−1])wascalculated
as:
Qriv.Big =
c×
(
HRIV − h)
if h> RBOTc×
(
HRIV − RBOT)
if h≤ RBOT (1)wherec[m2d−1 ]istheriverbedconductancecalculatedfrom
river length (approximated as the diagonal cell length), river bedwidthanddepthandriverbedresistance(assumed1day), h[m] isgroundwater head,andRBOT [m]is theriverbottom calculated for bankfull conditions (taken as a rule of thumb happeningevery1.5year).Surfaceelevation(DEM[m])istaken asthereferencelevelhere.
2. Smaller rivers are defined as drains; water can only leave the groundwater system through the drain when heads get abovethedrainagelevel.Inthiscase,thedrainagelevelswere taken equalto the DEM.The drainage wasthen calculated as Qriv.Small[m3d−1]:
Qriv.Small=
c×
(
DEM − h)
if h> DEM0 if h≤ DEM (2)
3. Qriv.Big and Qriv.Small quantify flow between streams and aquifers and are the main components of the baseflow, Qbf, which is negativewhen wateris drained fromthe aquifer.At the 5 resolution however, the main stream is insufficient to represent local sags, springs, and streams higher up in the mountainous area. We assume that groundwater above the floodplainlevel canbe tappedby localspringsthat are repre-sentedbymeansofalinearstorage-outflowrelationship.From theDEMandthestorageestimatedinS3(inthePCR-GLOBWB run) we calculate the amount of storage ([L], aquifer depth timesporosity).Fromthisfollowsthelocaldrainageflux.Tobe consistent withtheRIV andDRNpackages, thislinear storage termisalsonegativewhenwaterisdrained.TotaldrainageQbf [m3d−1]isthuscalculatedas:
Qb f =
(
Qriv.Big+ Qriv.Small)
−(
J· S 3, f lp)
(3) whereS3,flp [m3]is thegroundwater storageabove the flood-plain(obtainedfromPCR-GLOBWB)andJ[d−1]istherecession coefficient parameterized based on Kraaijenhof van der Leur (1958):J =
π
T 4 SyL2(4)
where T [m2d−1] is transmissivity,Sy [-] is thespecific yield
(for each hydrolithologicalunit; Table 1 and L [m]is the av-erage distance between streams and rivers as obtained from thedrainagedensity(vanBeek etal., 2011). Globallythe local drainagefluxsumsupto50%ofthetotalbaseflow.
2.2.Modelruns
PCR-GLOBWB andMODFLOWwereonlycoupledone-way.This means we first run PCR-GLOBWB at a daily timestep for 1960– 2010, andforce MODFLOW withmonthly averaged values ofthe
followingoutputs:surfacewaterlevelsarepassedtotheRIV pack-age, net recharge (groundwater recharge minus capillary rise) is passed to the RCH package, and groundwater abstractions are passedtotheWELLpackage.As wedonot knowtheactual loca-tionsofthepumpingwells,wesimulateabstractionswitha pump-ing wellpositioned inthemiddleofa MODFLOWcell (that coin-cideswithaPCR-GLOBWBcellwithgroundwaterabstraction).Also thelinearstorageterm(i.e.J· S3,flp)ispassedtotheWELL pack-age.After that,MODFLOWwasrunover theperiod1960–2010at monthlytimestepswiththesefluxesandboundaryconditions im-posed.
Thistype of couplingensures that the sameamount ofwater passesthroughthegroundwatermodelasthroughthe groundwa-terzone(Fig.1Astore3)ofPCR-GLOBWBandalsoyieldsthesame global average flux between surface and groundwater. Thus, the couplingensures fullterrestrialbalance closure.Also, itimplicitly includes capillary rise, which is calculated in PCR-GLOBWB (van Beek et al., 2011) andis included in the net recharge. However, itdoesnot includefeedback effectsofgroundwater levels onthe portioningofsurfacefluxesasa fullcouplingwoulddo,nordoes itexplicitlyportraytheeffectsofgroundwaterabstractionson sur-facewaterlevels.
2.3. Aquiferparameterization
Theaquiferparameterizationisbasedonavailable global data-sets on lithology (Global Lithology Map (GLiM), (Hartmann and Moosdorf, 2012) and permeability (Gleeson et al., 2014; 2011), combinedwithan estimateofaquiferthicknesses(de Graafetal., 2015). A detailed description of the aquifer parameterization is given in(de Graafet al., 2015), a summary is givenbelow since thefocusofthisstudyistoextendtheaquiferparameterizationby categorizingworld’saquifersinconfinedandunconfinedsystems. 2.3.1. General:aquiferpermeabilityandtransmissivity
Aquifer permeabilities were defined by lumping the lithology classes defined by Hartmann and Moosdorf (2012) to 7 com-bined hydrolithologies(adoptedfromGleeson etal.,2011), repre-senting broadlithological categories withsimilar hydrogeological characteristics. These units were subdivided further on the basis oftextureforunconsolidatedandconsolidatedsedimentary rocks (Table 1), resulting in a global map representing regional scale permeability. The polygons ofthe lithologicalmap, delineating a hydrogeological unit, were subsequently gridded to 30 (approx. 1x1 km at the equator) and aggregated as the geometric mean at 5 resolution (approx. 10x10 km at the equator). To calculate transmissivities(T[m2d−1])aquiferthicknessesare needed.These
thicknesses were estimated since no global data-set of observed aquifer thicknessis available. Using the assumption that produc-tive aquiferscoincidewithsedimentarybasinsandsediments be-low river valleys, the distinction is made between (1) mountain ranges,and(2) sedimentbasinsrepresentingthe aquifers. Subse-quently, aquifer thicknesses were estimated by relating these to terrainattributes(e.g.curvature)aftercalibrationoftheserelations withreportedaquiferthicknessesfromU.S.groundwatermodeling studies(seedeGraafetal.,2015foranextensivedescriptionofthis procedureto delineateaquifersystems andestimatethicknesses). Theseestimatedthicknessesrepresenttheproductiveaquifersuntil thefirstimpermeablebasis,meaningthatpermeableand(thinner) confiningunitsarelumpedtogetherintoonebigaquifersystem. 2.3.2. Delineationofconfininglayers
In thisstudywe delineateconfining layers andcategorize the aquifers of the world in confined and unconfined aquifers. The categorization is done using information on grain sizes of un-consolidated sediments in the lithologicalmap (GLiM)and
addi-Table 1
Lithologic and hydrolihologic categories.
Lithologic categories a Hydrolithologic categories b log k μgeo [m 2 ] b σ[m 2 ] b Sy [m/m] c
Unconsolidated sediments Unconsolidated (mixed or unmixed) −13.0 2 .0 0 .235
c.g. unconsolidated −10.9 1 .2 0 .360
f.g. unconsolidated −14.0 1 .8 0 .110
Siliciclastic sediments Siliciclastic sedimentary −15.2 2 .5 0 .055
c.g. siliciclastic sedimentary −12.5 0 .9 0 .100
f.g. siliciclastic sedimentary −16.5 1 .7 0 .010
Mixed sedimentary rocks Carbonate −11.8 1 .5 0 .140
Carbonate sedimentary rocks Evaporites
Acid volcanic rocks Crystalline −14.1 1 .5 0 .010
Intermediate volcanic rocks Basic volcanic rocks
Acid plutonic rocks Volcanic −12.5 1 .8 0 .050
Intermediate plutonic rocks Basic plutonc rocks pyroclastics metamorphic
water bodies Not assigned – – –
Ice and Glaciers
a Hartmann and Moosdorf (2012) bGleeson et al. (2011) , log k μ
geo is the geometric mean logarithmic permeability; σis the standard deviation; f.g. and c.g. are fine-grained and coarse-grained, respectively.
c S
y is the storage coefficient, average per category.
Fig. 2. A) Map of defined confining layers for unconsolidated sediments. B) Insets showing the quality of the input data used for the lithological map in more detail. Country and state borders are clearly visible. C) Histogram showing the global percentages of confining layers assumed for the minimum and maximum scenario, and per continent.
tional informationin the sedimentproperty description ofGLiM. Fig. 2A shows the resulting map of confining layers. Some clear spatialinconsistenciescan be seen, e.g. fortheUS overthe High Plainaquifer,ortheGerman-Polishborder(inconsistenciesare dis-cussed in Hartmann and Moosdorf (2012)). We interpreted the GLiM lithologies of ‘fine-grained unconsolidated sediments’ and ‘mixed-grainedunconsolidatedsediments’aspotentialnear-surface confining units. About 30% of the mixedgrained unconsolidated sedimentshavelayers ofclay,siltortillinthesediments descrip-tion. We assume, based on the uncertainty in the sediment
de-scriptions,thatthis30%representstheminimumpercentageof un-consolidatedmixedgrainedsedimentswithaconfiningunit.
Additionally, we classified regions of the world where we know confining layers are present but are not fully represented in GLiM; the coastal zones. We reconstructed the inland extent ofHolocenefine-grainedcoastalplainsthat overlayolder,coarser sedimentsfromthePleistocene.A reconstructionofsea-levelrise (e.g.ToscanoandMacintyre,2003)andcoastalsedimentation(e.g Syvitski et al., 2005) implies that 6000 years ago along many Holocenecoastalplainsthecoastlinewasmoreseawardscompared
Fig. 3. A) 2-dimensional schematization of coastal confined aquifer classification. The star indicates the apex. B) Spatial extent of a coastal confined aquifer.
tothe modern situation.Coastlines haverecededtens of kilome-ters since in many deltas (e.g. Stanley and Warne, 1994). With sea-levelrisingovertheHolocene,thedownstreamgradientalong theriverflattened,andsedimentsaccumulatedinthedownstream part.Theapexofthedeltaanddepositionareasindicatethepoint whereaccumulation originated.Thisaccumulationpoint is recon-structedusingtheprofilecurvaturealongtheriver,withits occur-rencewherethecurvaturechanges fromconvextoconcave start-ingfromthecoast. Theaccumulationpointswereidentified glob-allyforlargerrivers(>25 mwide)withtheir outleton,ornear thecoast(i.e.within50 km).Pointsathighelevations,clearlynot relatedtosedimentationinthecoastalzone,wereexcluded.
Fig. 3A shows a schematic cross-section of a coastal aquifer witha confininglayer ontop.Thickness oftheconfininglayer at thecoastisestimatedusingoceanbathymetry(elevationtakenone 5 grid-cell outof thecoast).Thethicknessatthecoastline is in-terpolatedalong theriverusingtheriver gradientuntilthe apex, wherethicknessisminimal.Theinterpolationisdoneforall iden-tifiedlarge riversandthicknessofthecoastalplainsisderivedby interpolation between the rivers, bounded by the elevation con-tours onwhich theapex are situated(Fig. 3B. Followingthis ap-proach,approximately 11% of the globalcoastline is classified as confinedafterinterpolation.
2.3.3. Confinedaquifers:aquifertransmissivityandstorativity Confinedaquifersaredefinedhereasthoseaquifersoverlainby finegrained confininglayers. We realize that theterm ‘confined’ aquifersmaynot entirelycover thedifferent hydrogeological set-tings that are found in reality. First, aquifers covered with fine-grainedsediments maystill allow vertical flow through the con-fininglayersandare thusinfact leakyorsemi-confined aquifers. Second,duetothevaryingthicknessandpermeabilityofthe con-fininglayers certain areas maybe in fact unconfined, depending onthe spatialscaleunder consideration.Third, inareas with ex-cessivegroundwater pumpingthe drawdown close to wells may causepartoftheoriginally(pressurized)groundwaterheadtofall belowthetopoftheaquifer,causingittobe partiallyunconfined. Allthesecasesmeanthatwhatwedefinehereasconfinedaquifers
Fig. 4. Aquifer parameter settings used for unconsolidated sediments for confined aquifer systems. D is thickness [m], D T is the total estimated aquifer thickness, D c
is estimated thickness of the confining layer. For the confining layer horizontal and vertical conductivities ( kh and kv [md −1 ]) were taken similar to fine grained un- consolidated sediments (see Table 1 ), for the confined aquifer kh and kv were taken similar to coarse grained unconsolidated sediments. S is the storage coefficient, de- fined as specific yield (confining layer, values Table 1 ) or storage coefficient (con- fined aquifer).
maybebettertermedaspartiallyandsemi-confinedaquifers. Hav-ing noted this limitation, for the sake of briefness of terminol-ogy, we call all aquifers covered withfine grained confining lay-ers‘confined aquifers’, includingthe hydrogeological settings de-scribedabove.Finally,wenotethatintheimplementationin MOD-FLOWwe useaglobaltwo-layer modelwherefortheareaswith unconfinedaquifersthetoplayerisgiventhesamehydraulic prop-ertiesastheunderlyingaquifer.
The parameter settings used for the unconsolidated confined aquifersandtheconfininglayersareshowninFig.4.Forallother rocks and sediments the values from Table 1 are used. In ab-senceof anyother information,verticalconductivity wasinitially assumedto bethe sameasthe horizontalconductivityandnext, a calibration procedure was used to calibrate the storage coeffi-cient, the horizontal conductivity aswell as the anisotropy ratio kh/kv assuming kh ≥ kv(see Table 1). Aquifer storage capacities were parameterized using specific yields for unconfined aquifers (Table 1) andfor confined aquifers a storage coefficient of 0.001 (Heath, 1983) wasassumed. As mentioned before, the estimated aquifer thicknesses represent productive aquifers until the first impermeable basis. This means that aquifers separated by semi-permeablelayersarelumpedintoonebigaquifersystem.Asa re-sult,headdeclinesmaybeunderestimated.Byincludingconfining layerswherefine-grainedsedimentpropertiesarepresentand us-ing an anisotropy ratio kh/kv (with kh ≥ kv), we expect to im-plicitlycorrect forneglectingsemi-permeable layers andlumping multi-aquifersintoonesystem.
Thickness ofcoastal confining layers is estimatedusing ocean bathymetry andinterpolation (as described in 2.3.2). For 80% of thecoastalconfinedgrid-cellstheestimatedconfininglayer thick-nessisapproximately10%ofthetotalestimatedaquiferthickness. Basedonthisfindingandbylackofanyotherinformation,we de-cidedtoapply athicknessof10%ofthetotal estimatedthickness toallconfininglayersnotlocatednearthecoast.
2.4. Analyzingtheeffectofschematization
We analyzedtheeffects ofincorporatingconfined aquifer sys-tems on simulated groundwater head fluctuations, groundwater flow patterns, groundwater-surface water interactions and stor-agechanges.Threescenarioswere formulateddescribingdifferent spatialdistributions ofconfined and unconfined aquifer systems: (1)unconfined systems only (following de Graafetal., 2015), (2) confined systems for fine grained unconsolidated sediments and mixedgrainedunconsolidatedsediments withfinegrainedlayers
(30%)(in GLiM)andcoastalconfinedregions(minimumconfining scenario),and(3)confinedsystemsforfinegrainedandall mixed grained unconsolidated sediments and coastal confined regions (maximumconfiningscenario).Theparametersettingsofconfined and unconfined systems are used asdescribed above (Fig. 4 and Table1).
2.5. Calibrationandvalidationofsimulatedgroundwaterheads Inordertocalibrateandvalidatethegroundwatermodel, time-serieswereselectedfromtheUS(availablefromtheUSGS:www. usgs.gov/water)andEurope(Rhine–Meusebasin:Sutanudjajaetal., 2011). The used time-series have a record covering at least five yearsandincludeseasonalvariation; fortheUS∼28000stations, forEurope ∼6000 stationswere available. An earliersteady-state version ofthemodel(de Graafetal.,2015)wasalreadyvalidated onall observationsprovided byFanetal.(2013)(givinga totalof 65303cellswithobservations),whichalsoconsistofobservations from Spain, Brazil and Australia. However, many of these obser-vationsare not timeseries butsingleortime-average values.We thereforeusedtheseinourprevioussteady-stateversion(deGraaf etal.,2015)butnottovalidatethepresenttransientsimulation.In casemultiple stationswere positioned inone grid cell,we chose nottoaveragethesebutusethemastheyare,becausethe differ-ent stations generallyhaveobservations over differenttime peri-ods. Also, stations havedifferent surfaceelevations orare placed in differentaquifers(i.e. deeporshallow). We randomlyselected halfoftheobservationsforcalibrationandhalfforvalidation(split sampleapproach).
Limited calibration was used to better fit the groundwater model to the head observations. The calculation times (approx-imately 2 weeks for a 1960–2010 run) were too long to allow for an full-fledged automated calibration (e.g. Doherty, 2015). As afirststep,we ranthegroundwatermodelinsteadystate to cal-ibrate thehorizontal conductivityandthe anisotropy ratioofthe confining layers. A single globalpre-factor was used varying the horizontalandverticalconductivitiesbetween0.01and100times the initial value, changing the anisotropy ratio kh/kv (a total of 9x9 =81runs).Using theparametersofthebestperforming run (the parametersettingwiththeminimumrootmeansquareerror RMSE)we then furthercalibratedthemodelperforming six tran-sientrunswhilechangingthestoragecoefficientbetween0.1and 5 timestheinitial valueandallowing thehorizontal conductivity and anisotropy ratio to be changed +/−10%. This calibration was performedusingthemaximumconfiningconditions,asweseethis asthe mostrealisticscenario. Thebest performing parameterset (with minimum RMSE)found was alsoused using theminimum confiningconditions,andtheresultsofbothrunswerecompared. Themodelperformancewasvalidatedbyevaluating simulated-to observed water table head time-series (the other half of the head observationlocations not usedfor calibration). Heads(with reference to ocean level) instead of depths were used as heads measurepotentialenergyandarethereforemoremeaningfulthan depths. The evaluation focused on mean monthly values, timing of fluctuations (given by the cross-correlation coefficient R), and amplitudes.Thelatteriscompared bythe(relative)inter-quantile rangeerror,QRE,calculatedas:
QRE = IQm 7525− I Q7525o IQo 7525 (5) where,IQm
7525 andIQ7525o aretheinter-quantilerangesofthe
mod-eledandobserveddatatimeseries. Thecross-correlationiscalculatedas:
R= smo
smso
(6)
wheresm andsoare thesamplestandard deviationsofthe mod-eled (m) and observed (o) samples, and smo is the sample co-variance between the two. It was evaluated which stations per-formedgoodagreement intiming(R> 0.5)andfluctuation mag-nitude(QRE < 50%). Results are shownin scatter plots. Addition-ally,wecomparedaveragedsimulatedheadstoaveragedheadsof theobservedtimeseriesandawidersetofobservationsincluding observationsofSpain, Brazil, andAustraliaprovided byFan etal. (2013).
Apartfromtheseglobalvalidationmeasureswealsoevaluated theperformance of themodel inthree areaswheregroundwater depletion is severe (as follows fromthe globalanalysis, sections 2.7and3.5): TheHigh Plainsaquifer(USA),The CentralValleyof California (USA) and India. For these areas, we compared on an aquifer-by-aquiferbasis the averageheaddecline rates(shown in thesupplementaryonlineinformation).
2.6.Groundwaterflowpatterns
The effect of considering confined aquifers on groundwater flow patterns is analyzed by simulating groundwater flow paths. Groundwaterflowpathsandtraveltimesweresimulatedforheads averaged over the simulation period, comparing two model se-tups: unconfined and confined aquifer systems (the unconfined and maximum confining scenarios). Flow paths were calculated withparticletrackinginMODPATH (Pollock, 1994),usingcell- to-cellfluxdensities[volume/time/area]asinput.InMODPATHaflow pathiscomputedbytrackingtheparticlethroughthesubsoilfrom one cell to thenext, from the location ofrecharge(at the water table)towardthepointofextractionordrainage(strongandweak sinks).
FollowingSchallerandFan (2009)thespatialredistribution of localrecharge wasevaluated using the ratio between groundwa-ter drainageandrecharge:Qgw/Rgw. The ratiowasestimatedfor sub-basinsthatwedefinedusingstream-orders,startingfromthe secondstream-orderat5.Thelocaldrainagenetworkwasdefined atthe5grid-cell,andforeverycellastreamispresent.Itfollows that,foraratioof1,allwaterrechargedinthebasinisalsodrained in the basin, or cross basin- boundary groundwater inflows and outflowsareapproximatelyequal.Ifpartofthewaterrechargedin acatchmentflows toa neighboringcatchment anditisnot com-pensated by cross-boundaryinflow, the catchment’s groundwater drainage is less than groundwater recharge, thus Qgw/Rgw < 1. The catchment isthen classified asa ‘netgroundwater exporter’. Ifmorewaterisdrainedthanrechargedinacatchment,the catch-mentisclassifiedasa‘netgroundwaterimporter’andQgw/Rgw> 1.Deviations from1 are primarily a function of geologyand to-pography,whileclimateandbasinsizeinfluencethemagnitudeof thesedeviations(SchallerandFan,2009).Thespatialpatternofnet importers andexporters is expectedto change by including con-finedaquifersystems.
2.7.Analyzingtheglobaleffectsofhumangroundwateruse
Weused thetwo-layer transientglobalgroundwatermodel to analyzetheeffectsofhumangroundwaterwithdrawalson ground-waterheads.The groundwaterabstractionstaken fromthe 1960– 2010 run with PCR-GLOBWB (including human water use) were subsequently used as input for the MODFLOW model through theWELL-package. Forconfinedsystems allabstractions were lo-cated in the confined aquifer, for unconfined systems abstrac-tions were located in the lower layer, which has the same con-ductivityand storagecoefficient as the top layer. The MODFLOW modelwasadditionallyforcedwithPCR-GLOBWB riverlevelsand net rechargeoutputs (recharge minus capillaryrise) that include the effects of abstractions and corresponding return flows over
Table 2
Parameter values used in the calibration process.
Prefactors Parameter values Number of discrete values ‘best’ value fkh ∈ { 10 −2 , .. . , 10 2 } kh = kh ori · f kh 9 khori · 1
fkv ∈ { 10 −2 , .. . , 10 2 } kv = k vori · f kv 9 kvori · 0.1
fSc ∈ { 0 . 1 , 1 , . . . , 4 , 5 } Sc = Sc ori · f Sc 6 Scori · 3
The subscript ‘ori’ refers to the original values used in the model, as presented in Table 1
themodelperiod1960–2010.Wesubsequentlyanalyzed thehead changeovertheperiod1960–2010 fromthesimulations withthe globalgroundwatermodel.ThePCR-GLOBWB runsshowthatover thepastdecades,totalwaterdemandincreasedgloballyfrom ap-proximately 900 km3y−1 for 1960 to 2000 km3y−1 for 2010
(Wadaetal.,2010).Totalgroundwaterabstractionsincreased glob-allyfromapproximately460 km3y−1for1960to980 km3y−1for
2000 (de Graaf etal., 2013). Thesevalues are in thesame range aspreviouslyreportedratesbytheInternationalGroundwater As-sessmentCentre:321 km3y−1 for1960to734 km3y−1 for2000
(www.un-igrac.org).
3. Results
3.1.Delineationofconfininglayers
Fig. 2A shows the spatial distribution of the confining layers. Fortheminimumconfiningscenarioabout6% ofthetotalaquifer areais classifiedasconfined,i.e.belongingto eithercoastal con-fined, or fine grained and layered mixed grained unconsolidated sediments (in GLiM, Hartmann and Moosdorf, 2012). This 6% is assumedtobe theminimum extentofconfined aquiferssystems worldwide. For the maximum confining scenario about 20% of thetotalaquiferareaisclassifiedasconfined,belongingtoeither coastalconfinedorfinegrainedandmixedgrainedunconsolidated sediments.Thisscenarioisassumedtorepresentthemaximum ex-tentofconfinedaquifersystems.Therelativedistributionper con-tinentdoes notdiffer much betweenthetwo scenarios (Fig.2C). Combiningcoastal andnon-coastal confinedaquifers,respectively 5.39x106km2to17.34x106km2oftheworld’saquifersare
clas-sifiedasconfined.
3.2.Calibrationandvalidation
Table2showstheparameterrangesofhorizontal(kh)and ver-ticalhydraulicconductivities(kv)andstoragecoefficients(Sy)used inthecalibrationandtheresulting‘best’parameterset(with min-imum RMSE). It was found that the model performance is best whentheratiobetweenkh andkvis1:0.1andwhenSywas mul-tipliedby 3.With theseparameters the model performance was evaluatedforthethreescenarios(unconfined,minimum,or maxi-mumconfining)focusingonmeanmonthlyvalues,timing of fluc-tuations and amplitudes of groundwater heads. Simulated head fluctuationswere comparedtotheobservedheadtime series(for EuropeandUS)thatarenotusedinthecalibration.Notethat the locationsoftheseheadtime seriesare oftenbiasedtowards river valleys,coastalribbons,andproductiveaquifers.
Scatterplots,plottingthegrid-cellminimum|Qre|(Eq.(5))and maximumR (Eq.(6)), forthe threescenarios are given inFig. 5. Overall thesescattersshow goodagreement in timing (R > 0.5), which(slightly) improves when confined aquifer systemsare in-cluded,seenasmorepointsabovethe0.5-line.Also,overall ampli-tudeerrorsaresmall(|Qre|<50%).Whenincludingconfining lay-erstheamplitudeerrorsincreaseslightly,seen asmorescatterto the50%-lineandabovethisline.Thisslightincreaseinamplitude errorsisparticularlyfoundinareaswithconsiderablerecharge un-der unconfined conditions and/or where groundwater levels are
disconnectedfromthelocaldrainageintheconfinedscenarioand where themodel mostlikely underestimates groundwater heads. Next, we evaluated for how many head observations (averaged overthe gridcell)our simulatedheadsshow goodagreement on the two criteria(clusterI inFig. 5), howmany performgood on one criteria(good onamplitude: clusterII, good ontiming: clus-ter III), and howmany perform bad onboth criteria (clusterIV). Results are given in the table of Fig. 5. For the unconfined sce-nario26%ofthesimulatedheadtime-seriesshowgoodagreement, comparedto29%forthemaximumconfinedscenario.Examplesof time-seriesoftheobservedandsimulatedheadsforwhichtiming andamplitudearewellcapturedaregiveninFig.6.Insteadof plot-tingactualheadvalues,weplottedtheanomalies relatedtotheir meanvalues.Theexamplesshowthatthemodelisabletocapture bothtimingandamplitudeoftheobservedheadsreasonablewell, andthatincludingconfininglayersimprovestheestimates,seenin particularforthe amplitudeestimate inthe CentralValley exam-ple.
The scatter andstatistics (Fig.7) of temporal meansimulated valuesagainst mean observedvalues forthe maximumconfining scenarioshowthataveragedvaluesarereproducedwell(veryhigh R2 and
α
closeto1)andare slightlybetterwhenonlyheadloca-tionsinareaswithsedimentsandsedimentarybasinsare consid-eredinthecomparison.Modelperformanceoftheother scenarios is comparable (see the table of Fig. 5. However, seen from both scattersthe trend is fora generalunderestimation of groundwa-terheads(meaningsimulatedgroundwaterlevelsaredeeperthan observed). This appears especially for the higher elevated areas, whereshallowwatertablesaremostlikelytobesampledbutare notcapturedbyourmodelasaresultofthegridresolution.TheR2
ofsedimentbasinsisthereforeslightlybetterthanwhenmountain areasare includedintheevaluation, howevergroundwater heads insedimentbasinsareunderestimated.
3.3. Globalmapofgroundwaterdepth
Fig. 8 (top) shows the long-term average groundwater depth simulated for the maximum confining scenario simulated under natural conditions (no groundwater pumping). General patterns in groundwater depths can be identified, andare similar for the three scenarios and comparable to previous results by de Graaf et al. (2015). At the global-scale, sea level is the main control of groundwater depth and throughout the entire coastal ribbon shallowgroundwatertablesare found.Theseareasexpand inland where thecoastal ribbon or coastal plainsmeets themajor river deltas,like theMississippi, Indus, andlarge wetlandareasof e.g. South America. At the regional scale, recharge is the main con-troltogether withtopography. Topographyinfluencesforexample theheadsintheflatareasofcentral Amazonandthelowlandsof SouthAmericathatreceivegroundwaterfromthehigherelevated areas. High, flat recharge regions also result in shallow ground-water tables, e.g. the tropical swamps of the Amazon. Regions withlowrechargeratesshowdeepgroundwatertables;thedeserts standoutwheregroundwater,ifpresent,isnowdisconnectedfrom thelocaltopography.Also,formountainregionsdeepgroundwater tablesare simulated.In theseareaslocal aquifersin sedimentary pockets in mountainvalleys are smaller than the grid resolution
Fig. 5. Scatter plots plotting amplitude error (| QRE |) against correlation ( R ) for the three scenarios, the maximum confined scenario shows the best agreement to the observed head time series. When more than one observed time-series was available in a grid-cell, the result of the highest performance value is plotted in the scatter. In the top figure cluster areas are indicated where simulated heads show good performance both on timing and amplitude error (I), good performance of timing (II), good performance on amplitude error (III), and bad performance on both (IV). Percentage of simulated head time-series per cluster for the three scenarios are given in the table, as well as model performance on average groundwater heads.
(<10 km)andarethereforenotcaptured.Asaresult,groundwater heads inthese regions are likelyunderestimated (de Graafetal., 2015).
The differencesingroundwaterdepth betweenthe twomodel layers(top:confininglayer,bottom:confinedaquifer)aresmall. In-sets ofparts ofthe USandEurope illustratethis (Fig.8 bottom). Forunconfined systems,wheretop andbottomlayers haveequal conductivities,anygroundwaterdepthdifferencesareexpectedto be small. For confined systems, where conductivities are differ-ent for the top confining layer and bottom confined aquifer, the
groundwaterheadsintheconfinedaquifers areoftenhigherthan those inthe overlying layers. This is forexample seen along the Dutchcoast,Poriver,Rhone,andalongtheGulfofMexico, includ-ingtheMississippi.However,whengroundwaterrechargedoesnot seep through the confining layer groundwater tables can exceed thehydraulicheadsoftheconfinedaquifer.Thisisseenfor exam-pleintheMississippiembaymentandtheGaronneregion(France).
Fig. 6. Examples of groundwater head anomaly comparison between measured data (red) and simulated data from the different scenarios; unconfined (yellow), minimal confined (green), maximum confined (blue) for a location in the Central Valley (USA) and a location in the Netherlands. We selected these two locations to show the difference between an unconfined (Central Valley) and a confined (Netherlands) system, for two large aquifer systems. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 7. Comparing average observed heads against simulated heads for the maximum confining scenario, for mountain ranges (‘all’, in general deep groundwater levels) and sediment basins (‘sediments’, in general shallow groundwater tables). R 2 is coefficient of determination and αthe regression coefficient of the line passing through the origin.
3.4.Groundwaterflowpatterns
Flowpathssimulationswereperformedusingaveraged ground-waterheads (1960–2010) forthe unconfined andmaximum con-finingscenarios.WefocusthediscussiononpartoftheUS(Fig.9) asthisregion containsshallowanddeepgroundwater tablesand confined and unconfined aquifer systems, and has been the fo-cusofrelatedprevious studies (Gleeson etal.,2011; Schallerand Fan, 2009). Confined systems are e.g. found for the High Plain aquifer, along the Gulf of Mexico, and Mississippi embayment (Fig. 9 left). The simulated flow paths show how groundwater rechargeis redistributedover time andspace(Fig. 9middle and
right).Longerflowpaths,crossingcatchmentsboundaries,provide additionalrechargeto theimporting catchment,thereby support-ing waterbudgets inthat catchment. Thedifference betweenthe unconfined andconfiningscenario(maximum)isevident (Fig.9), e.g. seen forthe coastal lowland aquifer, Mississippi Embayment aquifer, and the High Plain aquifer. Travel distances and travel times are generally shorter if confining layers are present. This resultsfrom thelower permeabilityofthe upperlayer, providing highergroundwatertables,higherdrainagedensitiesandasa con-sequencemorelocalgroundwatersystems.Occasionallyflowpaths endupintheloweraquiferwheretheybecomedisconnectedfrom thesurfacewatersystemoverlongerdistancesresultinginlonger
Fig. 8. Top) Simulated average water table depths (averaged for the period 1960–2010) in meters below the land surface on the upper aquifer under natural conditions (without human water withdrawal) for the maximum confining scenario and the best parameter set. Bottom) head different; heads in bottom minus top layer. For water under pressure in the confined aquifer (bottom layer) a positive value is shown. When heads in the confining layer are higher, negative values are shown.
Fig. 9. Defined confining layers (zoom of Fig. 2 ) and simulated travel times of groundwater flow paths for a part of the US for the unconfined and maximum confined scenario. Flow paths are overlain by major rivers.
Fig. 10. The ratio of groundwater baseflow to recharge (both in mm year −1 ) distinguishing groundwater importing ( Qbf / Rch > 1) and exporting ( Qbf / Rch < 1) sub-basins. A) Frequency of sub-basins found at a given Qbf / Rch ratio: comparing the results of maximal confining and unconfined scenarios. B) Zooms for Qbf / Rch ratio show the spatial distribution. White areas present missing values where recharge is negligible.
traveltimes.Theflow linesinthecaseofunconfinedaquifersare much longer with longer travel times as groundwater tables are generallydeeperandgroundwatersystemsarelargerinsize.
The difference in spatial distribution is also shown by the ratio between groundwater baseflow and groundwater recharge (Qbf/Rch) per a sub-basin (Fig. 10). In this study the sub-basin isdefined by stream-order, starting fromthe second stream- or-der.Aratioof0.1means 90%oftherechargedwaterisexported, a ratio of 2 or larger means at least half of the groundwater drainagecomes from neighboring catchments. The histograms in Fig.10showtheeffectofconfiningunitsonthedistributionofnet groundwater importerand exporter sub-basins.The results show that the changein distribution ofgroundwater exporterand im-portersub-basinsissmallwhenconfininglayersareintroduced. 3.5.Analyzingtheglobaleffectsofhumangroundwateruse
Fig. 11 shows, for the maximum confining scenario, a global mapofthecumulativedepletion(mmpergridcell)of groundwa-ter removedfrom storageby pumping.The well-known hotspots ofgroundwaterdepletion appear, e.g. HighPlains aquifer, Califor-nia,Indusbasin,andSaudi-Arabia(e.g.deGraafetal.,2013;Richey etal.,2015;Scanlonetal.,2012).Asanadditionalmeansof valida-tionwe compared thesimulated headdrops bygroundwater de-pletion withreported rates in a number of well-known stressed aquifers(Gleesonetal.,2012;Richeyetal.,2015):theCentral Val-ley,Highplainsaquifer,andIndia(seeFig.1intheSupplementary Information).Theseresultsshow thatdepletion estimatesin Cali-fornia,theSouthernGreatPlainsandIndiaarereasonabletogood (of similar order of magnitude and spatial pattern), while those in the northern part of the Great Plains are severely over- esti-mated.Possible explanations for thismismatchare an underesti-mationofaquiferpermeabilityandstoragecoefficient, overestima-tionofgroundwater abstractionratesandthe underestimation of enhancedsurfacewaterinfiltration.
Wecomparedthesimulateddepletiontrendsforthemaximum confiningscenarioandunconfinedscenario,calculatedusinghead declinesandstoragecoefficients(Fig.12).Wealsocomparedthese estimatestothedifferencebetweenrechargeandgroundwater
ab-Table 3
Previous published depletion cu- mulative over 1960–20 0 0 (as most studies did not simulate until 2010) compared to this study.
Global total km 3 Pokhrel et al. (2012) a 18 ,960 Wada et al. (2010) a 80 0 0 Wada et al. (2012) a 50 0 0 Konikow (2011) b 20 0 0 Döll et al. (2014b ) a 2240 This study b 2703 a Flux-based b Volume-based
straction (in caseof deficit) per grid-cell, which is the way that depletionrateswereestimatedusingflux-basedmethods(Pokhrel etal., 2012;Wadaetal., 2010).The totaltrendforthemaximum confiningscenarioshowsaglobalincreaseofdepletionwithsome periods of recovery between 1960–2000, after which we see an accelerationof depletion until2010. The recovery periods (1965– 1969,1977–1981)mustbetheresultofclimatevariability-related increases in rechargeand associated increasedcapture from sur-face water, aswe observe an exponential increase inabstraction minusrechargeoverthesameperiod.Thesuddenincrease in de-pletion after 1998 is most likely an interplay between climate (e.g.the1982-83strongElNiño)andincreasedsurfacewaterand groundwaterwithdrawals relatedtosocio-economic trends(Wada etal.,2010). Figure2(SupplementaryInformation)showsthat ef-fectofhydraulicpropertiesonthegroundwaterdepletionvolumes (notevolumes,not heads)isconsiderable,whichmakesestimates ofgroundwaterdepletionbyvolume-basedmethodsrather uncer-tain.
Usingthemaximumconfiningscenarioandthebestparameter set,thetotalgroundwater depletionover 1960–2000isestimated as2703 km3andover1960–2010at7013 km3(resultinginan
av-eragedepletion rateof431 km3 over 2000–2010).The depletion
estimateiscomparabletopreviouspublishedestimatesofvolume andflux-basedapproaches(seeTable3).Fig.12showsthatthe
un-Fig. 11. Cumulative groundwater depletion (in mm over 1960–2010) resulting from human water use for the maximum confining scenario and the best parameter set. For zoons of head declines see Fig. 1 .
Fig. 12. Trend of depletion globally, estimated for the maximum confining and unconfined scenario, and calculated as the cumulative deficit between recharge and abstrac- tion for grid-cells where abstraction is larger than recharge.
confined scenarioshowsingeneralthesametrend.However, the estimatedgroundwaterdepletionislarger.Theestimateddepletion differenceis4803 km3.Thisdifferencecanbeexplainedbythe
in-creaseinrivercaptureundertheconfined scenario.Thepresence of a confining layer results in a larger area of influence ofhead decline.
Despitethelowerpermeabilityoftheconfininglayer,thislarger area ofinfluence causesa largerdecrease inbaseflowora larger increase inriverbed infiltration.Indeed,thegroundwaterdrainage rateoftheconfinedaquifermodelis∼50 km3y−1lowerthanthat
ofthe unconfined aquifermodel, whichadds to3000 km3 over
1960–2010.
The estimated depletion calculated as the deficit between groundwater recharge and abstraction for cells where more wa-terisabstractedthanrecharged(flux-basedmethodasusedine.g. Pokhreletal., 2012;Wada etal., 2010)showsmuchbigger deple-tionvolumesof26,700km3 over1960–2010(withan average
de-pletion rateof 323 km3 over 2000–2010).Again, the difference
in estimated depletion can be explained by the increase in sur-facewatercapture, whichis notaccounted forbycalculating the recharge-abstraction difference but is included in the
groundwa-termodel.Thelarge differenceconfirmstheneedtousea lateral groundwatermodel,accountingforgroundwater-surfacewater in-teractions,whensensitivityofaquiferstostoragechangesis stud-ied.
4. Conclusions
Thispaperpresentsaglobal-scalegroundwatermodel simulat-ing lateralflows andhead fluctuationscaused by changes in cli-mateorhumanwateruseovertheperiod1960–2010.Itisthefirst global-scaletransientgroundwater model that includes a param-eterizationofunconfinedandconfinedsystems(presentedintwo modellayers),andsimulateslateral-flowsandgroundwater-surface waterinteractions inthesesystems. Thisresultsin morerealistic estimates of groundwater head fluctuations and storage changes thanbefore.
One ofthe main contributions ofthisstudyis the parameter-ization of the world‘s aquifer systems, including information on their vertical structure. The model’s aquifer parameterization is basedon globaldata-sets of surfacegeologyand hydraulic prop-ertiesandtopography-basedestimates oftheverticalstructure of theaquifer systems.Onlyglobally available datasets are used in order to keep the methods readily applicable to data-poor envi-ronments.Theworld’saquifersareclassifiedintoconfinedand un-confinedsystemstounderstandaquifersensitivitytogroundwater abstractionsandtoproperly projectfuturegroundwater level de-clines.
Results show that including confining layers (estimated 6– 20% of the total aquifer area) improve estimatesof groundwater head fluctuations and timing and change flow path patters and groundwater-surface water interaction rates. Model performance onaverageheadsandtimingoffluctuationsisslightlybetterwhen confinedsystemsare considered,andbest forthemaximum sce-nario(20%ofthetotalaquiferareaisconfined).Groundwaterflow pathswithinconfininglayersareshortercomparedtopathsinthe aquifer, while flows within the confined aquifer can get discon-nectedfromthelocal drainagesystemduetothe low conductiv-ityof the confininglayer. Thischange in groundwater hydrology is reflected in head fluctuations and flow magnitudes. Although modelperformanceintermsofheadsisonlyslightlybetterwhen confinedaquifersareconsidered,theestimateddepletionratesare closertopreviousvolume-basedestimates(e.g.Konikow,2011)and show that capture of surface-groundwater interaction are simu-latedmorerealisticallythanbefore.
Lateralgroundwaterflowsbetweenbasinsaresignificantinthe model,especially forconfinedaquiferareaswherelong pathsare simulatedcrossing catchmentboundaries,thereby supporting wa-terbudgetsofneighboringcatchmentsoraquifersystems.
The estimated global depletion over the period 1960–2010 is 7013 km3. This value is in the range of earlier published
val-ues(e.g.8000 km3 (Wada etal.,2010),5000 km3 (Wada etal.,
2012)).Hotspotregions ofdepletion arefoundforintensively irri-gatedareas,like theHighPlainsaquifer (simulatedat1390 km3
for2000) or CentralValley (simulated at48 km3 for2000).
Al-thoughdepletion isstronglyoverestimatedfortheNorthernHigh Plainscompared toreported data(∼240 km3 reported for2000
by USGS (Scanlon et al., 2012)), simulated temporal fluctuations mimic measured data well. The same holds for Central Valley, whereweslightlyunderestimatetotaldepletion(whichisreported at ∼53 km3 for 2000). The comparison of depletion volumes
for confined, unconfined, and estimated as the deficit between rechargeandabstraction,showsincreasedrivercapturelowersthe estimateddepletion volume. This confirmsthat a lateral ground-watermodel, withrealistic parameter settings forits aquifers, is preferredwhensensitivityofaquiferstogroundwaterabstractions isstudied.
Ofcoursethereareseverallimitationstoourapproach.Firstly, the grid-resolution (5, approx. 10x10 kmat the equator) is still toocoarse tocapturelocalaquifers inhigherandsteeper terrain. Groundwaterheadsareunderestimatedfortheseareas.Note how-ever,thatperchedwatertableswithinthesoilareaccountedforin thePCR-GLOBWB.
Secondly, our groundwater model is not fully coupled to our hydrologicalmodel,meaningthegroundwater-surfacewater inter-action is still simplified: two-way groundwater-surface water ex-change andthe resultinggroundwater recharge isfirst simulated inPCR-GLOBWB.Thesimulatedgroundwaterrecharge,andthe cal-culatedriverlevelsareimposedontheMODFLOWRIVpackage af-terwards.Althoughthisapproachcanbeexpectedtopreservethe samegroundwater-surfacewaterfluxesinthegroundwatermodel, theeffects ofgroundwater pumpingon thesurfacewatersystem aremoreintricateandwouldpreferablyrequireatightercoupling. Thus,thismodelingstudyisthenextsteptowardsafullycoupled globalsurfacewater- groundwatermodel.Infutureworkthe two models willbe fullycoupledtoincorporate feedbackeffectson a dailytimestepbasis.
Thirdly,themostobviouslimitationisthattheaquifer parame-terizationusedinthisstudyisarelativesimpleapproximationof thetrue3D-hydrogeologicalsetupofaquifersystemsthat contain multiplestackedaquifersandaquitards.Ourmodelthus produces only a first-order estimate of head fluctuations, storage loss and average head declineover the entire aquifer set,most likely un-derestimatingdepletionrateswhenabstractioncomesfromdeeper confinedaquifersystems.Hence,theglobalhydrogeologicalmodel, although at the edge of what can be expected when using only global information, leaves much room for improvement. A con-certed effort of the hydrogeological community to use informa-tionfromregionalonhydrological-and groundwatermodels hold-ing moredetailedinformationofe.g.presence ofconfininglayers orconductivityvalues,aswellasobservationstoupdatetheglobal two-layer modelwouldlikely leadto thelargest improvementin globalgroundwaterassessments(Bierkensetal.,2015).Also,more detailedinformationonthepreciseverticalhydrogeological struc-tureoftheaquifersaroundtheobservationwellscreens,usedfor calibrationandvalidationofthemodel,wouldimprovethe calibra-tion results. Knowledge on the sensitivityof aquifers to changes in climate and human water use is vital and needed to ensure sustainablewateruseparticularlyforthesemi-aridregionsofthe world,wheregroundwaterabstractionswillintensifyduetoan in-crease in drought frequency and duration combined with popu-lation growth,expandingirrigation areas,andrising standardsof living.Thereforeitisimportantnotonlytofocusonlyonareas al-readyexperiencinggroundwaterdepletion,buttoprovideaglobal overviewrevealing where andwhen in futurenew areas experi-encingwaterscarcitywilldevelopinthenearfuture.
Acknowledgments
Thisresearch wasfunded by theNetherlands Organizationfor Scientific Research (NWO) under the program ‘Planetary bound-aries oftheglobalfreshwatercycle’.This workwasdone onthe DutchnationalsupercomputerCartesiuswiththesupportof SURF-sara.Wewouldliketothankseveralcolleagues,amongothersprof. dr.PetraDöll(GoetheUniversity)andprof.YingFan(Rutgers Uni-versity)andtwoanonymousreviewers,fortheirconstructive com-mentsonanearlierversionofthispaper.
Supplementarymaterial
Supplementary material associated with this article can be found,intheonlineversion,at10.1016/j.advwatres.2017.01.011
