The footprint of Asian monsoon dynamics in the mass and energy balance of a Tibetan glacier

Introduction Conclusions References


Introduction
Changes in the Asian monsoon climate and associated glacier responses have become predominant topics of climate research.Their importance is founded on, e.g.(i) the effect of monsoon activity on the livelihood of millions of people (e.g.Tao et al., 2004), (ii) the effect of glacier change on regional water supply and global sea level rise (Kaser et al., 2006(Kaser et al., , 2010)), and (iii) the role of the monsoon in gobal teleconnections, which was discussed by Webster et al. (1998) and continues to be a research focus (e.g.Li et al., 2010;Park et al., 2010).
Understanding the direct influence of atmospheric conditions on glacier mass requires quantitative knowledge of the surface energy balance (SEB) and its relation to the mass balance (MB) of a glacier.Compared to mid and high latitudes, there are few detailed SEB/MB studies for the high Asian mountains.This scarcity was emphasized more than a decade ago by Fujita and Ageta (2000) and it has only slightly improved due to the difficulty of field data collection.A number of such studies are available for the Himalaya at the border region of Nepal and China (Kayastha et al., 1999;Aizen et al., 2002), including glacier-wide analyses.Other examined regions are the interior of the Tibetan Plateau in the Tanggula Mountains (Fujita and Ageta, 2000) and the plateau's north margin (Jiang et al., 2010), where distributed quantifications have been conducted, and the maritime southeast of the Tibetan Plateau, where point SEB studies were performed (Xie, 1994;Yang et al., 2011).
In the Nyainqêntanglha Mountains on the south-central plateau section, Caidong and Sorteberg (2010) modeled the SEB/MB of Xibu Glacier.However, the authors themselves noted that their study only revealed basic features, since no direct measurements from the glacier were available.For Zhadang Glacier in the same mountain chain, Zhou et al. (2010) investigated the runoff variability in relation to local air temperature and precipitation and made conceptual SEB considerations.In this paper we calculate the glacierwide SEB/MB of Zhadang Glacier in detail with the help of in-situ measurements.We therefore extend the knowledge of SEB/MB processes on the Tibetan Plateau, which meets the need to better understand the heterogeneous glacier changes on the plateau (e.g.Yao et al., 2012).However, we also pursue a specific application: to investigate the SEB/MB link to Asian monsoon dynamics.The study goals are therefore (i) to reveal the SEB/MB terms and their space-time variability; and (ii) to analyze these terms in relation to specific periods that reflect both intra-and interannual monsoon variability (the "monsoon footprint" on the glacier), since glacier fluctuations in the high Asian mountains have traditionally been linked to Asian summer monsoon conditions (e.g.Fujita and Ageta, 2000;Yang et al., 2011).
Our study is unique due to the combination of a multiyear data set (vs. single-season investigations), the use of distributed analysis (vs. point studies), and, most especially, the focus on the linkage of SEB/MB to multi-temporal monsoon dynamics (vs.focus on the monsoon onset, i.e. on interannual variability alone).Despite the heterogeneous glacier changes on the Tibetan Plateau, which result from interacting large-scale atmospheric flows and local relief factors (He et al., 2003;Fujita and Nuimura, 2011;Scherler et al., 2011), a common feature is the strong sensitivity to climatic conditions in summer (wet season).This was found by all previous SEB studies and stems from the regional climatic conditions, which show peaks in the annual air temperature and precipitation cycle at the same time.Of particular importance is the precipitation phase in response to air temperature variations (e.g.Kayastha et al., 1999;Fujita and Ageta, 2000;Caidong and Sorteberg, 2010), which has a double effect on MB through changes in accumulation by solid precipitation, and through changes in ablation by albedo.Thus, our goals build on the working hypothesis that local conditions in the monsoon season (boreal summer) govern the mass variability of Zhadang Glacier, which was also postulated from runoff observations downstream of the glacier (Kang et al., 2009;Zhou et al., 2010).

Methods, data basis and climatic setting
A physically-based SEB/MB model (Sect.2.3) is the main tool employed in this study.The forcing data are mainly provided by on-site measurements (Sect.2.1) and also partly by atmospheric model output (Sect.2.2), while measurements not used for forcing serve to evaluate the MB model's performance.We rely on large-scale reanalysis and satellite data Fig. 1.Zhadang Glacier (grey shading) in the digital terrain model used for mass balance modeling.The locations of available measurements are indicated, and contours are in meters a.s.l.(100 m spacing).Note that the glacier is debris-free.The inlay shows the atmospheric model domain and the two nested domains as smaller rectangles (Maussion et al., 2011).Horizontal resolution increases towards the innermost domain that contains the Zhadang area: 30, 10, and 2 km grid spacing.

