A study of Bangladesh's sub-surface water storages using satellite products and data assimilation scheme

2018·5.1 MB·English

by  M. Khaki

ScienceoftheTotalEnvironment625(2018)963–977
ContentslistsavailableatScienceDirect
Science of the Total Environment
journal homepage: www.elsevier.com/locate/scitotenv
A study of Bangladesh’s sub-surface water storages using satellite
products and data assimilation scheme
M.Khakia,*,E.Forootanb,M.Kuhna,J.Awangea,F.Papac,d,C.K.Shume,f
aSchoolofEarthandPlanetarySciences,DisciplineofSpatialSciences,CurtinUniversity,Perth,Australia
bSchoolofEarthandOceanSciences,CardiffUniversity,Cardiff,UK
cLEGOS,UniversitédeToulouse,IRD,CNES,CNRS,UPS,Toulouse31400,France
dIndo-FrenchCellforWaterSciences(IFCWS),IRD-IISc-NIO-IITMJointInternationalLaboratory,Bangalore560012,India
eDivisionofGeodeticScience,SchoolofEarthSciences,TheOhioStateUniversity,Columbus,OH,USA
fInstituteofGeodesyandGeophysics,ChineseAcademyofSciences,Wuhan430077,China
H I G H L I G H T S G R A P H I C A L A B S T R A C T
• WeassimilateGRACEdatatoimprove
a hydrological model estimations
overBangladesh.
• SquareRootAnalysis(SQRA)filteris
usedfordataassimilation.
• Weestimatesub-surfacewaterstor-
agechangeswithinthecountry.
• Climate impacts on the water stor-
agesareinvestigated.
• Independent in-situ measurements
areusedtoevaluatetheresults.
A R T I C L E I N F O A B S T R A C T
Articlehistory: Climatechangecansignificantlyinfluenceterrestrialwaterchangesaroundtheworldparticularlyinplaces
Received25July2017 thathavebeenproventobemorevulnerablesuchasBangladesh.Inthepastfewdecades,climateimpacts,
Receivedinrevisedform21December2017 togetherwiththoseofexcessivehumanwaterusehavechangedthecountry’swateravailabilitystruc-
Accepted23December2017 ture.Inthisstudy,weusemulti-missionremotelysensedmeasurementsalongwithahydrologicalmodel
Availableonline5January2018
toseparatelyanalyzegroundwaterandsoilmoisturevariationsfortheperiod2003–2013,andtheirinter-
actionswithrainfallinBangladesh.Toimprovethemodel’sestimatesofwaterstorages,terrestrialwater
Editor:R.Ludwig.
storage(TWS)dataobtainedfromtheGravityRecoveryAndClimateExperiment(GRACE)satellitemission
areassimilatedintotheWorld-WideWaterResourcesAssessment(W3RA)modelusingtheensemble-based
Keywords: sequentialtechniqueoftheSquareRootAnalysis(SQRA)filter.Weinvestigatethecapabilityofthedata
Bangladesh
assimilationapproachtouseanon-regionalhydrologicalmodelforaregionalcasestudy.Basedonthese
Groundwaterstorage
estimates,weinvestigaterelationshipsbetweenthemodelderivedsub-surfacewaterstoragechangesand
Dataassimilation
Hydrologicalmodelling remotelysensedprecipitations,aswellasaltimetry-derivedriverlevelvariationsinBangladeshbyapply-
GRACE ingtheempiricalmodedecomposition(EMD)method.Alargercorrelationisfoundbetweenriverlevel
heightsandrainfalls(78%onaverage)incomparisontogroundwaterstoragevariationsandrainfalls(57%on
average).Theresultsindicateasignificantdeclineingroundwaterstorage(∼32%reduction)forBangladesh
between2003and2013,whichisequivalenttoanaveragerateof8.73±2.45mm/year.
©2018ElsevierB.V.Allrightsreserved.
1. Introduction
South Asia, and in particular Bangladesh, is amongst the most
*Correspondingauthorat:DepartmentofSpatialSciences,CurtinUniversity,Perth,
water vulnerable regions of the world exhibiting an increase in
Australia.
E-mail address:[email protected](M.Khaki). droughtsandfloodsduetoclimatechange(McCarthyetal.,2001).
https://doi.org/10.1016/j.scitotenv.2017.12.289
0048-9697/©2018ElsevierB.V.Allrightsreserved.
964 M.Khakietal. /ScienceoftheTotalEnvironment625(2018)963–977
Groundwater is the main source of drinking and irrigation water significantroleinbetterunderstandingsurfaceandsub-surfacepro-
(almost 90%) in the country (Islam et al., 2013). Any consider- cessesrelatedtowaterredistributionwithintheEarthsystem(e.g.,
ablechangeinclimatewill,therefore,affectBangladesh’savailable Huntington,2006;Chenetal.,2007;Kuscheetal.,2012;Forootan
water, which is stored in different forms including aquifers, soils, et al., 2014; van Dijk et al., 2014). In particular, Shamsudduha
surface waters as rivers, lakes, man-made reservoirs, wetlands etal.(2012)showedahighcapabilityofGRACEmeasurementsfor
and seasonally inundated areas (Papa et al., 2015). Understanding studying water storage variations in the Bengal Basin. A growing
theinteractionbetweenprecipitation(mainlyprovidedduringthe numberofstudieshavealsoassimilatedGRACETWSinordertocon-
Monsoonseason)andwaterstoragechangesisimportanttorelate strainthemassbalanceofhydrologicalmodels(e.g.,Zaitchiketal.,
climatevariabilitytohydrology.Anin-depthunderstandingofthis 2008;Thomasetal.,2014;vanDijketal.,2014;Eickeretal.,2014;
interactioncanbemoredifficultinBangladeshduetothechanging Tangdamrongsubetal.,2015;Reageretal.,2015;Khakietal.,2017c;
behaviorofmonsoonalprecipitation(WangandDing,2006)aswell Schumacheretal.,2018).
asthelackofknowledgeontheirinfluenceonthehydrologyofthe The present study aims at assimilating GRACE TWS into the
region(Shahid,2010;Rafiuddinetal.,2010). World-Wide Water Resources Assessment (W3RA) hydrological
Groundwater accessibility has made Bangladesh an agro-based model (van Dijk, 2010) to analyze groundwater and soil moisture
country with the main product being rice, making it one of the changeswithinBangladesh.Whilethemainfocusisongroundwa-
world’slargestriceproducer(AbdullahAzizetal.,2015).Theexces- ter and soil moisture, surface water as an important water source
sive groundwater usage during the last two decades has resulted inBangladeshisalsostudiedsincesomesurfacewatersources(e.g.,
in serious problems of both rapid falling of groundwater levels lakes and rivers, except major ones) are not modeled in W3RA.
and the deterioration of its quality (Qureshi et al., 2015). Ground- Moreover, since GRACE TWS reflects the summation of all water
water depletion has been reported by Shamsudduha et al. (2009) compartments, for the first time, we use three different scenarios
between 1985 and 2005 within different regions in Bangladesh to account for surface water storage changes before data assimi-
suchasnorth-central,northwestern,andsouthwesternpartsofthe lation (see details in Section 3.1). The main reason for using the
country. This has also been shown by Shamsudduha et al. (2012) W3RAmodeltoperformourinvestigationsistorelyonthephys-
for the period of 2003 to 2007. Moreover, Sengupta et al. (2013) icalprocessesimplementedinthemodelequationstoconsistently
reportedthatgroundwaterin63(outof64)districtsofBangladesh separate GRACE TWS (since both model and observation errors
areseriouslycontaminatedwitharsenic,whichispartlyattributed are considered) into different water compartments that includes
toitsdepletion.Anumberofstudiesattributethedropingroundwa- groundwaterandsoilmoisture.Ashydrologicalmodelsareusually
terlevelsince1972totherainfalldecreaseandincreaseinhuman better resolved than GRACE data during the assimilation proce-
water usage (see, e.g., Mainuddin, 2002; Ahmed, 2006; McBean et dure, observations are downscaled, and therefore, higher spatial
al.,2011;Deyetal.,2011;Adhikaryetal.,2013).TheGroundwater resolution estimations of water storages will be available within
MonitoringSurveyReportofBangladeshAgriculturalDevelopment the study region (see, e.g., Schumacher et al., 2016). Here, we use
Corporation(BADC)andInstituteofWaterModeling(IWM)showed theensemble-basedsequentialtechniqueoftheSquareRootAnal-
athree-meterdropofgroundwaterlevelsinDhaka(SumonandAbul ysis (SQRA) filtering scheme (Evensen, 2004) to assimilate GRACE
Kalam,2014).Knappettetal.(2016)claimedthatanexcessextrac- TWS into W3RA. SQRA is preferred over the traditional ensemble
tioncausedthegroundwaterleveltodeclinemorethan1mnearthe Kalmanfiltersinceitoffersahighercomputationalspeed,simplicity,
BurigangaRiver,whichpassesinthesouthwestoutskirtsofDhaka andindependenceofobservationperturbations.Besides,Khakietal.
resultingininsufficientresourcesavailablefortherapidlygrowing (2017a)showed that this methodis highlycapable of assimilating
population. GRACETWSdataintoahydrologicalmodel.
Soil water storage variation is another important factor that Afterdataassimilation,weinvestigatetheconnectionsbetween
worsens the situation and affects agriculture. Furthermore, a the estimated groundwater and soil moisture storages (from
considerable amount of surface water from rainfall is consumed improvedmodel)andbothsurfacewaterlevelvariationsandrainfall
by human and thus is not able to recharge the groundwater (e.g., from multi-mission satellite remote sensing data over Bangladesh.
Kanoua and Merkel, 2015; Qureshi et al., 2015; Alimuzzaman, SatelliteradaraltimetryproductsofJason-1and-2,andEnvisatare
2017),whichcanaggravatetheconditionsmentionedabove.Apart usedinthisstudytoprovide19virtualrivergaugestationsforthe
from efforts by these studies, a comprehensive study is missing period 2003 to 2013 distributed across Bangladesh. Since satellite
to account for both groundwater and soil moisture variations and altimetrywasinitiallydesignedforoceanstudies(FuandCazenave,
theirconnectionstoclimatevariabilityandchangeovertheentire 2001),itsobservationsoverinlandwaterbodiesmustbecarefully
Bangladesh. post processed (Birkett, 1998; Calmant et al., 2008; Khaki et al.,
In this regard, hydrological models are important tools for 2015).Therefore,theExtremaRetracking(ExtR)technique,proposed
simulating and predicting sub-surface water storages with high byKhakietal.(2014),isappliedtoretracksatellitewaveformdata
saptio-temporalresolutions(e.g.,WooldridgeandKalma,2001;Döll toimproverangeestimationsandconsequentlyderivebetterwater
et al., 2003; van Dijk et al., 2013). However, imperfect modeling levelestimations.
of complex water cycle processes, data deficiencies on both tem- Further, we apply the statistical method of empirical mode
poral and spatial resolutions (e.g., limited ground-based observa- decomposition (EMD, Chen et al., 2007) to explore connections
tions),anduncertaintiesof(unknown)empiricalmodelparameters, betweenthegroundwaterandsurfacewaterfromthemodel,rainfall
inputsandforcingdatacausesomedegreesofdeficienciesinthem data from the Tropical Rainfall Measuring Mission (TRMM), and
(Vrugt et al., 2013; van Dijk et al., 2011; van Dijk et al., 2014). retracked surface water heights. EMD is an efficient approach to
Theselimitationsareaddressedthroughdataassimilation,whichisa extract cyclic/semi-cyclic components and is preferred over the
techniquethatincorporatesadditionalobservationsintoadynamic classicaltechniquessuchastheFourieranalysis(Chenetal.,2007;
modeltoimproveitsstateestimations(Bertinoetal.,2003;Hoteit Pietrafesaetal.,2016).
etal.,2012).Thetechniquehasbeenwidelyappliedandvalidatedin Theremainderofthisstudyisorganizedasfollows:inSection2,
thefieldsofoceanography,climate,andhydrologicalscience(Garner the study area, and datasets are presented. Section 3 provides
etal.,1999;ElbernandSchmidt,2001;Bennett,2002;Moradkhani a brief overview of the data assimilation filtering methods, the
et al., 2005; van Dijk et al., 2014; Reager et al., 2015). Several ExtR retracking method as well as the EMD approach. Results and
studies indicate that terrestrial water storage (TWS) derived from discussionarepresentedinSection4,andthestudyisconcludedin
the Gravity Recovery And Climate Experiment (GRACE) can play a Section5.
M.Khakietal. /ScienceoftheTotalEnvironment625(2018)963–977 965
2. Studyareaanddata that the full covariance field is well representative of the GRACE
data’serrorstructure.
2.1. Bangladesh Degree1ofStokes’coefficientsarereplacedwiththoseestimated
by Swenson et al. (2008) to account for the change in the Earth’s
Bangladesh is located in the Bengal Basin, where the Ganges, centerofmass.Degree2andorder0(C20)coefficientsarereplaced
Brahmaputra, and the Meghna rivers converge. The average bythosefromSatelliteLaserRangingsolutionsduetounquantified
temperature of the country ranges from 17◦C to 20.6◦C during largeuncertaintiesinthisterm(e.g.,ChengandTapley,2004;Chen
winter and 26.9◦C to 31.1◦C during summer (Rajib et al., 2011). etal.,2007).Colored/correlatednoisesintheL2productsarereduced
Thus, it is placed in the sub-tropical region with a humid, warm, usingtheDDK2smoothingfilterfollowingKuscheetal.(2009).This
andtropicalclimate,whichisdominatedbyasubtropicalmonsoon smoothing causes some degree of signal attenuation (Klees et al.,
originating over the Indian Ocean, which carry warm, moist, and 2008) and moving anomalies from outside the region (e.g., Bay of
unstableair(Ahmed,2006;Khanduetal.,2017).Anaveragedrought Bengal)(Chenetal.,2007;Khakietal.,2018).Tomitigatethisissue,
frequency in the country is reported to be equivalent to 2.5years following Swenson and Wahr (2002), we apply a Lagrange multi-
(Adnan,1993;Hossain,1990)whenrainfall,asthemostimportant plierfiltertodecreaseleakageerrorsovertheentireBangladesh.This
water supply, drops by almost 46% (Dey et al., 2011). The annual filter uses a basin averaging kernel method expanded in terms of
precipitation ranges from less than 1500 to ∼5000mm and varies sphericalharmonicsandsubsequentlycombinedwithL2potential
overdifferentpartsofthecountry,e.g.,1276mmand1337mmin coefficientstoimprovetheGRACEestimates(seedetailsinSwenson
thecentralandwesternregions,respectively(see,e.g.,Hasanetal., and Wahr, 2002). The L2 gravity fields are then converted to 1◦×
2013;Islametal.,2014). 1◦ TWS fields following Wahr et al. (1998). Note that the GRACE
dataprovidechangesinTWSwhileW3RAproducesabsoluteTWS.
Accordingly,themeanTWSforthestudyperiodistakenfromW3RA
2.2. W3RAhydrologicalmodel andisaddedtotheGRACETWSchangetimeseriestoobtainabsolute
valuesandmakethemcomparablewithmodeloutputs(Zaitchiket
The globally distributed 1◦×1◦ World-Wide Water Resources al.,2008).
Assessment (W3RA) model is used to simulate water storage over
Bangladesh. W3RA is a daily grid distributed biophysical model
2.3.2. Satelliteradaraltimetry
developedin2008bytheCommonwealthScientificandIndustrial
SatelliteradaraltimetrydataofJason-1and-2,i.e.,20-Hzsensor
ResearchOrganisation(CSIRO).Themodelsimulateswaterstorage
geographic data records (SGDR), and Envisat, i.e., 18-Hz SGDR
and flows to monitor, represent, and forecast terrestrial water
products are applied in this study. The data includes 260cycles of
storages(vanDijk,2010;Renzulloetal.,2014).Themeteorological
Jason-1covering2002–2008,166cyclesofJason-2covering2008–
forcing data sets for the model include minimum and maximum
2013, and 113cycles of Envisat covering 2002–2012. Jason-2 is a
temperature, downwelling short-wave radiation, and precipitation
follow-on mission of Jason-1 with a similar temporal resolution
fromPrincetonUniversity(seedetailinSheffieldetal.,2006).Effec-
of ∼9.915days and the ground cross-track resolution of ∼280km
tivesoilparameters,includingwaterholdingcapacity,andsoilevap-
(overtheequator),withthesamecharacteristicsasTopex/Poseidon
oration,relatedgreennessandgroundwaterrecession,andsaturated
altimetry mission (Benada, 1997; Papa et al., 2010a). Jason-1 and-
areatocatchmentcharacteristicsarethemodelparameters(vanDijk
2 data are obtained from the Physical Oceanography Distributed
et al., 2013). The model states used in this study include the top,
ActiveArchiveCenter(PO.DAAC)andAVISO,respectively.Addition-
shallow,anddeeprootsoillayers,groundwaterstorage,andsurface
ally,EnvisatRA2productswitha35daysrepeatcycle(30daysfor
water storage in a one-dimensional system (vertical variability).
neworbitafterOctober2010)arederivedfromESA(Table1).
MoredetailedinformationonW3RAcanbefoundinvanDijketal.
Altimeter ranges should be corrected for atmospheric impacts
(2013).
such as ionospheric, tropospheric, and electromagnetic effects
(Benada, 1997). We apply geophysical correction, including solid
earthtide,poletide,anddrytropospheric(Birkett,1995)tocorrect
2.3. Remotelysensedobservations
theranges.TheExtRpost-processingtechnique(Khakietal.,2014)
is applied on waveforms to retrack datasets and improve range
2.3.1. GRACE
measurements.Theretrackedaltimetrydataarethenusedtobuild
TheGRACElevel2(L2)monthlyStokes’coefficientsuptodegree
virtual time series for 19 different points (Fig. 1) located on the
and order 90 and their full error information (2003–2013) are
satellitegroundtracksanddistributedthroughoutthestudyarea.At
obtainedfromtheITSG-Grace2014gravityfieldmodel(Mayer-Gürr
eachvirtualpoint,severalpointsbelongingtothesamesatellitecycle
etal.,2014)andusedinthedataassimilationprocess.Themonthly
are considered, and the median value of the retracked altimetry-
fullerrorinformationoftheStokes’coefficientsisusedtoconstruct
based water levels is computed to address the hooking effect
an observation error covariance matrix for the GRACE TWS fields
(Frappartetal.,2006).Whileasatelliteispassingaboveawaterbody,
(Eicker et al., 2014; Schumacher et al., 2016). Note that different
itsemittedsignalislockedoveraspatiallylimitedpartofthewater,
GRACE products from various centers can lead to different results
whichcanresultinanerror.Thehookingeffect,causedasresultof
dependingontheirdataprocessingstrategies(Shamsudduhaetal.,
off-nadirrangemeasurements,leadstolessaccurateheightestima-
2017). However, for the sake of data assimilation, in addition to
tions (Seyler et al., 2008; Boergens et al., 2016). Afterwards, time
GRACE observations, we also need full error information associ-
series of water level variations from Jason-1 and -2 are combined
atedwiththeobservations.Schumacheretal.(2016)andKhakiet
with those of Envisat products to produce monthly surface levels.
al. (2017b) show that it is important to consider GRACE full error
Details of the datasets, model, and pass numbers of the altimetry
covariancematrixtoconductdataassimilationexperiments.Amore
missionsusedinthisstudyarepresentedinTable1.
comprehensive analysis of different GRACE products has already
beenperformedinarecentlypublishedpaperofSchumacheretal.
(2018). Their results indicate that while using the full covariance 2.3.3. Precipitation
matrixinthedataassimilationprocedure,differencesbetweenthe We use precipitation data of the Tropical Rainfall Measuring
GRACEproductsdonotsignificantlychangetoaffectthefinalresults. Mission Project (TRMM-3B43 products; version 7, (TRMM), 2011;
Therefore,weonlyuseITSG-Grace2014dataforwhichwearesure Huffman and Bolvin, 2012) to assess climate variability over
966 M.Khakietal. /ScienceoftheTotalEnvironment625(2018)963–977
Table1
Asummaryofthedatasetsusedinthisstudy.
Description Platform Detail Dataaccess
Terrestrialwaterstorage GRACE GRACElevel2(L2) https://www.tugraz.at/institute/ifg/downloads/gravity-_field-_models/itsg-_grace2014/
(TWS)
Altimetry-derivedlevel Jason-1 Passnumbers90and231 http://podaac.jpl.nasa.gov
height
Jason-2 Passnumbers90and231 http://avisoftp.cnes.fr/
Envisat Passnumbers337,438,795,896, http://envisat.esa.int/dataproducts/ra2-_mwr/
and982
Precipitation TRMM-3B42 Dailyaccumulatedprecipitation http://disc2.gesdisc.eosdis.nasa.gov/data/TRMM_L3/TRMM_3B42_Daily.7
Hydrologicalmodel W3RA TheCommonwealthScientificand http://www.wenfo.org/wald/data-_software/
IndustrialResearchOrganisation
(CSIRO)
Surfacewaterstorage Satellite-derivedsurfacewater Papaetal.(2015)
storageintheGanges–Brahmaputra
RiverBasin
In-situmeasurements BWDB http://www.ffwc.gov.bd/
Bangladesh. Incorporating more microwave sounding and imagery Bolvin,2012;FlemingandAwange,2013).Thedatasets,validatedby
recordsaswellasimplementingbetterprocessingalgorithmshave Khanduetal.(2017)overthestudyregionshowedpromisingper-
caused a large improvement in this version of data (Huffman and formance.Thegridded(0.25◦×0.25◦)precipitationproducts(2003
Fig.1. Thestudyareaisrepresentedbyblacksolidline.Thefigurealsocontainsthelocationsofvirtualstationsforsatellitealtimetrytimeseriesandvariousin-situstations.
M.Khakietal. /ScienceoftheTotalEnvironment625(2018)963–977 967
to 2013) are converted to 1◦×1◦ and used to investigate their respectively.X¯f,whichrepresentstheforecastensemblemean,and
connectiontowaterstoragechanges. themodelstateforecasterrorcovariance(Pf)arederivedby,
2.4. Surfacestoragedata X¯f = 1 (cid:4)N (X), (3)
N i
For the objective of data assimilation, considering that many i=1
surfacewatersources(indifferentforms,e.g.,lakesandriversexcept
fsehwoumldajboerorenmeso)vaerdenfrootmmGodReAlCeEdiTnWWS3dRaAt,as.uTrofactheiswaetnedr,swtoeraugsees Pf = N−1 1(cid:4)N (cid:2)Xif −X¯f(cid:3)(cid:2)Xif −X¯f(cid:3)T. (4)
satellite-derived surface water data in the Ganges–Brahmaputra i=1
River Basin (as the main source of surface water in Bangladesh)
provided by Papa et al. (2015). The data is based on a multi- The model state contains N different vectors of the model state
satellite approach that combines surface water extent from the variables.NotethatAf = [Af ...Af]istheensembleofanomalies,
1 N
Global Inundation Extent from Multi-Satellite (GIEMS, Papa et al., t(cid:2)he deviation o(cid:3)f model state ensembles from the ensemble mean
2006; Papa et al., 2010a; Prigent et al., 2012) and level height Af =Xf −X¯f . In the second update step, SQRA computes the
i i
variationsofwaterbodiesfromEnvisatradaraltimetrytoestimate ensembleperturbationsthrough.
surfacewaterstorage(Frappartetal.,2012)coveringtheperiodfrom
2003 to 2007. Since the study period is 2003 to 2013, canonical (cid:5)
correlationanalysisisappliedtoextendthedatafrom2007to2013. Aa=AfV I−STSHT, (5)
Satellite derived river height fluctuations of Section 2.3.2 that are
distributedacrossthestudyareaareusedintheprocessofextending where S and V are result from singular value decomposition of Af
thesurfacewaterstorageofPapaetal.(2015).MoredetailsonCCA (Af =USVT).Hisarandomorthogonalmatrix(e.g.,therightsingu-
areprovidedinSection3.4. larvectorsfromasingularvaluedecompositionofarandomN×N
matrix) for ensemble redistribution of the variance reduction (cf.
2.5. In-situmeasurements Evensen,2004;Evensen,2007;Khakietal.,2017a).
To evaluate the performance of data assimilation, 198 ground-
3.1.2. AssimilationofGRACEdata
water stations and 12 soil moisture stations (see Fig. 1) are used.
To assimilate GRACE TWS into the model, we use a summa-
These data are provided by the Bangladesh Water Development
tion of model’s vertical water compartments (e.g., soil moisture,
Board(BWDB)andInstituteofWaterModelling(IWM)inTheAsian
groundwater,andsurfacewater)at13gridpoints.Thissummation
Development Bank (2011). Fig. 2 shows the sample products of
is then updated by the GRACE TWS at the same location at every
different groundwater stations, as well as soil moisture variations
assimilationstep(wheneveranewobservationisavailable).Initial
measured at various depths. Specific yields ranging from 0.01 to
ensemble members are generated by perturbing the meteorolog-
0.2 (Shamsudduha et al., 2011; BWDB, 1994) are used to convert
ical forcing fields following Renzullo et al. (2014). In this regard,
well-waterlevelstovariationsingroundwaterstorage.Detailsofthe
the three most important forcing variables; precipitation, temper-
datasetsusedinthisstudyareoutlinedinTable1.
ature, and radiation are perturbed using Monte Carlo sampling of
multivariate normal distribution (with the errors representing the
standarddeviations)toproduceanensemble(with72membersas
3. Method
suggestedbyOkeetal.,2008).Theperturbedmeteorologicalforcing
datasetsarethenintegratedforwardwiththemodelfrom2000to
3.1. Dataassimilation
2003toprovideasetofstatevectorsatthebeginningofthestudy
period.
3.1.1. Filteringmethod
Two widely used tuning techniques of ensemble inflation and
Thesquarerootanalysis(SQRA)schemefortheensembleKalman
localization are applied to enhance the assimilation performance
filter(EnKF),presentedinEvensen(2004)isusedtoassimilatethe
especiallywhenalimitedensemblesizeisassumed.Ensembleinfla-
GRACE TWS into W3RA. SQRA, which is a deterministic form of
tion uses a small coefficient (i.e., 1.12 in our study; Anderson,
ensemble-based Kalman filters, uses a statistical sample of state
2001)toinflatepriorensembledeviationfromtheensemble-mean
estimatesandunliketraditionalKalmanfilteringmethod,doesnot
to increase their variations and alleviate the inbreeding problem
need an observation perturbation (Burgers et al., 1998; Sakov and
(Anderson et al., 2007). For localization, the Local Analysis (LA)
Oke, 2008; Khaki et al., 2017a). Instead, by introducing a new
scheme (Evensen, 2003; Ott et al., 2004; Khaki et al., 2017b) is
sampling scheme, SQRA uses unperturbed observations without
applied.LArestrictstheimpactofagivenmeasurementintheupdate
imposinganyadditionalapproximationslikeuncorrelatedmeasure-
step to the points located within a certain distance (3◦ following
menterrors(Evensen,2004).TheupdatestageinSQRAincludestwo
Khakietal.,2017b)fromthemeasurementlocation.Wealsoimple-
stepsstartingwithupdatingtheensemble-meanas,
mentthreedifferentcasestodealwithsurfacewaterstorageduring
(cid:2) (cid:3) dataassimilation.
X¯a=X¯f +K y−HX¯f , i=1...N, (1)
• Case1:AssimilatingtheGRACETWSdataafterremovingsur-
facestoragesintothemodelstatesexceptforthesurfacewater
(cid:2) (cid:3)
K=Pf(H)T HPf(H)T+R −1, (2) compartment.
• Case2:Addingsurfacewaterstoragetomodelsurfacewater
compartmentandusingtheGRACETWStoupdatethesumma-
where‘f(cid:4) standsforforecast,‘a(cid:4) foranalysis,andNistheensemble tionsofallwatercompartments.
number. In Eq. (1), X¯a is the mean analysis state, K represent • Case 3: Assimilating the GRACE TWS to update the sum-
the Kalman gain, and y is the observation vector. The transition mations of all water compartments (including surface water
and observation covariance matrices are indicated by H and R, compartment).
968 M.Khakietal. /ScienceoftheTotalEnvironment625(2018)963–977
Fig.2. (a):In-situgroundwaterlevelvariationsofvariousstations(A:21.790◦latitudeand91.870◦longitude,B:24.450◦latitudeand88.598◦longitude,C:26.190◦latitudeand
88.750◦longitude,D:25.590◦latitudeand88.900◦longitude,E:22.110◦latitudeand92.100◦longitude).(b):SoilmoisturevariationsatdefferentdepthsbelongtoRajshahiin-situ
station(24.3670◦latitudeand88.7000◦longitude).
InSection4.1,theresultsofallthecasescenariosarecompared y. The process extracts canonical coefficients u and v such that
witheachotherandevaluatedagainstin-situgroundwatermeasure- X = xTu and Y = yTv (X and Y are canonical variates) possess
ments. a maximum correlation coefficient (Chang et al., 2013) using the
followingfunction,
3.2. Empiricalmodedecomposition(EMD)
E[XY]
Theempiricalmodedecomposition(EMD)proposedbyChenet P= (cid:6) (cid:7)
sqrt E[X2]E[Y2]
al.(2007)isusedforanalyzingmultivariatedatasetsofthisstudy. (cid:8) (cid:9)
EMDestablishesdifferentfrequenciesandtrendswithintimeseries, = (cid:6) (cid:8)E uTxy(cid:9)Tv(cid:8) (cid:9)(cid:7)
which are called Intrinsic Mode Functions (IMFs), by considering sqrt E uTxxTu E vTyyTv
liosctahlatosacislliagtnioanlsis(Fcolamndproisnedetofalf.a,s2t0o0s4c;illRatililoinngs,s2u0p0e3ri)m. pTohseedidbeya = (cid:6) uTCxyv (cid:9)(cid:7), (6)
sqrt uTCxxuvTCyyv
slowoscillations(Flandrinetal.,2004).Thus,EMDdecomposesany
complicated data set into a finite and often small number of IMFs
hiddenintheobservations(Huangetal.,1998).WeapplyEMDon where Cxx and Cyy are covariance matrices of x and y respectively
allavailabletimeseriesofthisstudytoextractdifferentfrequencies andtheobjectiveinabovefunctionistomaximizethecorrelation
and also to find local trends for a better understanding of their P. Once the coefficients are calculated, they can be used to find
interrelationships. the projection of x and y onto u and v as canonical variates with
maximumcorrelation.Here,xcontainsthevectorsofsurfacewater
3.3. Retrackingscheme storagesfromPapaetal.(2015)ateachgridpointandyincludesriver
heightsvariationsfromsatelliteradaraltimetryinasametemporal
Inthisstudy,weusetheretrackingmethodtoimprovealtimetry scaleastheformerdata(2003to2007).Afterperformingcanonical
estimations of river height variations. The retracking process is correlationanalysis,thecomputedcanonicalcoefficientofuandv,
essential since complex waveform patterns are usually observed andanewsetofvariablesy(from2007to2013)areusedtoestimate
overrivers.Tothisend,ExtremaRetracking(ExtR)post-processing thecanonicalvariateofx.Thecombinationofsurfacewaterstorages
technique(Khakietal.,2014,2015)isused.TheExtRisathreestep (x) using the extracted u from the first part has the maximum
filter that starts by applying a moving average filter to reduce the correlationtothealtimetry-derivedriverheightsvariability.Hence,
randomnoiseofthewaveforms.Itthenidentifiesextremumpoints this coefficient vector can be used to transform river heights into
ofthefilteredwaveforms,andfinally,extractsthemainleadingedge surfacewatersateachgridpoint.
amongsttheestablishedextremumpoints.Themethodisappliedto
processdifferenttypesofwaveformsandimprovelevelestimations
asdemonstratedinKhakietal.(2014,2015).Thefilterisemployed 4. Results
here to retrack satellite radar altimetry data and then removing
surfacestoragefromTWS(seeSection3.4).Fig.12showsriverlevel 4.1. Dataassimilation
fluctuationsfordifferentpartsofBangladesh(Fig.3a)andtheentire
areaofthecountry(Fig.3b). Before discussing groundwater and soil moisture variations
withinBangladesh,theeffectofdataassimilationonterrestrialwater
3.4. Canonicalcorrelationanalysis(CCA) storagetimeseriesanditscapabilitytoimprovemodelsimulations
are investigated. Fig. 4 shows average TWS time series within
Canonical correlation analysis (CCA) seeks to find the linear Bangladeshbefore(model-freerun)andafterdataassimilation.The
relationshipbetweentwosetsofmultidimensionalvariablesxand figurealsocontainsGRACETWStimeseries.Itcanbeseenthatdata
M.Khakietal. /ScienceoftheTotalEnvironment625(2018)963–977 969
Fig.3. Averageriverheightvariationtimeseriesfromsatelliteradaraltimetryfordifferentparts(a)andfortheentirearea(b)ofBangladesh.Theaverageerrorforeach
measurementispresentedaserrorbars.
assimilation largely reduces misfits between model-free run and of the in-situ measurements. This is also done for outputs of the
observations. WaterGAPGlobalHydrologyModel(WGHM;moredetailsonDöllet
To assess whether data assimilation (e.g., in Fig. 4) can result al.,2003,MüllerSchmiedetal.,2014),aswellasestimatedstorages
in better water storage estimates, in-situ groundwater and soil byvanDijketal.(2014),indicatedherebyW3,whomergedGRACE
moisture measurements are used for validation. Time series of observationswithanensembleofhydrologicalmodeloutputs.The
groundwater and soil moisture anomalies are generated for each comparison between these products and data assimilation results
station. Groundwater and soil moisture results from all the three allowsustobetterinvestigateanyachievedimprovements.Forthis
assimilationcases(cf.Section3.1.2)arespatiallyinterpolatedusing purpose,theRMSEandcorrelationsbetweenin-situandestimated
thenearestneighbour(theclosestfourdatavalues)tothelocation timeseriesarecalculated.
970 M.Khakietal. /ScienceoftheTotalEnvironment625(2018)963–977
Fig.4. AverageTWSchangetimeseriesfromdataassimilation(case1),model-freerun,andGRACE.
Table2summarizestheaverageRMSEandcorrelationforeach etal.,2017).Here,wecompareW3RA’stoplayerestimationswith
of the three data assimilation case. From Table 2, it can be seen in-situof0–10cm,toplayerplusshallow-rootsimulationswithin-
that the groundwater results are more correlated to in-situ mea- situof0–20cm,andsummationofthetop,shallow,andaportion
surements after the application of each assimilation case (0.81 on of deep-root soil layers with 0–50cm in-situ measurements. Note
average),0.39largerthanmodelsimulationswithoutapplyingdata that WGHM and W3 outputs are provided at a single aggregated
assimilation (model-free run). An average RMSE improvement of layerandcorrespondinglyarecomparedwithin-situsoiltimeseries
51.16% (at 0.95 confidence level) in case 1 shows a significant at the depth 0–50cm. Table 3 shows that the highest correlation
influence of the data assimilation scheme, approximately 4.44% improvements, 18.31% (on average) for all layers and 25.25% for
and 39.11% larger than cases 2 and 3, respectively. It is also evi- the deep layer. Case 2 also represents considerable improvements
dentfromTable2thatdataassimilationresults,especiallycases1 slightly smaller than case 1, still 11.57% larger than case 3, 6.97%
and2outperformgroundwaterestimatesofWGHM.Notethatthe larger than WGHM, and 9.25% larger than W3. Both Tables 2 and
providedW3doesnotincludegroundwaterandthereforeweuseit 3 demonstrate a high capability of data assimilation in improving
for soil moisture comparison only. Table 2 emphasizes that model model simulations of different compartments. These tables also
groundwaterestimationscanbesuccessfullyimprovedwithrespect indicateabetterperformanceoftheimplementeddataassimilation,
to the in-situ measurements if they are fine tuned by GRACE data specificallycases1and2,comparedtoWGHMandW3.
throughtheassimilationespeciallyforcases1and2. To better analyze the differences between each assimilation
Furthermore,correlationanalysisiscarriedoutbetweenin-situ case, we compare their RMSE during 2007. In 2007, a major
soilmoisturemeasurementsatvariousdepthsanddataassimilation flooding(followingENSOrains)occurredacrossSouthAsiaaffecting
resultsfromdifferentscenarios,aswellassoilmoistureestimatesof Bangladesh (Gaiha et al., 2010). This phenomenon can help us to
WGHMandW3(Table3).In-situmeasurementsatdifferentdepths monitor performances of each case in such an extreme situation
are compared with different layers from data assimilation results. and their ability to distribute observed TWS between all water
Forthispurpose,in-situsoilmoisturetimeseriesof0–10cm,0–30 compartments. Groundwater estimates from each case and in-situ
cm, and 0–50cm depths are compared to the model top, shallow, measurements are used to calculate RMSE for each assimilation
and deep soil moisture layers. While the model soil moisture of case (Fig. 5), where the least errors are estimated by cases 1 and
top layer corresponds to the thickness between 5 and 10cm, the 2. Assimilating the GRACE TWS without considering surface water
model shallow and deep-root soil layers broadly represent 10–21 storage within the area (case 3) causes larger errors especially in
cmand3–6msoilthicknesses(seealsoRenzulloetal.,2014;Tian AprilandSeptember.Thelargesterror,however,isobtainedforthe
Table2
Statisticsofgroundwatererrors.Foreachcase,theRMSEaverageanditsrange(±XX)atthe95%confidenceintervalispresented.Improvementsindataassimilationresultsare
calculatedwithrespecttothegroundwaterstoragesfromthemodelwithoutimplementingdataassimilation.
Improvement(%)
Assimilationscenario Correlation RMSE(mm) Correlation RMSE
Case1[RemovedsurfacestoragesfromGRACETWS] 0.86 35±5.65 51.16 57.36
Case2[AddedsurfacestoragestoW3RAsurfacewater] 0.82 39±5.18 48.78 52.92
Case3[Nosurfacestoragecorrectionapplied] 0.75 68±7.72 44.02 18.25
WGHM 0.79 57±5.37 46.83 30.89
Mode-freerun 0.42 83±9.29 – –
M.Khakietal. /ScienceoftheTotalEnvironment625(2018)963–977 971
Table3 showmorenegativegroundwatervariationswithindifferentpartsof
Averagecorrelationsimprovements(at95%confidenceinterval)betweenin-situand Bangladeshthansoilwaterstoragevariations(seeFig.6).Bothwater
soilmoistureestimateswithrespecttomodel-freerun.
compartmentsindicatelargersignals(intermsofamplitude)inthe
Filter 0–10cm 0–20cm 0–50cm centralandnorthwesternpartsofBangladesh.Positivesoilmoisture
Case1 10.42 19.27 25.25 variationsarefoundinthecentertowardeastandnorthwithinthe
Case2 11.10 17.88 24.48 studyperiod,especiallyfortheassimilationandWGHMmaps.Larger
Case3 5.25 8.34 12.91 groundwater variations are also captured in the same area. While
WGHM – – 17.51
assimilation results show negative groundwater changes over the
W3 – – 15.23
central,eastern,andtoalesserdegreesouthernparts,WGHMonly
indicatesnegativevariationsinthesouthernandeasternparts.Fig.7
indicatesthatsmallerwaterstoragevariationsinthenorthwestern
model-freerun.Hereafter,weusetheresultofdataassimilationfor andnortheasternpartsofBangladeshduring2003–2013.
case1sinceitperformedslightlybetterthancase2andsignificantly The average time series of soil moisture and groundwater
betterthancase3intermsoftheRMSE(seeFig.5). storagesfromdataassimilationareshowninFig.8aandb,respec-
The model’s water storage variations computed by assimilat- tively. We estimate spatial averages for all the time series at grid
ing GRACE TWS data into W3RA are presented in Figs. 6 and 7. pointsforthisfigure.Fig.8ashowsslightdeclinesinthesoilwater
Temporal averagesof soilmoistureand groundwaterstoragevari- storage after 2007, which can be related to variations of surface
ations for each grid point from data assimilation, WGHM, and W3 water storage in the same period. The correlation between the
in the study area are displayed in Figs. 6 and 7, respectively. We surfacewaterstorageandsoilmoisturetimeseries(afterremoving
findlargecorrelationsbetweenassimilationresultsandWGHM(0.76 seasonaleffects)isfoundtobe0.92(fora95%confidenceinterval),
onaverageforsoilmoistureand0.82onaverageforgroundwater) 34% higher than the correlation between groundwater and soil
and W3 (0.71 on average for soil moisture) outputs. The results moisture.Thisindicatesthatastrongerconnectionexistsbetween
Fig.5. ComparisonbetweenRMSEachievedfromimplementingeachdataassimilationscenarioaswellasmodel-freerunduring2007.Incase1,surfacestoragesisremoved
fromGRACETWS,incase2,surfacestoragesisaddedtoW3RAsurfacewater,andcase3referstothedataassimilationwithnosurfacestoragecorrection.
Fig.6. Spatialdistributionofaveragesoilwaterstoragevariationsfromdataassimilation,WGHM,andW3.
972 M.Khakietal. /ScienceoftheTotalEnvironment625(2018)963–977
Fig.7. SpatialdistributionofaveragegroundwaterstoragevariationsfromdataassimilationandWGHM.
the surface water storage and soil moisture over the area. Annual each variable are presented in Fig. 9. The spatial distribution of
variationsofgroundwaterstorages,however,showalargerdecline precipitation within Bangladesh indicates larger rainfall in south-
incomparisontosoilmoisturestoragevariations,especiallybetween eastern parts. TWS distribution in EOF2 follows the same pattern.
2008and2012.Asignificantdecreaseingroundwaterstorageisseen LargewaterstoragesarecapturedbyEOF3inthenorthwest.More
inFig.8bwithanaveragerateof8.73±2.45mm/year,showingan informationcanbeextractedfromprecipitationandTWStimeseries.
overall∼46%reduction.Thedeclineinwateravailabilitycanbedue The first three principal components (PC1, PC2, and PC3) of each
toover-extractionofgroundwaterresourcessincesuchadecreaseis datasetareshowninFig.10.Largeprecipitationimpactsarefound
notseeninprecipitation(seeSection4.2formoredetails). in2003and2007.Anegativeanomalyinrainfallisfoundin2010
and 2012, as well as in the period between 2005 and 2007 (PC1).
4.2. Statisticalanalyses TWStimeseriesdemonstratedeclinesbetween2005and2007,and
particularlyafter2009followingadropevent(Mondoletal.,2017).
First, the precipitation and TWS over Bangladesh is analyzed. Results,however,showanincreasein2007inagreementwithENSO
To explore the climate variability and its relationship with water rainfall. The overall average of TWS variations during the study
storages, precipitation will be compared to the data assimilation periodisnegative(∼11.48±3.19mm/year)fortheentirecountry.
results.Principalcomponentanalysis(PCALorenz,1956)isapplied A similar trend, however, is not observed in precipitation even
on GRACE TWS and precipitation time series at each grid point thoughthereisashorterperiodnegativedeclinein2005andafter
to explore their spatio-temporal variations. The first three most 2010. Fig. 10 illustrates that although in some cases a variation in
dominant empirical orthogonal functions (EOF1, EOF2, EOF3) for precipitationresultsinachangesinTWS,continuousTWSreduction
Fig.8. Averagesoilmoisturestorage(a)andgroundwaterstorage(b)timeseriesfromassimilation,WGHM,W3,andin-situmeasurements.