Application and evaluation of a snowmelt runoff model in the Tamor River basin, Eastern Himalaya

17 Pages·2013·2.34 MB·English

HYDROLOGICALPROCESSES
Hydrol.Process.(2013)
PublishedonlineinWileyOnlineLibrary
(wileyonlinelibrary.com)DOI:10.1002/hyp.10005
Application and evaluation of a snowmelt runoff model in the
Tamor River basin, Eastern Himalaya using a Markov Chain
Monte Carlo (MCMC) data assimilation approach
Prajjwal K. Panday,1* Christopher A. Williams,2 Karen E. Frey2 and Molly E. Brown3
1WoodsHoleResearchCenter,149WoodsHoleRoad,Falmouth,MA,02540,USA
2 GraduateSchoolofGeography,ClarkUniversity,950MainStreet,Worcester,MA,01610,USA
3 BiosphericSciences,NASAGoddardSpaceFlightCenter,Greenbelt,MD,20771,USA
Abstract:
Previousstudieshavedrawnattentiontosubstantialhydrologicalchangestakingplaceinmountainouswatershedswherehydrologyis
dominatedbycryosphericprocesses.Modellingisanimportanttoolforunderstandingthesechangesbutisparticularlychallengingin
mountainousterrainowingtoscarcityofgroundobservationsanduncertaintyofmodelparametersacrossspaceandtime.Thisstudy
utilizesaMarkovChainMonteCarlodataassimilationapproachtoexamineandevaluatetheperformanceofaconceptual,degree-day
snowmeltrunoffmodelappliedintheTamorRiverbasinintheeasternNepaleseHimalaya.Thesnowmeltrunoffmodeliscalibrated
usingdailystreamflowfrom2002to2006withfairlyhighaccuracy(averageNash–Sutcliffemetric~0.84,annualvolumebias<3%).
TheMarkovChainMonteCarloapproachconstrainstheparameterstowhichthemodelismostsensitive(e.g.lapserateandrecession
coefficient) and maximizes model fit and performance. Model simulated streamflow using an interpolated precipitation data set
decreases the fractional contribution from rainfall compared with simulations using observed station precipitation. The average
snowmeltcontributiontototalrunoffintheTamorRiverbasinforthe2002–2006periodisestimatedtobe29.7±2.9%(whichincludes
4.2±0.9% from snowfall that promptly melts), whereas 70.3±2.6% is attributed to contributions from rainfall. On average, the
elevation zone inthe 4000–5500m range contributes the mosttobasinrunoff,averaging 56.9±3.6% of all snowmelt input and
28.9±1.1% of all rainfall input to runoff. Model simulated streamflow using an interpolated precipitation dataset decreases the
fractional contribution from rainfall versus snowmelt compared with simulations using observed station precipitation. Model
experimentsindicatethatthehydrographitselfdoesnotconstrainestimatesofsnowmeltversusrainfallcontributionstototaloutflow
butthatthisderivesfromthedegree-daymeltingmodel.Lastly,wedemonstratethatthedataassimilationapproachisusefulfor
quantifyingandreducinguncertaintyrelatedtomodelparametersandthusprovidesuncertaintyboundsonsnowmeltandrainfall
contributionsinsuchmountainouswatersheds.Copyright©2013JohnWiley&Sons,Ltd.
Supportinginformationmay be foundintheonline versionofthis article.
KEYWORDS snowmelt; hydrology;SRM; APHRODITE; MODISsnow; Himalaya;MCMC;data assimilation
Received 11January2013; Accepted 6August 2013
INTRODUCTION tropospherictemperaturesfromthelongestavailablerecord
of microwave satellite measurements reveals widespread
Growingevidenceindicatesthatthelow-latitudemountain-
pre-monsoonal and mean annual warming over the HKH
ous cryosphere is undergoing change owing to recent
region at a rate of 0.21±0.08°C/decade (Gautam et al.,
warming(IPCC,2007;Yaoetal.,2012).Geographicareas
2010). Regional climate projections by the Intergovern-
wherethewatercycleisdominatedbysnowmelthydrology
mental Panel on Climate Change (IPCC, 2007) indicate
areexpectedtobemoresusceptibletoclimatechange,asit
CentralAsiamaywarmbyamediantemperatureof3.7°C
affectstheseasonalityofrunoff(Barnettetal.,2005;Adam
by the end of the 21st century, with warming over high
etal.,2009).TheHinduKush-Himalaya(HKH),regarded
altitudesacrosstheHKH(Pandayetal.,inreview).Thereis
asthe‘watertowers’ofAsia,isonesuchcriticalregionasit
an increasingly large body of evidence of a glacio-
nourishes several of the major Asian watersheds through
hydrological response along the east–west transect of the
perennialriverssuchastheGanges,IndusandBrahmaputra
HKH corresponding to the these climatic changes, with
(Messerli et al., 2004). Most recent research analyzing
glaciers in the Eastern Himalaya exhibiting retreat and
negative mass balance and glaciers in the Karakoram and
northwestern Himalaya exhibiting a positive mass balance
*Correspondenceto:PrajjwalK.Panday,WoodsHoleResearchCenter,
over the last few decades (Bolch et al., 2012; Yao et al.,
149WoodsHoleRoad,Falmouth,MA,02540,USA.
E-mail:[email protected] 2012).Suchclimate-drivenresponsesofmountainousriver
Copyright©2013JohnWiley&Sons,Ltd.
P.K.PANDAYETAL.
hydrology may pose significant challenges for this region. observationstoconstrainmodelswiththegoalofproviding
Despiteitsregionalimportanceandvulnerabilitytoclimate better understanding of the underlying processes
change, there is uncertainty associated with the physical (McLaughlin et al., 2005). The Bayesian approach in
processes that control runoff generation at high altitudes hydrologicalapplicationsallowstosystematicallyincorpo-
(Bolch et al., 2012; Pellicciotti et al., 2012) and subse- rate prior information on model parameters in the data
quently with the rates and magnitude of climate change assimilationframework(Zobitzetal.,2011)andtocompute
impactsonsnowcoverandsnowmelthydrology. probabilitydistributionsofmodelledoutputsandcharacter-
Water resource management and the evaluation of ize uncertaintyinthe model,usuallythroughMonteCarlo
impacts of climatic change require quantification of simulation (Mishra, 2009). This provides a means of
streamflow variability. Hydrological models provide a objectively and quantitatively diagnosing the parameters
framework to investigate these relationships and subse- andprocessestowhichamodelismostsensitiveandfrom
quently use a scenario approach to evaluate projections of whichonecanlearn.Fewstudieshaveemphasizedtheneed
hydrologicalvariables(Maureretal.,2009).Theassessment for parameterization in the Himalayan region and have
of hydrological impacts of climate change is particularly performed extensive calibration (Konz et al., 2007;
challenginginmountainousenvironmentsastheyconstitute Immerzeel et al., 2011; Pellicciotti et al., 2012) and
complex hydrological systems owing to extreme heteroge- validation with additional data sets such as glacier mass
neity in vegetation, soils, topography and spatially and balancesandremotelysensedsnowcover.Pellicciottietal.
temporallyvaryingsnowcoverandsnowmeltpatterns(Gurtz (2012) also emphasized that an assessment of the internal
et al., 2005; Panday et al., 2011). Continued efforts in consistencyofthemodelisimportantowingtotheproblem
hydrologicalmodellingofsnowandglaciermeltintheHKH ofmultiplicityofparametersetsorequifinality.
regionmostlythroughconceptuallylumped,semi-distributed In addition to parameter uncertainty, sensitivity of
(ratherthanphysicallybasedmodels)haveprovidedbaseline hydrological models to input data such as precipitation
patternsofspatiotemporalpatternsofprecipitationandrunoff also needs to be evaluated particularly in complex
at very large scales (Bookhagen and Burbank, 2010; topography and high elevations. Hydrological models in
Immerzeel et al., 2010). However, poor accessibility and mountainousregionswithscarceobservationsareusually
scarcityofadequatenetworkofhydrometeorologicalstations forced with precipitation data located at lower altitudes,
at high altitude regions still remain a major impediment to which can introduce added uncertainty. Satellite-derived
runoff modelling. As a result, only a few studies have precipitation estimates and interpolated rain-gauge obser-
exploredsnowmeltrunoffmodelsatthecatchmentorriver- vations have been used in other studies as inputs to a
basinscaleinthisregion(Braunetal.,1993;Shresthaetal., snowmelt runoff model (SRM) instead of gauged precip-
2004;Akhtaretal.,2008;Tahiretal.,2011). itation. Bookhagen and Burbank (2010) used improved
Parameterizationisusuallychallengingowingtoscarcity Tropical Rainfall Measuring Mission rainfall estimates to
of data in mountainous regions such as the HKH, and a characterize the spatiotemporal distribution of rainfall,
soundunderstandingofthesnowmelthydrologyneedstobe snowfall and evapotranspiration across the Himalayan
emphasized prior to assessment of water resources in a region, whereas Tahir et al. (2011) successfully used
changing climate. Several catchment-scale and basin-scale interpolateddailyprecipitationdatatosimulatesnowmelt
modelling studies have addressed the glacio-hydrological hydrologyintheHunzaRiverbasin.
processes in the Himalayan region (Konz et al., 2007; A thorough assessment of hydrologic processes,
Akhtaretal.,2008;ButtandBilal,2011;Immerzeeletal., patternsandtrendsrequiresminimization ofuncertainties
2011;Tahiretal.,2011).Amajorityofthesestudieshave in modelling and available input data. The objective of
adopted parameter values from the literature without thisstudyistoexaminespatialandtemporalvariationsin
rigorouscalibrationandwithoutquantifyinghowparameter snowmelt contributions to runoff in the Hindu Kush-
uncertaintyinfluencesmodelledresults.Transferofparam- HimalayanregionandTamorRiverbasinbyusingadata
etersacrossspaceandtimeinhydrologicalmodellingcanbe assimilation approach that both combines available input
problematic,particularlyinconceptualmodelsofsnowmelt data with a degree-day based SRM and estimates model
whenparameterseitherdonothaveaphysicalmeaningor parameters by using a Markov Chain Monte Carlo
are not easily measurable (Sun et al., 2012). In order to (MCMC) technique. Additionally, we investigate and
addresstheissuesofparameterizationandmodeluncertain- evaluate the performances of the SRM when driven by
ty, data assimilation methods are typically applied in observed precipitation as well as interpolated gridded
hydrology to integrate different datasets, combine model precipitation data. We demonstrate how a data assimila-
outputs with observations, update hydrological models, tion approach can be coupled with a SRM to investigate
quantifymodeluncertaintiesandsensitivityand/orimprove the issues of model parameter sensitivity and uncertainty
modelaccuracy(BevenandFreer,2001;Clarketal.,2006). and improve performance in a data-scarce basin in the
Data assimilation can be viewed as a method to use Eastern Himalaya.
Copyright©2013JohnWiley&Sons,Ltd. Hydrol.Process.(2013)
APPLICATIONANDEVALUATIONOFASNOWMELTRUNOFFMODELINTHEHIMALAYA
METHODS
Study area
This geographic domain of this study is the upstream
portionoftheTamorRiverbasinintheEasternHimalaya
of Nepal (Figure 1). The Tamor basin covers an area of
4274km2, with an elevation range of 422 to 8505m
(Figure 2). The Tamor River, which is one of the major
tributaries of the Sapta Koshi River, flows between high
peaks such as Kanchenjunga (8586m) and Everest
(8848m) on gorges reaching depths of 2000m. Precip-
itation for the region is highly seasonal and concentrated
duringthesummermonsoonmonths(JunetoSeptember)
(Figure 3). The upper reaches of the basin are snow
covered and glaciated, contributing melt water to
streamflow.
Figure2. Area-elevationcurvefortheTamorRiverbasin
Hydrometeorological data
acquired from the Consultative Group on International
Hydrometeorological measurements for the Tamor
AgriculturalResearch–ConsortiumforSpatialInformation.
River basin in the Eastern Himalaya were acquired from
Thebasinwasdividedintofourzones,andthezonalmean
the Department of Hydrology and Meteorology, Nepal.
hypsometricelevation(Figure2)wasusedastheelevation
Daily temperature data at Taplejung station (27.35°N,
towhichbasestationtemperatureswereextrapolatedforthe
87.67°E at 1732m elevation) and daily precipitation data
calculationofzonaldegreedays(TableI).
measured at Lungthung station (27.55°N, 87.78°E at
1780m elevation) in the Taplejung district were used for
ModerateResolutionImagingSpectroradiometer(MODIS)
the hydrological modelling (Figure 1). Discharge data
measured daily from 2002–2006 for the Tamor River at snow cover
Majhitarstation(27.15°N,87.71°Eat533melevation)were Satellite-derived snow cover is the most efficient and
usedforcalibration.TheTamorbasinwasdelineatedatthe readily available input for snowmelt runoff modelling.
Majhitar station using the 90m digital elevation data Snow cover was derived using the MODIS/Terra Snow
Figure1. LocationandtopographyoftheTamorRiverbasinintheeasternNepaleseHimalaya
Copyright©2013JohnWiley&Sons,Ltd. Hydrol.Process.(2013)
P.K.PANDAYETAL.
AsianPrecipitation–HighlyResolvedObservational
DataIntegrationtowardsEvaluation(APHRODITE)
precipitationdata
Spatiallyrepresentativegriddedprecipitationdataarealso
used to drive a snowmelt runoff model, as valley stations
may not be representative for basin precipitation. The
APHRODITEofthewaterresourcesprojecthascompileda
daily gridded precipitation dataset covering 1951–2007
periodbyanalyzingrain-gaugeobservationdataacrossAsia
(Yatagaietal.,2009).TheAPHRODITEdatausedinthis
research (APHRO_MA_V1003R1) have a 0.25×0.25°
gridded spatial resolution, which was interpolated from
meteorological stations throughout this region. Whereas
satellitedataconsistentlyunderestimateoroverestimatethe
amountofprecipitation,particularlyinthemonsoonseason
Figure 3. Average monthly streamflow (m3/s) and total monthly (Duncan and Biggs, 2012), APHRODITE data have been
precipitation(mm)from1996to2006atMajhitarstation
shown to provide the best precipitation estimates in most
regionsexceptathighelevationswhereaccuracyofthisdata
setislimited(Andermannetal.,2011).Weusedthebasin-
TableI. Elevationzone,elevation range,mean elevationand
averaged APHRODITE precipitation as input to the SRM
zonalarea forthe four zonesof theTamor River basin
assumingittobeconstantacrosselevationzones.Themodel
Elevation Meanelevation Zonal wasrunusingAPHRODITEprecipitationdatatodetermine
Elevationzone range(m) ofzone(m) area(km2) thereliabilityofthedatasetforsnowmeltrunoffmodelling
intheTamorRiverbasin.
1 533–2500m 1685.9 1355.5
2 2500–4000m 3208.1 1111.3
3 4000–5500m 4764.3 1297.9 Snowmelt runoff modelling
4 >5500 5988.8 509.2
Modelling of stream discharge was based on the SRM
Totalbasin — — 4273.8
developed by Martinec et al. (2008). SRM is a conceptual,
deterministic, degree-day hydrological model that simulates
and forecasts daily runoff resulting from snowmelt and
Cover8-DayL3Global500mGrid(MOD10A2;Halletal.,
precipitationfromhydrometeorologicalinputdata.Themodel
2006)andupdatedeveryeightdaysusingasnowdepletion
was designed for streamflow modelling in basins where
curve.TheMODISsnowcoverisbasedonasnowmapping
algorithmthatusesat-satellitereflectanceinMODISbands snowmelt is a major contributor, and more recently, it has
4(0.545–0.565μm)and6(1.628–1.652μm)tocomputethe beenappliedforevaluationofimpactsofclimatechangeon
seasonal snow cover and runoff (Immerzeel et al., 2010).
Normalized Snow Difference Index (Hall et al., 2002). A
SRM has been applied to basins ranging from 0.76 to
pixel will be mapped as snow if the Normalized Snow
Difference Index≥0.4 and MODIS band 2 (0.841– 917444km2andacrosselevationrangesfrom0to8840min
0.876μm) >11%; however, if the MODIS band 4 over 100 basins in 29 different countries (Martinec et al.,
reflectanceis<10%,thepixelwillnotbemappedassnow 2008). This model was chosen in this study to assess
snowmelthydrologyowingtoitsminimaldatarequirements
even if the other criteria are met. The MODIS product
andtransferabilitytoavarietyofgeographicregions.Onthe
distributedthroughtheNationalSnowandIceDataCentre
basisoftheinputofdailytemperature,precipitationandsnow-
was acquired for the scene h25v06, covering the study
coveredareavalues,SRMcomputesthedailystreamflowas
region from 2002 to 2006. The MODIS images were
reprojected to Universal Transverse Mercator zone 45
projectionsystem,andtheTamorRiverbasinareawasthen Qnþ1 ¼∑½CsznaznðTnþΔTznÞSzn (1)
extracted toassessthe fraction ofsnowcover for different
altitudinal zones. Average snow cover was estimated by þCrzn PznÞAzð10000=86400Þð1(cid:2)knþ1Þ
interpolating linearly between the previous and next
available cloud free images for days with cloud cover þQn knþ1;
following Tahir et al. (2011). The conventional depletion knþ1 ¼xQyn; (2)
curveswerederivedforeachelevationzoneusingthe8-day
snowcoverdata,whichprovidetheinputtotheSRMasthe where Q (m3/s) is discharge at day n+1, C is the runoff
decimalportionofsnow-coveredareaforeachday. coefficienttosnowmelt(Cs )andtorain(Cr )foreachzone
n n
Copyright©2013JohnWiley&Sons,Ltd. Hydrol.Process.(2013)
APPLICATIONANDEVALUATIONOFASNOWMELTRUNOFFMODELINTHEHIMALAYA
(subscript z), a (cm/°C/day) is the degree-day factor, toaccountingforlossesviarunoffcoefficients,introducing
T+ΔT(°C)arethedegreedays,P(cm)istheprecipitation animplicitassumptionthatfractionalcontributionstorunoff
perday,A(km2)istheareaofthebasin,Sisthefractional areproportionaltotherelativeabundanceofeachinput.
snowcoverandkisthedischargerecessioncoefficient.The Equation(1)isappliedtoeachzoneofthebasin,andtotal
daily average streamflow on day n+1 is computed by discharge iscomputedas sumofallzonal discharges.The
addition of snowmelt and precipitation that contributes to SRM structure as applied inthis study assumes a timelag
runoff with the discharge on the preceding day. The between the daily temperature cycle and the resulting
snowmelt of the preceding day is a product of the degree- dischargecycleof18h.Theuseofthistimelagissupported
day factor (a), zonal degree days (T+ΔT) and the snow- bytheexpectedrelationbetweenlagtimeandbasinsizeas
coveredareapercentage(S).Thesnowmeltrunoffcoefficient shown in a SRM intercomparison report by the World
(Cs) and the zonal area (A) are further multiplied with the MeteorologicalOrganization(1986).Althoughthetimelag
previousproducttocomputethepercentagethatcontributes parameter can be adjusted to improve synchronization of
to runoff. Similarly, the product of the rainfall runoff simulated and observed peaks of average daily flows,
coefficient (Cr) and the zonal area (A) determines the calibration of the recession coefficients provides a similar
precipitationcontributiontorunoff.Therecessioncoefficient effect (Martinec et al., 2008). The parameters for runoff
(k)isalsoanimportantparameterasitdescribesthestructure coefficients, lapse rate, critical temperature, recession
of the falling limb of the hydrograph and correspondingly coefficient and degree-day factors reflect characteristics at
determinestheproportionofdailymeltwaterproductionthat abasinscaleandaredeterminedinthisstudywithanMCMC
appears immediately in the runoff. This is obtained with dataassimilationbymaximizingthefitbetweentheobserved
analysis of historical discharge data based on Equation 2, and modelled streamflows. The modelled streamflow is
whereparametersxandyarederivedbyanalyzingdischarge assessedforaccuracyusingtheNash–Sutcliffe(NS)statistic
onagivendayasafunctionofthevalueonthepreviousday. and the volume difference (D ), where volume difference
v
A critical temperature is used to determine whether a demonstrates bias. (Supplementary information, A)
precipitationeventisrainorsnowoverdifferentpartsofthe (Martinecetal.,2008).
basin. When precipitation is determined to be snow, its
contributiontorunoffdependsonwhetheritfallsonasnow- Model calibration using MCMC optimization
covered or previously snow-free area. Snowfall on snow- AdataassimilationtechniqueknownasaMCMCroutine
covered areas is assumed to become part of the seasonal (Zobitz et al., 2011) was used to explore the parameter
snowpack, its depth is not tracked further, and its spacesintheSRM.TheapplicabilityofBayesianmethods
contributiontorunoffisdeterminedbythesnowdepletion inhydrological studieshasincreased owing totheiruse in
curve. Daily snowmelt depth from snow-covered areas is parameter estimations and uncertainty assessments (Smith
determined on the basis of degree days and a factor that and Marshall, 2008). The MCMC algorithm is one of
convertsthistomassrelease,qualitativelyrepresentingthe severaldataassimilationtechniquesandusestheMetropolis
maximumpotentialmeltgiventheenergyinput.Incontrast, Hastings algorithm (Metropolis et al., 1953), which
snow falling as precipitation over snow-free areas is iteratively adjusts the model parameters to yield the best
accumulated until the next day warm enough to produce matchbetweentheobservedandmodelledstreamflow.The
melting at which time it is depleted on the basis of the output provides a set of accepted parameter values and
degree-days approach. During the winter half year information about the posterior distributions of the
(October–March), the model cancels any existing storage parameters consistent with the prior information, model
of preceding temporal snowfalls when snow cover is structureandobservations.AdditionaldetailsoftheMCMC
maximum such that snow on previously snow-free areas algorithmusedforoptimizationaredescribedmorefullyin
becomes part of the seasonal cover. When precipitation is Supplementary information, B and in Zobitz et al. (2011)
consideredtoberainduringthecoldseason,rainfallrunoff andBraswelletal.(2005).
contributestototalrunoffonlyfromthesnow-freearea,and Basedontheliteratureandrangeofparametervaluesused
therainfalldepthisreducedbytheratioofsnow-freeareato in other snowmelt runoff models in the HKH region
zonalarea.Atalaterstageduringtheseason,whichissetto (Immerzeeletal.,2010; ButtandBilal,2011; Tahiret al.,
April 1, rain falling on snow-covered areas is allowed to 2011), we first specified prior ranges for the parameters.
fully contribute to runoff such that rain from entire zonal TableIIshowstheparameterswiththepriorrangesthatwere
area is accounted. The SRM approach estimates fractional assigned seasonally or annually across different elevation
contributions to total discharge from (1) daily snowmelt zones. Seasonally varying parameters were optimized for
from the snowpack for each elevation, (2) the melt of the snowmelt and monsoon period (June–August) and the
precipitation delivered as snow but to areas described as snow accumulation and dry period (September–May).
snow-free at the time of snowfall and (3) zonal rainfall. Snowmelt and rainfall runoff coefficients were varied
Contributionsofeachtototalriverflowareestimatedprior seasonallyandbyelevationsinformedbyTahiretal.(2011)
Copyright©2013JohnWiley&Sons,Ltd. Hydrol.Process.(2013)
P.K.PANDAYETAL.
TableII. Initial parameter rangesusedfor MarkovChainMonte Carlodata assimilation
Parameterranges across eachelevation zone
Parameter
Zone1 Zone2 Zone 3 Zone 4
(533–2500m) (2500–4000m) (4000–5500m) (>5500m)
Snowrunoff coefficient June–August 0.3–0.8 0.3–0.8 0.3–0.8 0.3–0.8
September–May 0.3–0.8 0.3–0.8 0.3–0.8 0.3–0.8
Rainfall runoffcoefficient June–August 0.3–0.8 0.3–0.8 0.3–0.8 0.3–0.8
Sep–May 0.3–0.8 0.3–0.8 0.3–0.8 0.3–0.8
Degree-day factor(cm/°C/day) Annual 0.3–0.6 0.3–0.6 0.5–0.9 0.5–0.9
Lapserate (°C/100m) Annual 0.5–0.7
Critical temperature(°C) Annual 0.0–2.0
xcoefficient Annual 0.9–1.4
ycoefficient Annual 0.0–0.25
SnowmeltandmonsoonalmonthsareJune–August,andSeptember–Mayisthesnowaccumulationanddryperiod.Lapserate,criticaltemperatureandx
andycoefficientsarenotvariedacrosselevationzones.
andwereassignedequalpriorparameterranges.Although Model sensitivity and uncertainty
degree-day factors tend to vary with different seasons,
In addition to model calibration, we also carried out
climate and surfaces, under similar conditions they are several experimental runs. The model was allowed to
expectedtoincreasewithincreasingelevation(Hock,2003). calibrate via MCMC by setting the snowmelt runoff
Similar to previous studies (Butt and Bilal, 2011; Tahir coefficients to zero such that snowmelt contributions were
etal.,2011),whichutilizeadegree-dayfactorthatincreases neglectedinthesimulationofthehydrograph.Similarly,the
with elevation, we also increase the parameter range rainfall runoff coefficients were also set to zero such that
graduallyacrosselevationzones.Theelevationzonegreater rainfallcontributionswereexcludedinthefinalsimulationof
than5500m,whichareglacierandicecoveredareassumed the hydrograph. In addition to these experiments, we
tohaveahigherdegree-dayfactorasthevaluestendtobe conducted additional sensitivity and uncertainty analyses
considerably higher for ice compared with snow (Hock, viaensemblestreamflowsimulations(1000runs)toidentify
2003). Parameters such as lapse rate, critical temperature parametersthatcontrolmodelperformanceanduncertainties
and x and y coefficients were calibrated annually. Other in modelled streamflow and total input contributions.
studies using SRM have used a similar variation of Parameter sensitivity and uncertainty in streamflow in the
parameters.Althoughmostparameterscanpotentiallyvary time history were analyzed as 5% and 95% confidence
across elevations and seasons, increasing parameter sets intervals of the ensemble runs. First, we generated an
increase potential ranges of parameter values and thereby ensembleofmodelruns(‘bestensemble’hereafter)fromthe
increase uncertainty (Seibert, 1997). For this study, we posteriordistributionofacceptedparametersetsfromMCMC
implementedtheMCMCtechniqueusingthreechains,each that represent the most likely set of parameters. Second,the
ofwhichconsistedof10000iterations.Thebestparameter sensitivity to parameters was tested by perturbing the
valuesfromthesechainswereusedinthefinalsetof10000 distribution of a parameter by some amount and computing
iterations, which is used to determine the parameter thecorrespondingchangesinmodeloutputwhilepreserving
distributions. The range of parameter values was deliber- thedistributionofotherparameters.Thisisperformedinthis
atelysettobelargeatfirsttoexamineifrealisticparameter study by randomly drawing from the altered distribution of
values were selected by the MCMC technique. If the each parameter (e.g. ±10% from the mean value of a
posterior parameter distributions were less plausible, the parameter from the posterior distribution) to generate
parameter ranges were logically constrained by literature ensemble runs. This allows quantifying the effects of each
valuestobemorerealistic.Thesnowmeltandrainfallrunoff parameter or a combination of parameters on model output
coefficients were constrained in the ranges specified in and on total input contributions to streamflow by providing
Table II based on Tahir et al. (2011) and Butt and Bilal uncertaintyboundsontheoutputstatistics.Third,anensemble
(2011),degree-dayfactorrangeswerespecifiedonthebasis ofparametersetswasgeneratedbyrandomlysamplingfrom
ofHock(2003),andpriorrangesforlapserateandcritical the cumulative distribution curve constructed from the
temperaturewerespecifiedfromTahiretal.(2011)andButt posteriordistributionofeachparameter.Theensemblemodel
andBilal(2011).TheSRMwasrunviaMCMCwithprior simulation uses the full range of parameter variations and
parameterrangesasspecifiedinTableIIforthehydrological combinationsandtherebyprovidesoveralluncertaintiesinthe
years2002to2006. choiceofmodelparametersgeneratedfromMCMC.
Copyright©2013JohnWiley&Sons,Ltd. Hydrol.Process.(2013)
APPLICATIONANDEVALUATIONOFASNOWMELTRUNOFFMODELINTHEHIMALAYA
RESULTS maximum snow cover of ~40% with summer zonal snow
coveragefallingto~12%.Elevationsupto4000mhadon
Summary of meteorological parameters
average~9%annualmaximumsnowcover.
Long term annual temperatures (1970–2006) at the
Taplejungstationaveraged16.3°C,whereasaveragewinter
Model parameterization
and summer temperatures were 10°C and 21.3°C,
respectively. Mann–Kendall tests indicate significant Whenthepriorrangesoftheparametersweredeliberately
(p<0.05) positive trends in both annual maximum and set to be large, the model calibrated well with a good fit
average temperatures for the Taplejung station during this between the observed and modelled streamflow; however,
period. Average annual precipitation at the Lungthung thisalsoresultedinseveralmodelswithsimilarperformance
station was ~2484mm for the 1996–2006 period. On and physically implausible coefficients. Logically cons-
average,approximately76%oftheprecipitationisreceived training the prior parameter ranges provided models with
during the summermonsoon.Themonthly streamflowfor betterfitandparametersthatwerephysicallyplausiblewhen
theTamorRiverbasinattheMajhitarstationalongwiththe varied across seasons and elevations. The posterior
totalmonthlyprecipitationatLungthungstationindicatesa distributions of parameters after the MCMC optimization
fingerprint of monsoonal precipitation in the Eastern forallparametersforthe2002–2003hydrologicalyearsare
Himalaya from June to September (Figure 3). Average showninFigure5.Whenthemodelwasoptimizedusingthe
streamflow measured at the Majhitar station averaged MCMC approach, some of the most sensitive parameters
256.2m3/s annually during the same period, whereas the (suchasxandycoefficients,criticaltemperatureandlapse
monsoon months had an average streamflow of 573m3/s. rate) exhibited a narrow, unimodal distribution with their
Snow cover depletion curves for the Tamor River basin ranges significantly reduced from the initial parameter
indicatethatthehighestelevationzone(above5500m)had ranges. The x and y coefficients required to compute the
a maximum ~85% snow-covered area for the 2002–2006 recessioncoefficientwerethemostsensitiveparameters,as
period, dropping to ~73% during summer (Figure 4). The they determinedthe structureof hydrograph recessionand
elevation zone between 4000 and 5550m had an annual correspondingproportionofdailymeltwaterthatappeared
(a) (b)
(c)
Figure4. SnowcoverdepletioncurvesfromtheModerateResolutionImagingSpectroradiometer8-daymaximumsnowcover500-mresolutionproductinthe
TamorRiverbasinfrom2002to2006in(a)zone2,elevationrangefrom2500–4000m,(b)zone3,elevationrangefrom4000–5500mand(c)zone4,elevation
rangegreaterthan5500m.
Copyright©2013JohnWiley&Sons,Ltd. Hydrol.Process.(2013)
P.K.PANDAYETAL.
Figure5. Posteriorparameterdistributionsforthe2002–2003hydrologicalyearsashistograms
immediately in the runoff. This shows that the MCMC that contributes to runoff and are expected not to vary as
approach performs well in constraining the most sensitive much across types of precipitation. Some of the observed
parameters of the model. Once a narrow posterior differencesinthesecoefficientswithinasinglezonecould
distribution is established for the sensitive parameters bearesultofthetrade-offamongparametersafterthemost
during the MCMC approach, changing other parameters sensitiveparametershavebeenconstrainedwell.Inadequate
does not result in a significant increase in model precipitation data from higher elevations zones in the
performance.Therefore,someoftheremainingparameters mountainous regions further hindered our ability to
suchassnowmeltandrainfallrunoffcoefficientsexhibited constrain the runoff coefficients. Parameter estimates tend
broadposteriordistributions,andlargestandarddeviations to be stable across years, though Tcrit showed large
where good simulations of the observed streamflow were variation, but also had a large spread. Degree-day factor
obtained over broad ranges (Figure 5 and Table III). The estimatestendedtobethehighestforelevationsgreaterthan
posteriordistributionsofrunoffcoefficientsforbothrainfall 4000m (zones 3 and4), which had the highestfraction of
andsnowmeltarecomparableacrosstypesofprecipitation snow-coveredarea.Thepeaksinposteriordistributionsfor
and all simulated years. These coefficients represent the the degree-day factor correspondingly rise from zone 1 to
fraction of liquid water, from either snowmelt or rainfall, zone3elevationranges.
Copyright©2013JohnWiley&Sons,Ltd. Hydrol.Process.(2013)
APPLICATIONANDEVALUATIONOFASNOWMELTRUNOFFMODELINTHEHIMALAYA
Table III. Meanandstandard deviation ofparameters from theposterior distribution
Meanandstandard deviation ofparameteracross zones
Parameter
Zone 1 Zone2 Zone3 Zone 4
Year (533–2500m) (2500–4000m) (4000–5500m) (>5500m)
Cs June–August 2002 0.55±0.15 0.59±0.13 0.72±0.07 0.56±0.13
2003 0.46±0.13 0.44±0.11 0.47±0.12 0.52±0.15
2004 0.54±0.14 0.54±0.14 0.68±0.11 0.55±0.14
2005 0.55±0.14 0.55±0.15 0.69±0.11 0.55±0.15
September–May 2002 0.58±0.14 0.61±0.11 0.59±0.13 0.54±0.15
2003 0.60±0.14 0.56±0.13 0.71±0.07 0.56±0.15
2004 0.70±0.09 0.71±0.07 0.67±0.11 0.54±0.14
2005 0.59±0.14 0.43±0.11 0.40±0.0/8 0.54±0.14
Cr June–August 2002 0.73±0.05 0.74±0.05 0.69±0.09 0.55±0.15
2003 0.67±0.11 0.56±0.13 0.54±0.15 0.54±0.14
2004 0.71±0.05 0.68±0.09 0.59±0.14 0.55±0.14
2005 0.76±0.04 0.74±0.08 0.76±0.05 0.56±0.14
September–May 2002 0.69±0.07 0.68±0.09 0.65±0.10 0.55±0.14
2003 0.37±0.07 0.43±0.09 0.40±0.08 0.54±0.14
2004 0.71±0.07 0.70±0.08 0.54±0.14 0.56±0.14
2005 0.69±0.07 0.70±0.08 0.72±0.07 0.51±0.13
a Annual 2002 0.48±0.08 0.50±0.0 0.61±0.08 0.68±0.102
2003 0.42±0.08 0.41±0.08 0.79±0.0 0.69±0.11
2004 0.53±0.05 0.53±0.06 0.81±0.05 0.70±0.`2
2005 0.45±0.09 0.44±0.09 0.80±0.08 0.71±0.12
LR Annual 2002 0.64±0.01
2003 0.51±0.01
2004 0.68±0.01
2005 0.54±0.01
Tcrit Annual 2002 1.26±0.28
2003 1.16±0.49
2004 1.83±0.13
2005 0.94±0.62
x coefficient Annual 2002 1.34±0.03
2003 1.26±0.01
2004 1.29±0.02
2005 1.11±0.01
y coefficient Annual 2002 0.08±0.01
2003 0.06±0.01
2004 0.07±0.01
2005 0.03±0.01
Cs,snowmeltrunoffcoefficient;Cr,rainfallrunoffcoeficient;a,degreedayfactor;LR,lapserate.
The SRM is usually calibrated manually where known months. Across all simulated years, the observed and
ranges ofparameter values aremaintained while changing modelled streamflow show a linear association with some
otherparametersbytrialanderror.Thistypeofcalibration biasduringhighdischargeperiods[Figure6(b,d,fandh)].
can be inefficient in exploring the parameter space and ThemodelsimulatedstreamflowshaveahigherNSstatistic
resulted in an average NS statistic of 0.64 and volume in the years where the summer peak discharge is lower
differenceof10%acrossallsimulatedyearsforthisstudy. relativetootheryears(TableIVandFigure6).
ModelcalibrationusingMCMCoptimizationhadagreater
overall improvement with NS statistic of ~0.84 averaged Model sensitivity and uncertainty
acrossallsimulatedyears.Figure6(a,c,eandg)showsthe Theresultsoftheexperimentrunsindicatethatwiththe
modelledstreamflowagainstobservedstreamflowusingthe snowmelt runoff coefficients (Cs) set to zero, the model
MCMC approach for four hydrological years (utilizing optimized as well as when Cs values were also included
observedprecipitationattheLungthungstation).TheSRM through the MCMC technique. The rainfall runoff
is able to capture the daily streamflow in the annual coefficients (Cr) are parameterized to higher values, and
hydrologyofthebasin;however,themodelunderestimates themodelisabletosimulatethehydrographbyoptimizingx
high discharge events, particularly during the monsoonal andycoefficientswithmodelperformanceequivalenttothat
Copyright©2013JohnWiley&Sons,Ltd. Hydrol.Process.(2013)
P.K.PANDAYETAL.
Figure6. TimeseriesandscatterplotsofobservedandmodelledstreamflowusingobservedprecipitationatLungthungstationfor2002–2003(a,b),
2003–2004(c,d),2004–2005(e,f)and2005–2006(g,h)
withthebestestimatesofparametersasinTableIV.Thefirst accepted by the MCMC approach. Because the most
ensembleparametersets(bestensemble)randomlyselected sensitive parameters have narrow, unimodal distributions,
fromtheposteriordistributionoutputfromMCMCdidnot similarmodelfitandperformanceisachievedwithinabroad
exhibit any significant variations in model fit and perfor- distribution of other parameters. However, when the
mance. This is to be expected as the parameter sets are distribution of the most sensitive parameters was altered
Copyright©2013JohnWiley&Sons,Ltd. Hydrol.Process.(2013)