In-situ measurements
The pivotal observational data are the records from automatic weather station (AWS) 1 on Zhadang Glacier, which is situated at 5665 m in the ablation zone (Fig. 1).Characteristics and usage of available data are summarized in Table 1.AWS1 has been in operation since 2009, but unfortunately has two data gaps due to the extreme environment.However, three periods of sufficient data coverage exist for our study: 27 April-14 July 2009 (period 1), 1 October 2009-25 June 2010 (period 2), and 16 August 2010-15 September 2011 (period 3).This yields a total of 743 days, which is the time frame for the MB modeling and the analysis of the monsoon impact.All data were post-processed according to the instrument manuals, but also carefully screened following the procedure described in Mölg et al. (2009a) in order to detect and correct measurement errors.Basically, this procedure aims to identify periods when radiative sensor heating and/or riming of instruments may have occurred.Only 0.9 % of wind data and between 0 and 0.7 % of the other variables needed correction, which indicates high data quality.
Note that we use accumulation recorded by the sonic ranger (SR50) only for the months of October to April (Table 1), since air temperature (T a ) at AWS1 in the summer months clearly exceeds the melting point and therefore the SR50 does not capture the liquid fraction of precipitation.However, total precipitation is a required input for the model (Sect.2.3).Further, subsurface temperature measurements at AWS1 were also performed at depths less than 5 m, but were obviously affected by radiative heating.Thus, we only use data from the initial depths of ≈ 5.6 and ≈ 9.6 m for MB model evaluation and lower boundary condition definition, respectively (Table 1).Regarding the latter, data from ≈ 9.6 m depth show an almost constant temperature of 268.6 K (standard deviation 0.06 K).We initialize the MB model after midnight on the first day of each of the three periods, and include measured surface temperature as well as data from the sensors at ≈ 3 and ≈ 1.5 m depth since the heating problem during night is minimal.
The time-dependent height of the T a /relative humidity (RH) and wind speed sensors, as well as the depth of subsurface sensors, can be estimated from the SR50 record.Within period 3, however, surface height change is not available from 4 October 2010 to 30 June 2011 (68 % of period 3).This gap was filled with SR50 data from nearby AWS2 (Fig. 1), which has been operated since 2 October 2010 at 5566 m altitude (wind instrument from 21 May 2011 onward).AWS2 does not record ice ablation, as it is located slightly off the glacier margin on a moraine.Nevertheless, the data gap coincides with a period of weak ablation, which will be shown in the results section.Thus, we deem this correction to be a reasonable alternative to assuming a random sensor height.Wind speed and T a /RH from AWS2 (same instruments as at AWS1) were also used to examine vertical gradients, but only for intervals when AWS1 and AWS2 sensor heights agreed within 0.2 m and when there was high enough wind speed for sufficient natural ventilation of the T a /RH sensor (Georges and Kaser, 2002; n = 1629 h for wind; n = 1340 h for T a /RH).We find no gradients in wind speed and RH, as differences are within the instrument accuracy.Typical T a gradients were considered for defining MB model parameters (Sect.2.3).
Finally, data from a precipitation gauge (Geonor T-200B; sensitivity ≤ 0.1 mm) at AWS2 between 22 May 2010 and 15 September 2011 are used to (i) constrain atmospheric model output (Sect.2.2) and (ii) obtain winter precipitation for period 3 when SR50 data at AWS1 are unavailable.Also, several ablation stakes on Zhadang Glacier serve in the MB model evaluation.The locations of these measurement sites are presented in Fig. 1.Stake readings performed by personnel from the Institute of Tibetan Plateau Research are available for three intervals that overlap with each of the three model simulation periods once, as will be detailed in Sect.3.1.

Atmospheric modeling
We run the Advanced Weather Research and Forecasting (WRF) numerical atmospheric model (Skamarock and Klemp, 2008) with a domain that covers a large part of Asia and the Northern Indian Ocean (Fig. 1) for the years 2009-2011.Multiple grid nesting in the parent domain yields a local-scale spatial resolution over Zhadang Glacier of 2 km (Fig. 1), which reproduces the real terrain altitude well (see below).All details of the model configuration are presented in Maussion et al. (2011), e.g. the grid structure (their Fig. 1) and the chosen model options and forcing strategy (their Table 1).We initialize WRF every day from the Global Forecast System's final analysis (Maussion et al., 2011), which differs from the traditional regional climate model approach where the model fields evolve over the entire simulation time.WRF output here, by contrast, is strongly constrained by the observed synoptic weather patterns and thus represents a regional atmospheric reanalysis.
Running WRF and producing a high-resolution reanalysis data set are part of a larger project, but for the present study we are only interested in the simulation of two variables that www.the-cryosphere.net/6/1445/2012/The Cryosphere, 6, 1445-1461, 2012 are not available from measurements: local precipitation during the monsoon season and cloud cover fraction at Zhadang Glacier, both from the 2 km-resolution grid.The cloud cover is determined for a 50 km view field around AWS1 following the method of Mölg and Kaser (2011), which diagnoses saturated regions using the total condensate mixing ratio.Simulated precipitation is extracted from one of the four grid points surrounding AWS1 that has the least altitude difference (only 8 m).However, we constrain the WRF output with the ≈ 16 months of measurements from the precipitation gauge (Fig. 2a).Simulated precipitation shows an excellent correlation with the measured seasonal evolution (r = 0.99; Fig. 2a), a point that is vital for this study.Also, the simulated diurnal cycle in the summer monsoon months with the main and secondary maximum during night and early afternoon (Fig. 2b), respectively, agrees with observations on the Tibetan Plateau (Ueno et al., 2001).Thus, we have confidence in the simulated variability, but apply a scaling factor (K WRF ) of 0.56 for the amount of precipitation (Fig. 2a).This does not necessarily mean that WRF overestimates precipitation almost two-fold.Rather, it most likely reflects undercatch of the gauge (Goodison et al., 1998) orsince the scaled precipitation proves to be useful for the MB modeling (Sect.3.1) -loss of snow on the glacier by wind drift, which has been observed during field work but is not resolved in the MB model.Both processes cannot be quantified at our site and are therefore absorbed by K WRF < 1.
There is a detailed discussion in Sect.3.3 about the influence of K WRF on the results.Figure 3 presents the time series of all local variables that are used to force the MB model, and indicates their source (measurement or WRF output).We will repeatedly refer to details of Fig. 3 later, but in general the plots nicely illustrate the local climatic setting with the coincidence of summertime peaks in air temperature and the moisture variables (Fig. 3a-d).In principle, the MB model could also be forced with atmospheric model output for the AWS1 grid cell alone, and in this context the modeling period could be extended as well.However, in this paper we attempt to use as many insitu measurements as possible, so the study is limited to the three years 2009-2011 (Fig. 3).

SEB/MB modeling
The physical MB model is described in detail in Mölg et al. (2008Mölg et al. ( , 2009a)), thus we only give a brief summary.The model calculates the specific MB from mass gains by solid precipitation, surface water deposition and refreezing of liquid water in the snow pack, and from mass losses by surface melt, sublimation and subsurface melt.The following SEB equation lays the foundation for calculating these mass fluxes: These symbols signify incoming shortwave radiation (S ↓), broadband surface albedo (α), incoming and outgoing longwave radiation (L ↓ and L ↑), the turbulent sensible and latent heat fluxes (QS and QL), the heat flux from precipitation (QPRC), and the heat flux from the ground (QG).An energy flux has a positive (negative) sign when it induces an energy gain (sink) at the surface.The sum of these fluxes yields a resulting flux F , which represents the latent energy for melting (QM) if glacier surface temperature is 273.15K. QG consists of fluxes of heat conduction (QC) and penetrating shortwave radiation (QPS), where the latter -owing to the energy conservation principle -is always an energy sink at the surface.
At the upper boundary, the model is forced hourly by the six variables shown in Fig. 3.Note that Fig. 3 presents daily means, while the hourly forcing naturally encompasses a greater range than visible in the figure (e.g.hourly T a can rise to 10-15 • C in summer).As stated in Sects.2.1 and 2.2, these inputs are provided by local measurements or atmospheric model output at one point on the glacier, and are distributed over the entire glacier using vertical gradients (see below) if the model is run in distributed mode.A fixed bottom temperature is prescribed at the lower boundary in the ice, which we chose at 9 m depth based on the measurements in this region (Sect.2.1).The remaining subsurface layers are at 0.1, 0.2, 0.3, 0.4, 0.5, 0.8, 1.4, 2, 3, 5, and 7 m depth, which yielded a stable solution in all runs.For these model layers the initial temperature profile is specified from the available subsurface measurements (Sect.2.1) by linear interpolation.
A digital terrain model at 60 m resolution (re-sampled from the Shuttle Radar Topography Mission; Rabus et al., 2003) and the 2009 glacier extent (Bolch et al., 2010) constitute the topographic boundary (Fig. 1).A concise overview of the MB model, to illustrate which atmospheric forcing variables affect each energy and mass flux in the model, is provided in Table 1 in Mölg et al. (2009a).
Since the last model version (Mölg et al., 2009a) we added a few new features based on published work: (i) the local air temperature gradient can vary between day and night/morning to better capture the effect of katabatic wind development (Petersen and Pellicciotti, 2011); (ii) snow settling and compaction is simulated from the viscous fluid assumption (Sturm and Holmgren, 1998), and together with refreezing can contribute to snow densification.Specifically, we use Eqs.( 5)-( 8) in Vionnet et al. (2012); (iii) the transition solid-liquid precipitation is specified by a temperature range using linear interpolation (e.g.Mölg and Scherer, 2012), instead of abruptly by one air temperature threshold; (iv) a second stability correction for stable conditions is available so the turbulence damping factor can be calculated from either equation 11 or 12 in Braithwaite (1995); (v) the energy flux from precipitation may be included in the SEB (standard equation; e.g.Bogan et al., 2003); and (vi) a second parameterization of L ↓ from air temperature, water vapor pressure and cloud cover (Klok and Oerlemans, 2002) provides an alternative to the original formulation (Mölg et al., 2009b).
In light of the main goal to unravel the "monsoon footprint" in the glacier SEB/MB, we also attempt to quantify the model uncertainty.Thus, we perform a Monte Carlo simulation consisting of 1000 realizations, in which all important model parameters, including the vertical gradients for extrapolation of meteorological conditions, as well as selected structural uncertainties, are varied (Appendix A).This simulation is conducted for one point at the location of AWS 1, since the computational expense for the distributed model is too large.From the point results we select three combinations of parameter settings that are maintained for running the full distributed MB model (Appendix A), i.e. the final ensemble size is 3.One of these combinations, henceforth called reference run (REF), reflects the best match with the measured surface height change.The other two settings are used to spread the uncertainty range around REF (Appendix A).All error estimates in the remainder of the text are based on this method, unless otherwise noted.One advantage of the method is that the uncertainty does not simply increase with progressive simulation time like in classical sensitivity experiments where one parameter is changed and all others are held constant.Instead, in the present method the uncertainty range is mainly dependent on surface conditions (snow vs. ice), which is shown later in the model evaluation (Sect.3.1) and in Appendix A.

Characterization of monsoon dynamics
The specific component of the Asian monsoon system that affects the Tibetan Plateau is the Indian Summer Monsoon (ISM) (Ding, 2007).To relate glacier energy/mass fluxes to the ISM we focus on (i) interannual variability associated with the monsoon onset and cessation, and (ii) intra-seasonal variability tied to the active/break cycle -both of which are typical oscillations in the monsoon dynamics (Webster et al., 1998;Ding, 2007).To characterize these dynamics we follow Prasad and Hayashi (2007) in identifying active and break periods of the ISM, since the authors showed that their method satisfies both precipitation-and circulation-based indices used in previous research.A horizontal wind shear index (HWSI) is defined as the difference in the 850 hPa zonal wind (from NCEP/NCAR reanalysis data: Kalnay et al., 1996) between a southern (5-15 • N, 40-80 • E) and northern region (20-30 • N, 70-90 • E).Active periods are defined as days with HWSI > 1σ over 2009-2011, and break periods as days > 1.5 σ in the northern region's zonal wind (i.e.strong westerlies).
Figure 4a shows that establishment of the ISM circulation in the examined years occurred between 1 June (2011) and 21 June (2009), when the curves start to rise.The ISM ceased between 16 September (2011) and 19 October (2010), when the curves flatten.The former interval contains the mean onset date in terms of convection and rainfall in the Zhadang region (Webster et al., 1998).Thus, we define the yearly period of monsoon onset (monsoon cessation) from 1-21 June (16 September-19 October), which comprises the time windows used in the subsequent analysis (Sect.3.4).Only two years are available for the analysis of the cessation period, as AWS data in 2011 ends before 16 September.However, both 2009 and 2010 contain the interval 1-19 October and cover the respective cessation dates (6 October in 2009 and 19 October in 2010).
Between 22 June and 15 September (ISM core season) the succession of active periods and the remaining non-active (weak and break) periods characterizes the intra-seasonal ISM variability (Prasad and Hayashi, 2007).Consequently, we have 105 active days and 39 non-active (weak/break) days during the core season in our data set.Amongst the non-active periods, 10 days can be classified as break days, which will also be considered in the analysis (Sect.3.4).These amounts of days are typical of the active/break cycle in the monsoon region (Webster et al., 1998).
Late monsoon onset and few active periods in 2009 (Fig. 4a) are in agreement with the weak cumulative strength of the ISM circulation in this year, as seen in Fig. 4b.To draw robust conclusions about monsoon activity, it is desirable to supplement the circulation index with a convection indicator (Wang and Fan, 1999), for which we extract topof-atmosphere outgoing longwave radiation (Liebmann and Smith, 1996) over the ISM precipitation region.This second indicator yields the same result (Fig. 4c): 2010 and 2011 exceed 2009 in terms of ISM strength, which is particularly obvious from day 200-270 (ca.July-September) where 2009 shows markedly reduced convection (i.e. higher values in Fig. 4c than other years).

MB model performance
We first evaluate the model at the point scale against observations at AWS1 in Fig. 5.The observations of a rather stable surface height during the winter months and surface lowering in summer are well reproduced (Fig. 5a).Regarding the latter, the model also captures the differential ablation between 2009 (strong) and 2011 (weak) well.Figure 5a also illustrates that the model uncertainty is greatest when snow cover is removed during periods of strong ablation (first half of the 2009 curve).This is a reasonable finding since modeling snow ablation, which is complicated by variable snow density and refreezing processes, is more difficult than modeling ice loss (e.g.Mölg et al., 2009a).The root-mean-square difference (RMSD) in Fig. 5a is only a tiny fraction of the observed amplitude (≈ 3 m), and the explained variance is very high.
A similarly-good performance is evident for glacier surface temperature (T sfc ), where the model captures the variability between 245-250 K mean daily T sfc in winter and melting point in peak summer well (Fig. 5b).Measured T sfc is based on emitted radiation (Table 1), which typically leads to uncertainties over glacier surface on the order of 2-2.5 K (e.g.Greuell and Smeets, 2001;Mölg et al., 2008).The RMSD in Fig. 5b, however, is below 2 K and corroborates the model.Moreover, the temperature variability in the subsurface (Fig. 5c) indicates that penetrating shortwave radiation is simulated reasonably well, since no systematic model bias is evident (Hoffman et al., 2008;Mölg et al., 2008).Finally we consider global radiation (an important driver of the SEB) in Fig. 5d, which shows good agreement for variability but a rather high RMSD.As soon as the two curves are a 5-day smooth, however, the RMSD drops below 10 % of the mean measured global radiation, and 10 % is indeed a more realistic estimate for measurement uncertainty in remote field places than the nominal accuracy (e.g.Michel et al., 2008).Thus, incoming shortwave radiation generated by the MB model can be interpreted reliably for means over five or more days, an averaging period always used in this study for the monsoon analysis.
To evaluate the distributed output of the model, we calculate the mean MB over the available ablation stakes as well as the associated gradient of the vertical balance profile (VBP) in in simulation period 3 except for the last 13 days in September 2010 (which are missing in the model).Since ablation late in September is usually small, we neglect these days.However, for the second interval of stake readings (3 September 2009-17 May 2010) the model is missing the entire month of September, since AWS data in period 2 starts on 1 October 2009.As ablation in early September can still be large, we run the MB model for September 2009 with only atmospheric model output as MB model forcing in order to better evaluate the results for the second interval (without September 2009, modeled MB has a positive bias).The only discrepancy in Fig. 6a concerns the VBP in the 2010 interval, where the model shows a higher gradient.However, there is agreement for the net mass flux in the same interval, which is most important as the area-integrated mass is of primary interest in this study.In all other cases, model and measurements agree within the error bounds (Fig. 6a), which gives us confidence that ablation processes are well simulated.Correlation coefficients between the single stakes and the associated model locations range between 0.5 and 0.81 in the three intervals of Fig. 6a (significant at 5 % based on a two-sided t-test), so the model also captures the basic structure of observed spatial variability.Note that strong mass loss in summer 2009 affected the entire ablation area (Fig. 6a), not only the AWS1 site (Fig. 5a).

SEB/MB characteristics
The altitudinal dependence of MB is shown in Fig. 6b, where areas of positive modeled MB are confined to elevations > 5750 m.The most notable feature of the VBP is the change in slope between 5600 and 5700 m.The only other distributed MB model study for a nearby region, which was also checked against stake measurements, is for Xiao Dongkemadi Glacier in the Tanggula Mountains for October 1992-September 1993.There, Fujita and Ageta (2000) also found steepening of the VBP around 5600 m; however, it occurs above the equilibrium-line altitude (ELA).If we examine only period 3, when ablation on Zhadang Glacier was rather weak (Figs.5a and 6a), the "knee" in the VBP is also shifted along the x-axis and above the ELA (Fig. 6b).
The single year studied by Fujita and Ageta (2000) was also characterized by weak ablation and a slightly positive glacier-wide MB.In the pure modeling study by Caidong and Sorteberg (2010), the steepening of the VBP occurs at ≈ 5800 m.These findings suggest that a dual VBP gradient is a robust feature of glaciers in the central Tibetan Plateau, and is more determined by altitude (weak above 5600-5800 m) than by the MB in a specific year.The result also fits with theoretical work, as the shape/gradient in Fig. 6b is almost a perfect mixture of the typical mid-latitude and subtropical VBP (Kaser, 2001).The explanation of the shape can be found in the energy fluxes: the VBP in Fig. 6b mimics the profile of QM (not shown), i.e. the energy available for melting diminishes clearly above 5700 m.What processes drive the variability of the glacier-wide MB?In general (Fig. 7a; unit is W m −2 ), S ↓ (260.4) and L ↓ (200.9)dominate energy input, followed by QS (17.9),QC (4.5), and very small QPRC (0.2). Reflected shortwave radiation S ↑ (−187.7),L ↑ (−268.2),and QL (−10.9),QPS (−3.4), and QM (−13.7) are the energy sinks at the surface.A salient feature of the variability is the period of April-June both in 2010 and 2011, when S ↑ developed into an equally strong energy sink as L ↑ (i.e. when the two lower curves converge in Fig. 7a).The pattern was completely different in 2009, when anomalously low albedo weakened this process of energy removal.This, in combination with the highest L ↓ in the record, resulted in exceptionally high availability of melt energy in June-July 2009.Also noteworthy is a continuously negative QL on the monthly scale (Fig. 7a), which is not unexpected owing to the generally dry conditions in the central Tibetan Plateau.
The low reflectance of solar radiation in summer 2009 is clearly manifested in the MB components (Fig. 7b).First, surface melt is much higher than in the simulated months of the ablation seasons 2010 and 2011.Second, the high energy availability at the surface also caused more penetrating shortwave radiation and thus subsurface melt in 2009, which hardly occurs at other times.Due to the negative QL discussed above, sublimation is also evident in the MB record.It peaks in the months prior to monsoon onset (e.g.February-April 2010) when (i) T sfc rises after the winter minimum (Fig. 5b) but the atmosphere remains rather dry, which favors a large surface-air vapor pressure gradient, and (ii) higher wind speeds (Fig. 3e) drive turbulence.Melting is absent from November-March, and in total 3.7 % of the modeled grid cell-scale melt happens at air temperatures below 0 • C. The latter was previously detected as a typical feature of glaciers in High Asia (e.g.Aizen et al., 2002).Refreezing of liquid water can occasionally be larger than accumulation by solid precipitation, generally in spring/early summer (Fig. 7b).This feature affirms previous statements that refreezing in the snow is an evident process on Asian high-altitude glaciers (Ageta and Fujita, 1996;Fujita and Ageta, 2000).On average, 13 % of surface melt refreezes, which is less than the 20 % obtained for Xiao Dongkemadi Glacier (Fujita and Ageta, 2000).The glacier-wide MB estimate from the REF run over the 743 modeled days (−1.63 × 10 3 kg m −2 ) is composed as follows (same unit): solid precipitation (1.02), refreezing (0.32), deposition (0.03), surface melt (−2.52), sublimation (−0.28), and subsurface melt (−0.20).The dominance of melt over sublimation at this site fits into the regional-scale pattern of ablation characteristics suggested by simplified MB modeling (Rupper and Roe, 2008).We can also give two annual MB estimates for Zhadang Glacier, from 1 October to 30 September of the subsequent year for consistency with previous studies (Kang et al., 2009).Glacier-wide MB based on the model for 2009/2010 (2010/2011) is −154 ± 43 (−382 ± 41) kg m −2 .The model data gap from 26 June to 15 August 2010 can be accounted for fairly reasonably, since the initial condition for period 3 on 16 August 2010 (based on observations) is a snow-free glacier (i.e.all the snow at the end of period 2 is assumed to be lost in the gap for the MB estimate).Still, 2009/2010 is more likely an over-rather than under-estimation, since ice may have also ablated in this gap.For 2010/2011, the model ends 15 days earlier in September 2011, which must be neglected.The values found here are within the range measured from 2005-2008, −1099 to 223 kg m −2 per year (Kang et al., 2009).usual sensitivity studies that change the value of one internal model parameter.From the calculations, we find little influence of using variable surface roughness lengths, while QPS, topographic shading and refreezing are on the order of the model uncertainty (0.21 × 10 3 kg m −2 ).Also important are snow compaction/settling due to the resultant effects on snow density and subsurface melt.Suppressing the latter in the model leads to a saving of 0.34 × 10 3 kg m −2 in the MB, which is a higher absolute value than the subsurface melt term in the MB budget (−0.20 × 10 3 kg m −2 ; see Sect.3.2).This result indicates that the absence of subsurface melt has a feedback potential (mostly through the influence on snow depth and thus albedo).The strongest impact on MB is clearly from the stability correction of turbulent fluxes (Table 2).A strong feedback process also operates here, since the increase in absorbed shortwave radiation (5.3 W m −2 ) from the initially turbulence-driven acceleration of snow ablation is eventually larger than the increase in the turbulent flux itself ( QS + QL = 3.6 W m −2 ).The stratification of the surface layer, on the other hand, reduces the strength of QS and QL to 54 % of their value in neutral conditions.Thus, not accounting for stability effects can lead to a strong negative MB bias, and Table 2 suggests that any physical MB model for Tibetan glaciers must include stability correction.Also, empirically-based models driven by T a and precipitation alone should incorporate parameterizations of refreezing (e.g.Gardner et al., 2011) for Tibetan glaciers.Refreezing only shows an intermediate effect in Table 2, but due to the systematic seasonal importance in the MB (Fig. 7b) its neglect seems unwarranted as soon as seasonal variations are modeled and interpreted.Furthermore we present in Table 2 classical, static MB sensitivity to constant changes in the two atmospheric forcings T a and precipitation, using the typical perturbations of 1 • C and 10 %.First, in both T a and precipitation perturbations, the MB change is hardly affected by conditions in winter.Hence, our study complements other evidence that the "summer accumulation type glaciers" of Asia are extremely sensitive to atmospheric conditions in the warm season, as discussed in the introduction.Second, the general MB sensitivity is calculated from Table 2 as the mean of absolute MB for negative and positive perturbations before being converted to an annual value by multiplying with the factor 365/743, i.e. the number of days per year/days in model record.The annual value then amounts to (all units in the remaining paragraph in 10 3 kg m −2 = m w.e.) 0.47 per • C for the T a perturbations, and to 0.14 per 10 % for the precipitation perturbations.These numbers can be compared to the few available studies that also employed identical perturbations to a glacier-wide, SEB-based MB model run over at least one year.A typical glacier in the European Alps (Klok and Oerlemans, 2002) seems to respond slightly stronger (≈ 0.67 per • C and ≈ 0.17 per 10 %), while the extremely maritime and high-precipitation glaciers in New Zealand (Anderson et al., 2010) show a clearly higher sensitivity in the order of 2.0 per • C and 0.4 per 10 %.On the other hand, the high-altitude equatorial glaciers on Kilimanjaro (Mölg et al., 2009a) show a lower response in their dry climatic environment (0.24 per • C and 0.09 per 10 %) than Zhadang Glacier.Our calculations therefore support the generally accepted relation of increasing MB sensitivity as the climatic conditions become wetter and warmer (e.g.Fujita and Ageta, 2000).

We performed sensitivity runs for the REF configuration in
Finally, we vary K WRF , which is the only parameter used to scale MB model input (Sect.2.2).This is done for period 1 as it relies entirely on scaled summer precipitation as input (Fig. 3d), shows the highest mass amplitude (Fig. 5a), and hence is most sensitive to the scaling.In concert with K WRF we vary the density of solid precipitation, because (i) it is a free model parameter (Appendix A) and together with K WRF determines the actual precipitation height, i.e. the variable that enters the MB model (Mölg and Scherer, 2012), and (ii) this procedure helps to explore the physically meaningful range of the scaling (see below).Multiplication of K WRF by 1.25, 1.5 and 1.75 changes the glacier-wide MB in period 1  A1; Mölg and Scherer, 2012), we can certainly treat the K WRF • 1.75 case as the upper sensitivity of the scaling approach.However, 143 kg m −2 is lower than the modeled uncertainty for period 1 (217 kg m −2 ), so it can be ruled out that the scaling procedure alters the interpretation of MB model results.

Impact of monsoon variability
The final important question in this study is as follows: what is the role of the Asian monsoon in the SEB and MB processes discussed so far?We begin with the monsoon onset period 1-21 June.A series of systematic differences in the MB components, which are not affected by model error, is evident in Fig. 8a: the later the onset of the ISM, the (i) less accumulation by solid precipitation, (ii) more sublimation, (iii) more melt at the surface, and (iv) more subsurface melt as well.The overwhelming signal is the sharp increase in surface melt when the ISM commences late, as in 2009.and the higher latent energy flux of melting is clearly controlled by low albedo and reduced reflection of solar radiation.Thus, the timing of monsoon onset leaves a clear footprint on the glacier through the albedo effect, which leads to higher mass loss on the glacier the later the ISM is established in a particular year.The monsoon cessation period from 1-16 October is the theme in Fig. 9.The clearest signal again concerns surface melt (Fig. 9a), however, this time in the opposite way: the presence of the monsoon coincides with both increased surface melt and slightly higher sublimation, and therefore with higher ablation on the glacier.Figure 9b illustrates that the higher melt amount is mainly a response to more absorbed shortwave energy.As differences in S ↓, L ↓, and solid precipitation are all small (Fig. 9), the higher albedo in 2009, and therefore the MB in the cessation period, seems to be mostly a delayed response to the snow cover evolution in the preceding core monsoon season: mean modeled snow depth on 1 October is 0.26 m higher in 2009 than in 2010.Thus, while it is difficult to elaborate the ISM impact during its cessation phase from only two years of data, the most important conclusion from Fig. 9 concerns the amplitude of the MB response, which is almost an order of magnitude smaller than during monsoon onset (cf.Fig. 8a).Hence, in general the ISM footprint at the end of the monsoon season is of little importance.

T. Mölg et al.: Monsoon impact on a Tibetan glacier
Figure 10 finally turns to the monsoon core season.The amplitude of mass fluxes (Fig. 10a) is comparable to the monsoon onset phase (Fig. 8a), but a series of systematic differences between active and non-active ISM days is not evident.Surface melt is higher during break periods than during active periods, but this difference does not hold for active versus non-active days (Fig. 10a).Solid precipitation, albedo and L ↓ are higher during non-active and break periods than on active monsoon days (Fig. 10), a surprising finding in the context of the core season given that monsoon air masses carry abundant moisture.We discuss this result in the paragraph below from three angles, but wish to first note that introducing a lead/lag of one or two days to the analysis in Fig. 10 for the monsoon/local SEB and MB relation does not change the results, since active and non-active periods occur as clusters of several consecutive rather than isolated days (Fig. 4a).
First, Prasad and Hayashi (2007) analyzed the synoptic structure of the intra-seasonal ISM variability in detail.Their Fig. 4 shows that convective centers are clearly shifted away from the main ISM precipitation region (defined by Ding, 2007, as 65-105 • E, 5-27.5 • N) during nonactive periods, but little difference is evident over the Tibetan Plateau.An analysis of WRF output as well as satellitebased data (Fig. 11) supports the idea that the expected increase in precipitation during active monsoon days exists for the large-scale ISM area, but not for the Zhadang region.Thus, atmospheric variability over the ISM area does not seem to directly influence the Zhadang region.Second, Ueno et al. (2001) highlighted that precipitation on the Tibetan Plateau is dominated by weak but frequent events in the monsoon season that typically originate from mesoscale systems or local deep convective cells.Hence, the monsoon can be understood as the "background trigger", but local to mesoscale processes forced by the complex relief structure, as well as by the numerous lakes (e.g.Lake Nam Co close to the investigation site: Haginoya et al., 2009), modify the large-scale flow and cause a unique precipitation regime over the Tibetan Plateau.A recent study (Chen et al., 2012) also shows that one of the moisture source regions for the Tibetan Plateau in summer is situated on the plateau itself, and that local water recycling is evident in the same season.Third, it is well appreciated that mid-latitude flow impacts High Asia in the non-monsoon season, but less is known about mid-latitude intrusions during the monsoon.For instance, Ueno (2005) shows that mid-latitude flow anomalies affect the Tibetan Plateau as late as in mid May, so it cannot be ruled out that mid-latitude circulation also modifies ISM activity in the core season, especially when the monsoon onset occurs early.This should be investigated in more detail in future studies in the context of cryospheric changes in High Asia.
In summary, the moisture regime on the Tibetan Plateau during the ISM core season seems to be controlled by local and regional circulations, and/or by the fact that large-scale ISM variability does not advance as far north as onto the central Tibetan Plateau.Therefore, local conditions may show a poor correlation to large-scale flow at this time of the year, e.g. as also found for measured precipitation at Lhasa from May to September (Caidong and Sorteberg, 2010).In this vein, local SEB and MB on Zhadang Glacier do not follow a systematic pattern in relation to large-scale ISM variability during the monsoon core season (Fig. 10).

Conclusions
Our state-of-the-art, distributed and physical mass balance model (Mölg et al., 2009a) reproduces the available measurements on Zhadang Glacier well, and a novel combination of the Monte Carlo and ensemble approaches allows for a reasonable time-varying estimate of model uncertainty.Model forcing is mainly based on field measurements, but also on output from high-resolution atmospheric modeling -a strategy that should become more common for data-sparse regions (e.g.Van Pelt et al., 2012).From the glacier model's output it is evident that interannual variability in mass balance has its physical origin in late spring/early summer, when the energy sink by reflected shortwave radiation is much weaker in high-ablation years than in other years.Affected Table A1.Free parameters in the MB model.Base value (V) and uncertainty (U) are from the literature, or constrained by on-site field measurements (meas.-Sect.2.1) and atmospheric model output (atmo.-Sect.2.2).For assumptions (assum.), the uncertainty is chosen to be relatively large (20 %).Bold parameters indicate structural uncertainty, others parametric uncertainty (e.g.Tatang et al., 1997).For the Monte Carlo simulation, values are assigned randomly from a uniform distribution, except for case 13 (normal distribution) since the authors provide uncertainty as one standard deviation (Gromke et al., 2011).Uncertainty is the 95 % confidence interval of the gradient in the atmospheric model output.c Several roughness lengths are specified because the MB model uses a scheme for space/time-varying roughness, which is described in Mölg et al. (2009a).d Between the measured initial snow depth at AWS1 and the highest altitude on the glacier in the digital terrain model; there are no measurements other than at AWS1, but field experience clearly suggests increasing snow depth with altitude.e 0.2 K is the amplitude of the measurements at that level (Sect.2.1).f Optimal values of these three albedo scheme parameters are determined from measurements over period 2 (the only period with a reliable snow depth record), as described in Oerlemans and Knap (1998).The uncertainty reflects varying choice of the parameters that are held constant in the optimization procedure.
if the Monte Carlo simulation is done with the point model, vertical gradients influence the result due to the altitude difference between AWS1 and the associated grid cell in the digital terrain model (≈ 20 m).Also, these gradients clearly differ between the three chosen parameter combinations for the distributed runs (i.e.REF and two uncertainty settings), which allows us to capture uncertainty for the higher glacier reaches as well.

Fig. 2 .
Fig. 2. Precipitation in the atmospheric model.(a) Accumulated precipitation between 22 May 2010 and 15 September 2011 at the precipitation gauge close to AWS2 (5566 m) and simulated at the associated WRF grid cell (5533 m).A scaling factor of 0.56 for WRF equals the total measured sum.(b) Mean diurnal cycle in WRF at the same grid cell for the summer (May-September) and winter (October-April) over 2009-2011.

Fig. 3 .
Fig. 3. Daily means of (a) air temperature, (b) relative humidity, (c) cloud cover fraction, (e) wind speed, (f) air pressure, and (d) daily sums of all-phase precipitation at AWS1 between 27 April 2009 and 15 September 2011.Ticks on the x-axes indicate the first day of the respective month, and black (blue) signifies measurements (atmospheric model output).In panel (d), the data gaps due the measurement failure are indicated to distinguish them from zero precipitation, and the measurement source for surface accumulation in winter is indicated as well (SR50 or gauge).For consistency, SR50 derived precipitation (actual height) has been converted to w.e.values by a density of 200 kg m −3 .

Fig. 5 .
Fig. 5. Measurements versus MB model (REF run) at AWS1 with correlation and root-mean-square statistics.(a) Accumulated surface height change (error bars reflect model uncertainty as defined in Sect.2.3), where every curve starts at 0 at the beginning of the periods; (b) glacier surface temperature; (c) subsurface temperature, whereas observed depth varies between start and end from 5.6-3.2m in period 1, 5.3-5.1 m in period 2, and 4.4-4.1 m in period 3 (model values found by linear interpolation between layers; note a measurement gap within period 2 in May 2010); (d) global radiation.All data are daily mean values.
Fig. 6.(a) Mean mass balance (black) and vertical balance profile (VBP) gradient (grey) over the available ablation stakes (see Fig. 1) and the associated model locations (found by bi-linear interpolation) for the three intervals indicated at the bottom.Error bars are defined as 1σ for mass balance, and as 95 % confidence interval of the least-squares regression coefficient for the VBP gradient.(b) Modeled VBP over the entire simulation period (743 days) and period 3 (16 August 2010-15 September 2011).Model data in (a) and (b) are from the mean model (REF and two uncertainty settings, see Sect.2.3).

Fig. 7 .
Fig. 7. Glacier-wide mean monthly (a) SEB components with surface radiation terms shown as lines, albedo as dots, and the remaining fluxes as bars (see Sect. 2.3 for abbreviations; precipitation heat flux is not shown as it is negligibly small) and (b) MB components from April 2009 to September 2011 in the REF run.Note the y-axis break and variable scaling in (b) due to large melt amounts in 2009.

Fig. 10 .
Fig. 10.Glacier-wide (a) MB and (b) SEB components in the monsoon core season for active, non-active, and break composites in 2009-2011.Error bars reflect model uncertainty as defined in Sect.2.3.Note the different y-axis scaling in (b).

Table 1 .
Measurement specifications for AWS1 located at 5665 m a.s.l.The height (depth) values refer to the initial distances to the surface on 27 April 2009.Last column indicates the usage for the mass balance modeling: forcing (F), parameter setting (P), or model evaluation (E).The two radiation components yield the measured albedo.
a With ventilated radiation shield.b The months October to April.c Uses the principle of emitted radiation.
Table 2 because little is known about the importance of particular physical processes for the Tibetan glaciers.For these runs, certain structural components of the MB model (physical parameterizations) are deactivated, thus they differ from www.the-cryosphere.net/6/1445/2012/The Cryosphere, 6, 1445-1461, 2012

Table 2 .
Sensitivity of glacier-wide MB over the entire simulation period (with respect to REF settings: −1.63 × 10 3 kg m −2 ) to particular processes or (last four lines) changes in the forcing.Relative changes in parentheses.Winter (W) and summer (S) half years are defined as October-March and April-September, respectively.