IN SITU MEASUREMENTS

In a first step, we use observational data from an automatic weather station (AWS1) to force COSIMA (Table 1). AWS1 has been operated in the ablation zone of Zhadang glacier at 5665 m a.s.l. since 2009 (Figs. 1 and 2). Because of the extreme environment, two AWS failures occurred in late summer 2009 and 2010, respectively. Three periods are sufficiently covered by AWS data: (P1) 27 April 2009–14 July 2009, (P2) 1 October 2009–1 July 2010, and (P3) 16 August 2010–10 June 2012. To obtain a continuous data series from 2009 to 2012, these gaps were filled with atmospheric model data from the High Asia Refined analysis (HAR; Maussion et al., 2011, 2014). Except for the period 4 October 2010–30 June 2011 within P3, snow accumulation recorded by a sonic ranger (SR50) is directly used as model input. The SR50 is drilled into the glacier ice to record absolute elevation changes of both snow and ice surfaces. In the following, we refer to this surface elevation change as surface height change. To convert actual snow height to water equivalent, a density of 250 kg m^-3 is chosen for freshly fallen snow following Mölg and Scherer (2012). The missing period was characterized by weak ablation, and at that time the sensor at AWS1 was buried by snow. We used SR50 data from nearby AWS2 (Figs. 1 and 2, part g) which has been operated since 2010 at 5550 m a.s.l. in front of the glacier instead. Subsurface temperatures have been measured at AWS1 until a maximum initial depth of ∼8.5 m, where ice temperature is almost constant (4.7 °C) and only 11% of the annual air temperature amplitude remains. This value is used to define the lower boundary condition (Tb) of COSIMA. Further, temperature measurements at the surface and in ∼1.5 m, ∼4.5 m and ∼5 m depth served for initialization of the temperature profile and the latter also for model evaluation. Following Mölg et al. (2012), the actual height of the air temperature (T_air) and relative humidity (RH) sensor, as well as the depth of the subsurface sensors, was estimated from the SR50 record. From sensors in two heights, T_air and RH were linearly interpolated to the 2 m standard height to be consistent with the HAR data. T_air at AWS1 and AWS2 was used to calculate an annual mean linear vertical gradient (-0.007 K m^-1) that is applied in the distributed model version.

Longwave incoming radiation is used for model evaluation. It is not directly measured at AWS1 but calculated as a residual from measured net radiation (R_net), longwave outgoing radiation, and incoming (SW_in) and reflected shortwave radiation (SW_out; Table 1). Longwave outgoing radiation is obtained from the measured surface temperature by the Stefan-Boltzmann law.

The AWS data were processed according to the instruments' manuals (e.g., temperature correction for SR50, CS215, and IRTS-P, wind correction for NR-LITE, see Table 1) and screened in order to detect and correct measurement errors. Mölg et al. (2012) use the same AWS data but for a shorter period (April 2009 to September 2011). No further corrections for radiative sensor heating and/or riming were conducted, as Mölg et al. (2012) identified only 0.9% of wind data and between 0% and 0.7% of the other variables to be erroneous.

FIGURE 1.

Zhadang glacier (gray area) with glacier outline from 2009 and locations of available measurements and fields of view of the cameras. Glacier outline after Yu et al. (2013), Zhang et al. (2013), and Zhou et al. (2010). Contours are in meters a.s.l. (100 m spacing); coordinate system is UTM 46N, WGS84. The inset shows the location of the study site within China and the Tibet Autonomous Region (gray).

The amount of snowfall, calculated from the SR50 record, leads to the development of an unrealistically high snowpack when directly serving as model input. Therefore, we corrected the calculated snowfall by setting values <1 mm to 0, which is below the measurement accuracy of the sensor (Table 1).

TABLE 1

Measurement specifications for automatic weather station AWS1 located at 5665 m a.s.l. The last column indicates the usage for the mass balance (MB) modeling: forcing (F), parameter setting (P), or model evaluation (E). The two shortwave radiation components yield measured surface albedo. Net radiation is used to calculate incoming longwave radiation (LW_in) from measurements of all other radiation fluxes.

FIGURE 2.

Daily means of (a) incoming shortwave radiation, (b) air temperature, (c) relative humidity, (d) cloud cover fraction, (e) wind speed, (f) air pressure, and (g) daily sums of solid precipitation at automatic weather station 1 (AWS1) between 1 January 2009 and 30 June 2012. Measurements are indicated by black color; atmospheric model output from High Asia Refined analysis (HAR) is gray. In panel (g), the different location of measurements for surface accumulation between 4 October 2010 and 30 June 2011 is indicated (AWS2). The solid proportion of HAR precipitation is scaled by 0.56 in order to match measurements and converted to actual snow height applying a density of 250 kg m^-3.

A total of 25 ablation stakes in the eastern part of the glacier serve for model calibration and evaluation (Fig. 1). Measurements are available for 16 intervals within the ablation seasons between May 2009 and October 2011. The intervals range between four days and seven months. Snow pits were dug at AWS1 in 2009, 2010, and 2011 in order to measure vertical density profiles that serve for density assumptions within the model and for model evaluation.

TIME-LAPSE PHOTOGRAPHY

In addition to the ablation stakes, we use images of a terrestrial time-lapse camera system installed on a mountain ridge nearby the glacier for distributed evaluation of COSIMA. The system consists of two fixed focus (28 mm) Canon EOS 60D/50D cameras (cam1 and cam2), which were located about 1 km north of the glacier tongue. The field of view covers between 57% (cam2) and 60% (cam1) of the glacier area, respectively (Fig. 1). Images were taken automatically every four hours between 2010 and 2012 (00:00, 04:00, 08:00, 12:00, 16:00, 20:00 Beijing Time [BT]). Cam1 operated from 23 May 2011 until 18 June 2012, cam2 with interruptions between 22 May 2010 and 2 September 2012. For this study, we focused on daily images at 16:00 BT (14:03 local time) during the melting seasons in order to detect the spatial pattern and the mean altitude of the transient snow line. The timing of the images was chosen as a compromise between cloud cover formation in the afternoon and the low incidence angle of the sun in the morning (Farinotti et al., 2010; Huss et al., 2013). In order to orthorectify the sequence of images, nine ground control points (GCP; Fig. 1) were measured on rock moraine outcrops at the glacier edge between 2009 and 2011 by means of handheld GPS devices (Garmin).

HIGH ASIA REFINED ANALYSIS

HAR data ( http://www.klima-ds.tu-berlin.de/har/) for Zhadang glacier is available between January 2001 and December 2011. We use hourly HAR data with 10 km resolution (Maussion et al., 2011, 2014) both to fill the gaps in the AWS data from April 2009 to June 2012 and, in a second step, to continuously run COSIMA between October 2001 and September 2011. Hourly means of incoming shortwave radiation (W m-2), air temperature (°C, 2 m), relative humidity (%, 2 m), air pressure (hPa), wind speed (m s-1, 10 m), all-phase precipitation (mm), and cloud cover fraction are extracted from the grid cell that directly covers the glacier area. This grid cell has also the least altitude difference to AWS1 (-39 m) (Fig. 3). The altitude of the HAR grid cell (5611 m a.s.l.) is within the altitude range of Zhadang glacier. Therefore, no statistical downscaling is required. As the cloud cover fraction is not directly available from measurements, this is the only parameter continuously taken from HAR until December 2011.

FIGURE 3.

Monthly means or sums (precipitation) of meteorological variables from the HAR data set at the atmospheric grid cell containing Zhadang glacier (5611 m a.s.l.), October 2001–September 2011. The “COupled Snowpack and Ice surface energy and MAss balance model” (COSIMA) is forced with hourly values of these variables. The scaling factor of 0.56 is already applied on precipitation amounts.

Mölg et al. (2012) compared HAR precipitation with measurements from a precipitation gauge in front of the Zhadang glacier over ∼16 months and found a good correlation with the measured seasonal cycle but obtained a scaling factor of 0.56 for the amount of HAR precipitation. This factor reflects a possible HAR overestimation, as well as measurement errors from the precipitation gauge and/or the loss of snow on the glacier by wind drift. Snowdrift has been observed during fieldwork and is also captured by the time-lapse photography, but it is not integrated in the current MB scheme of COSIMA. Applying this scaling factor to HAR precipitation is suitable for AWS gap filling (Fig. 4) and leads to a good result for the solely HAR-forced COSIMA run (Fig. 5). Annual sums for total precipitation and snowfall from the HAR data set using the scaling factor of 0.56 are within the range of values given by Zhang et al. (2013) for Zhadang.

FIGURE 4.

Measurements and AWS-forced COSIMA results at AWS1. All data are daily mean values. (a) Albedo. (b) Longwave incoming radiation; note that LWin is not directly measured at AWS1 but calculated as a residual from measured net radiation, Ts, SWin, and SWout. (c) Accumulated surface height change; gray squares show average surface height over the three ablation stakes closest to AWS1. (d) Glacier surface temperature. (e) Ice temperature; the sensor depth of ice temperature measurements varies between 5.6 and 2.0 m during the period; in October 2009 the ice temperature sensor has been reinstalled [vertical black dashed line in (e)]; no ice temperatures are available after September 2011.

For the distributed COSIMA runs, altitudinal gradients for relative humidity and air pressure are calculated from four HAR grid cells surrounding AWS1 (Table 2). For wind speed and cloud cover, no distinct altitude dependency could be determined from both AWS measurements and HAR data. The precipitation gradient is derived from the model calibration process.

SEB/MB MODELING

We use a SEB model in combination with a multilayer subsurface snow and ice model (COSIMA) to compute the MB of Zhadang glacier. COSIMA is run on the SRTM (Shuttle Radar Topography Mission) DEM (digital elevation model; Rabus et al., 2003; 3 arc seconds, resampled to grid spacing of 45 m). The complete source code is provided open access ( https://bitbucket.org/glaciermodel/cosima). The two models are coupled in order to account for meltwater percolation, retention, and refreezing. The calculated surface meltwater serves as input for the subsurface model. COSIMA computes the MB at an hourly time step from the sum of accumulation by solid precipitation and refreezing of percolating meltwater in the snowpack, and ablation by melt and sublimation. Forcing parameters are incoming shortwave radiation (W m^-2), air temperature (K, 2 m), relative humidity (%, 2 m), air pressure (hPa), wind speed (m s^-1, 10 m), all-phase precipitation (mm), and cloud cover fraction. The model consists of several modules for, for example, solving the heat equation, calculating percolation and refreezing, and calculating densification. The modular structure allows replacing single parameterizations easily.

In the distributed COSIMA runs, total precipitation from the chosen HAR grid cell is extrapolated to the total altitude range of Zhadang glacier by applying a constant altitudinal lapse rate. This lapse rate is derived from the model calibration process using the ablation stakes. To determine the amount of solid precipitation that accumulates on the glacier surface, we apply a sinusoidal function (Möller et al., 2007) on total HAR precipitation. The function describes the transition between solid and liquid precipitation in a temperature range between +1 °C and +5 °C (Fujita and Ageta, 2000; Zhou et al., 2010; Mölg et al., 2012). The amount of solid precipitation is considered as added mass in the subsurface density profile with a density of 250 kg m^-3 (Mölg and Scherer, 2012). The redistribution of freshly fallen snow through wind drift is not integrated in the current MB scheme of COSIMA.

To calculate the energy available for surface melting (Qmelt) the SEB model within COSIMA solves the energy balance equation (e.g., Oerlemans, 2001):

For the application of the model at the location of AWS1, SW_in is the only SEB component that comes from direct measurements. In the distributed model run, SW_in is derived from a radiation model after Kumar et al. (1997) that computes clear-sky direct and diffuse shortwave radiation in response to geographical position, altitude, elevation, shading by the surrounding terrain, slope and aspect of the grid cell and albedo of the surrounding terrain. The radiation model runs on the SRTM DEM (3 arc seconds, resampled to grid spacing of 45 m). The spatially distributed potential SW_in at any pixel of the DEM (x,y; SW_in,x,y,pot) is then corrected for cloud cover through SW_in measured at AWS1 (or HAR SW_in) as follows: For each time step potential SW_in at the location of AWS1 (SW_{in,x_AWS,y_AWS,pot}) is set to measured SW_in (SW_{in,x_AWS,y_AWS}) or HAR SW_in (SW_{in,x_AWS,y_AWS,HAR}) at the location of AWS1. Then, at each time step we calculate the ratio (ri) of SW_in,x,y,pot to SW_{in,x_AWS,y_AWS,pot}:

The spatially distributed cloud-corrected shortwave irradiance SW_in at any pixel of the DEM (x,y) is calculated as

The parameterization of surface albedo (α) follows the scheme of Oerlemans and Knap (1998), where α is determined as a function of snowfall frequency and snow depth:

In these equations, t_snow is the time since the last snowfall, t* is a constant for the effect of aging on snow albedo, h is the snow depth and d* is a constant for the effect of snow depth on albedo. The free parameters in the albedo scheme are determined from measurements of SW_in and SWout at AWS1 between 2009 and 2011 and according to Mölg et al. (2012). Final values are albedo fresh snow (α_frsnow) = 0.9, albedo firn (α_firn) = 0.55, albedo ice (α_ice) = 0.3, t* = 6 days, and d* = 8 cm (Table 2). LW_in and LW_out are obtained by the Stefan-Boltzmann law. For LW_in, the atmospheric emissivity ϵ is calculated after Klok and Oerlemans (2002):

Turbulent heat fluxes Q_sens and Q_lat are calculated through the bulk aerodynamic method after Oerlemans (2001) between the surface and 2 m, using Tair, relative humidity (RH), and wind speed (u) data from AWS1 or HAR:

In these equations, ρ_air is air density, calculated from air pressure, T_air, and specific humidity at 2 m; c_p is specific heat capacity of air (1004.67 J kg^-1 K^-1); L_E is latent heat of evaporation (2.514 × 10⁶ J kg^-1); L_S is latent heat of sublimation (2.849 × 10⁶ J kg^-1); qair and qs are specific humidity at 2 m and at the surface, calculated from RH, p (air pressure), and saturation water vapor pressure. RH is assumed to be 100% at the surface. The bulk-transfer coefficients for sensible and latent heat (C_se and C_lat, respectively) are found to differ by less than 5% (Oerlemans, 2001). Therefore, we assume C_se = C_lat. C_se is calculated after Braithwaite (1995):

Similar to Mölg et al. (2012), the surface roughness length z₀ increases linearly from 0.24 mm for fresh snow (Gromke et al., 2011) to 4 mm for aged snow (Brock et al., 2006). If the surface is snow-free, we assume z₀ to be 1.7 mm (Cullen et al., 2007). A correction of turbulent fluxes for stable conditions is carried out after Braithwaite (1995):

After Braithwaite (1995), g is the acceleration of gravity (9.81 m s^-2) and Tair(C) is the air temperature in °C while T_air itself is absolute temperature in K.

In order to calculate Q_G, vertical profiles of subsurface temperature (T_sub) and density need to be known. These are calculated within the subsurface model of COSIMA. The model uses a vertical structure that consists of layers with an equal thickness of 0.1 m in the point model version and 0.2 m in the distributed model runs. Each subsurface layer is characterized by a temperature, density, liquid water content, and depth. In this layer structure with a total depth of 10 m, T_sub is solved from the thermodynamic heat equation (e.g., Anderson, 1976).

Q_G is the sum of the conductive heat flux (Q_C) and the energy flux from penetrating shortwave radiation (Q_PS). Q_PS is calculated following Bintanja and van den Broeke (1995). By

Q_C is determined from the temperature difference between the surface and the two uppermost subsurface layers of COSIMA and depends on the thermal conductivity (λ) of the medium (ice or snow). λ is calculated from the subsurface density after Anderson (1976):

The initial temperature profile is defined from the available subsurface measurements and linear interpolation with depth. Tests showed that model results for most SEB and MB components are the same for an initial temperature profile linearly interpolated between T_air (= T_s ) and T_b as for an initial temperature profile linearly interpolated between the subsurface measurements. The vertical profiles of the temperature distribution of the two initialization schemes converge after a few days. In this study, the temperature profile is initialized with T_air and T_b only. The density profile of the initial snowpack is calculated by a linear interpolation between 250 kg m^-3 for the uppermost snow layer and 550 kg m^-3 for the lowermost snow layer as estimated from the snow pits. The depth of the lowermost snow layer depends on the thickness of the snowpack.

First, the temperature profile is calculated without considering the effects of refreezing, subsurface melt, or Q_PS. In case a snow/firn pack is present, the model then calculates the densification of the dry snowpack through an empirical relation after Herron and Langway (1980):

TABLE 2

Free parameters in COSIMA for the application at Zhadang glacier.

Therefore, the amount of refreezing is limited by w and by T_sub, which cannot be raised above the melting point. Refreezing is largest in early summer, when percolating meltwater reaches cold winter snow layers, and in autumn, when capillary water freezes (van Pelt et al., 2012). When a wet snow layer is present in summer, internal water can refreeze at night. This leads to an increase in T_sub and ρin the respective layer.

When T_sub exceeds the melting point in the model, T_sub is reset to 273.15 K and the remaining energy is used for subsurface melt. This meltwater is added to the water content of the respective layer. If w exceeds w_i, the excess water percolates downward into the next layer until it eventually reaches the glacier ice and runs off. Surface runoff occurs when bare ice is exposed at the surface.

MODEL CALIBRATION AND UNCERTAINTY ASSESSMENT

The formulations within COSIMA contain several unknown or poorly constrained parameters. Mean values for these parameters were obtained from observations at the glacier site or taken from literature (Table 2). Mölg et al. (2012, 2014) conducted extensive uncertainty and sensitivity estimations for a point MB model of similar structure at Zhadang glacier. Mölg et al. (2012, 2014) conducted a Monte Carlo simulation for the AWS1 site on Zhadang glacier by varying 29 parameters of the MB model. When the snow cover is removed during periods of strong ablation, derived model uncertainty is largest. Mölg et al. (2012) explained this finding with the complexity of modeling varying snow density and refreezing processes compared to modeling ice loss. For the sensitivity analysis, Mölg et al. (2012) deactivated certain physical parameterizations within the MB model. The authors found the strongest impacts on MB from the stability correction of turbulent fluxes, from snow compaction and settling, and from subsurface melt. Moreover, Mölg et al. (2012) stressed the importance of incorporating refreezing parameterizations for glaciers on the TP. For further details, we refer to Mölg et al. (2012).

As the distributed COSIMA runs forced with AWS data start in late April, we assume an initial snow depth on the glacier that increases linearly with altitude. Based on snow depth measurements at the ablation stakes, we choose an initial vertical snow depth gradient of 0.0015 m m^-1. In the HAR forced distributed model runs starting at 1 October, glacier-wide initial snow depth is set to zero.

MB measurements at the 25 ablation stakes on the glacier between 2009 and 2011 are used to tune the altitudinal precipitation gradient, both for the AWS- and the HAR-forced distributed COSIMA runs. For calibration, we perform five model runs with each of the AWS- (AWS_r1-AWS_r5; Table 3) and the HAR-forced (HAR_r1-HAR_r5; Table 3) distributed models and randomly leave five stakes out for each model type. In each of the runs r1 to r5 the same five stakes are neglected in both the AWS- and HAR-forced COSIMA run. This results in mean precipitation gradients of 0.046 ± 0.004% m^-1 for AWS model runs and 0.040 ± 0.003% m^-1 for HAR model runs (Table 3). The lapse rates are applied in all following COSIMA runs for Zhadang glacier. These model runs are called AWS-forced or HAR-forced reference runs (AWS-REF, HAR-REF; Table 3). Calibrating AWS- and HAR-forced model runs over all available ablation stakes leads to precipitation gradients of 0.047% m^-1 and 0.042% m^-1 for AWS- and HAR-driven models, respectively.

To determine the uncertainty of COSIMA, we run the model with the calibrated precipitation gradient for the five stakes left out in the respective calibration run. For each of the 16 measurement intervals (j) and every calibration run, the results for MB are averaged over the five stakes. The model uncertainty per day (U_day ) is calculated after:

TABLE 3

Precipitation lapse rates obtained from the calibration runs r1 to r5 for both AWS- and HAR-forced COSIMA; the RMSE for each model run refers to the five stakes left out in each run. In each of the runs r1 to r5 the same five stakes are neglected in both the AWS- and HAR-forced model run. The respective mean lapse rate is applied in the AWS- and HAR-forced reference run (AWS-REF, HAR-REF).

IMAGE PROCESSING

In order to obtain quantitative information, for example, the transient snow line from the sequence of images taken by the two fixed time-lapse cameras, the images had to be geo-referenced and orthorectified (Corripio, 2004). This procedure, including manual snow-line mapping, was carried out applying the software RSG (Remote Sensing Software Graz). In a first step, one image per week and camera was manually orthorectified. The orthorectification requires the definition of a transfer function relating the two-dimensional pixels of the image to the three-dimensional orography as represented by the DEM (SRTM DEM, 3 arc seconds, bicubically resampled to 1 m resolution). The correct scaling functions were defined by adjusting the position of five to six selected GCP to the corresponding locations in the images. The image coordinates of the GCP and therefore the scaling functions changed due to slight instability of the cameras over time due to deformation of permafrost under the tripods. Applying a spline function, the weekly detected GCP positions were interpolated to all available camera images, and used for orthorectification in batch mode. The transient snow line was then derived by means of manual digitization from daily images taken at 16:00 BT. In case the image at 16:00 was unusable because of low clouds or fog, the snow line of a different image of the same day was selected. During the three ablation seasons considered in this study, only 8 days out of 247 were discarded completely.

MODEL EVALUATION THROUGH ATMOSPHERIC AND GLACIOLOGICAL MEASUREMENTS

First, COSIMA is forced with observed hourly data at AWS1 between 24 April 2009 and 10 June 2012 (see Table 1 and Fig. 2). The AWS-forced model performance is evaluated at point scale against observations of LW_in, albedo, surface height, T_s, and T_sub at AWS1 (Fig. 4). These observations did not serve as model input or for model calibration (see Table 1). The root mean squared error (RMSE) of 35.5 W m^-2 for LW_in and 0.096 for α are within the sensors' accuracy (Table 1) and confirm the applied parameterizations for the calculation of LW_in and α within COSIMA. LW_in is not directly measured at AWS1 but calculated as a residual from measured net radiation (R_net), T_s, SW_in, and SW_out. The uncertainty of LW_in is calculated as additive error from R_net (20%), SW_in (5%), SW_out (5%), and T_s (0.3 K, corresponds to 3% for Zhadang glacier; Table 1) to ∼21%. SW_in and SW_out yield measured α. The considered sensors' accuracy is 7% of daily totals (Table 1). At lower values of measured LW_in during winter, model values are generally too large (Fig. 4, part b). LW_in is calculated from T_air, cloud cover (N), and water vapor pressure (calculated from RH; see Equations 6 and 7). On the TP, these parameters are lowest during winter, resulting in a lower LW_in in winter compared to the summer months. Higher modeled LW_in implies that T_air , N, and/or RH taken from HAR are higher than the respective conditions at AWS1 in winter. However, T_air and RH measured at AWS1 and HAR data agree very well. Further uncertainties in modeled LW_in are introduced through the applied parameterization of the atmospheric emissivity. The uncertainty of LW_in can also vary seasonally depending on, for example, changing surface properties and emissivity. Therefore, discrepancies in the seasonal cycles of measured and modeled LW_in may originate either in the model or in the observations.

FIGURE 5.

(a) Observed and modeled mass balance (MB) over the available ablation stakes (see Fig. 1) for three intervals. The calculation of the error bars is explained in the text. Dates are given as dd/mm/yyyy. (b) Observed and modeled MB for the individual ablation stakes and their individual measurement intervals. COSIMA is forced with AWS data and HAR data.

The AWS-forced model can reproduce the seasonal pattern of an invariant or smoothly increasing surface height in winter and surface lowering in varying intensity in summer very well (Fig. 4, part c). An RMSE of 0.049 m for elevation changes of the snow and ice surface is a convincing result. Daily mean measured T_s ranges between -29.7 °C in winter and 0 °C in summer (Fig. 4, part d). Model results are in strong correlation with this variability and an RMSE of 2.12 K is within the measurement uncertainties that range in the order of 2–2.5 K due to the dependency on emitted radiation as also mentioned by Mölg et al. (2008, 2012). Similar good results are evident for T_sub (Fig. 4, part e). To evaluate model performance, we selected measured T_sub at an initial depth of 5.6 m due to the effect of radiative heating of the sensors closer to the surface. The actual measurement depth of T_sub changes according to the varying surface height and the eventual reinstallation of the sensor (Fig. 4, part e). These processes are taken into account for the calculation of modeled T_sub from the subsurface temperature profile. An RMSE of only 0.15 K indicates a reasonable parameterization of Q_G and confirms the good performance of the AWS-forced COSIMA at point scale.

The evaluation of the HAR-forced COSIMA at point scale was performed in the same way as described for the AWS forced model runs. RMSE values are 0.35 m for surface height change, 2.87 K for T_s, and 0.16 K for T_sub. As expected, the AWS-forced model can reproduce the respective parameters to a higher degree than the model forced with HAR data. The main reason for this might be frequency and intensity of precipitation that cannot be exactly reproduced by the 10 km HAR data, especially at smaller scales. The downscaling of the 10 km HAR data to the 45 m SRTM DEM grid using altitudinal lapse rates does not account for dynamic effects of a high-resolution terrain structure. Uncertainties in the timing of precipitation events and in the amount of precipitation from HAR were already mentioned by Maussion et al. (2011). Generally, precipitation processes in high-mountain regions are complex (Barry, 1992). The simulation of precipitation on the TP using an atmospheric model therefore is not trivial (Maussion et al., 2014). Furthermore, permanent weather stations on the TP are scarce and limited to lower altitudes (Qin et al., 2009). This and the difficulty in measuring precipitation in mountainous regions limit the availability of validation data for atmospheric model output.

To evaluate the distributed model output for both AWS-REF and HAR-REF runs, the MB over the available ablation stakes is calculated (Fig. 5). Stake readings are available for 16 intervals within the ablation seasons between May 2009 and October 2011 (Fig. 4, part c). Considering all available intervals leads to an RMSE between measurements and results of the AWS-forced MB model of 0.34 ± 0.05 m w.e. (R = 0.83). The RMSE calculated with daily values is ∼0.007 m w.e. The RMSE between measurements and results of the HAR-forced COSIMA is 0.36 ± 0.03 m w.e (R = 0.73). The overall scattering in the HAR-forced model is larger, once again probably due to uncertainties in modeled precipitation frequency and intensity.

The uncertainty ranges for the model results presented in Fig. 5, part a are calculated from Equation 17. The uncertainties for the measurements are calculated from the assumption of a ±0.05 m uncertainty in every stake reading and a ±20% uncertainty of the respective mean measured MB value due to the uncertainty involved in the calculation of the snow density. For all three intervals, model and measurements agree within the error bounds (Fig. 5, part a). This gives confidence that both accumulation and ablation are well simulated.

The comparison of individual stake readings for all 16 intervals with the respective AWS-REF and HAR-REF model results (Fig. 5, part b) allows a detailed temporal and spatial analysis of the model performance. For approximately 10 stakes, measured MB is around 0.4 to 0.5 m w.e., whereas modeled MB ranges between 1.0 and 2.2 m w.e. (Fig. 5, part b). These “outliers” stem from the longest considered measurement interval from 3 September 2009 to 17 May 2010. AWS- and HAR-forced model runs in these cases overestimate mass loss in the lower glacier regions compared to stake measurements. Discrepancies might be caused by additional snow accumulation in the lower regions of the glacier tongue by wind drift. High measured wind speeds at AWS1 during the winter months (Fig. 2) support this assumption. Changes in the air temperature lapse rate could also be responsible for the overestimated mass loss by COSIMA.

The location of the ablation stakes is strongly biased to the eastern part of the glacier (Fig. 1). This implies that the derived precipitation lapse rates are not valid for all glacier regions. Consequently, the derived model uncertainties might underestimate the total uncertainty because the glacier-wide spatial variability is not well captured.

MODEL EVALUATION THROUGH TIME-LAPSE PHOTOGRAPHY

We used transient snow-line altitudes inferred from the ortho-images for validation of the distributed AWS-REF COSIMA run. In 2012, available AWS data limit the MB model run to 10 June. For days where a snow line exists, we compare the mean altitude of the mapped snow line at 16:00 BT to the respective model output. The modeled snow line consists of those pixels that form the border between snow-covered and snow-free pixels and that lie within the field of view of the camera (Fig. 6). The respective mean altitude is calculated from the underlying SRTM DEM (45 m). In some images, the glacier area within the field of view of the cameras (Fig. 1) was either totally snow-covered or completely snow-free. The latter is the case only in 2010 after the snow line exceeds 5700 m a.s.l.

The AWS-forced COSIMA reproduces the observed evolution of the mean daily snow line altitude very well. The overall RMSE between mapped and modeled snow line altitude is 20.8 m (Fig. 6). This value is supposed to describe reasonable boundaries for the MB model uncertainty. Uncertainties in the manually mapped snow lines can occur from errors in the image processing, the subjective assessment of the spatial snow line pattern, and the low viewing angle of the cameras over some parts of the glacier. To quantify the uncertainty in the manual snow-line mapping, we apply an accuracy of ±5% of the daily snow-covered glacier area adopted from Paul et al. (2013). Paul et al. (2013) estimated the accuracy of glacier outlines derived from remote sensing and from manual mapping. The authors give an uncertainty of 30% of the glacier area for debris-covered glacier parts and 5% for clean ice with sufficient contrast to the surrounding regions. Considering the overall high contrast between glacier ice and snow, we assume an accuracy of ±5% of the snow-covered glacier area for the determination of the snow line at Zhadang glacier. This uncertainty can be translated into an error in the inferred daily snow-line altitude. Daily errors range between 20 and 32 m (Fig. 6) with an average uncertainty of ±26.3 m.

FIGURE 6.

Mapped and modeled mean transient snow line altitude in the ablation season in 2010, 2011, and 2012 at 16:00 Beijing Time (BT), respectively. The uncertainty ranges correspond to the calculated root mean squared error (RMSE) of COSIMA (20.8 m) and the estimated errors in daily mapped snow line altitude (20–32 m; see explanations in the text). COSIMA is forced with AWS data.

FIGURE 7.

Mapped and modeled temporal evolution of the mean transient snow line altitude at 16:00 BT for the three ablation seasons; the uncertainties correspond to the calculated RMSE of COSIMA (20.8 m) and the estimated errors in daily mapped snow line altitude (20–32 m; see Fig. 6). COSIMA is forced with AWS data.

On average, the snow line retreats from the lowest parts of the glacier tongue between May and June (Fig. 7). In 2010, surface melt starts relatively late, in the beginning of June, but is very strong until September (Fig. 4, part c). This is clearly visible in the snow line evolution (Fig. 7, part a). In 2010, the ablation zone develops not until 11 June, but then the mean altitude of the snow line rises quickly until the glacier is completely snow-free. The following months are characterized by intermediate snowfall events and strong melt resulting in rapid changes in transient snow line altitudes. In September the permanent winter snow cover builds up in the course of storm events in autumn both in 2010 and 2011. In 2011, the snow line is visible on 3 June for the first time, around one week earlier than in 2010, and rises very slowly, repeatedly interrupted by snowfall events (Fig. 7, part b). The uppermost areas of the glacier remain snow-covered throughout the whole summer. The rather small ablation in 2011 is also evident in measured and modeled surface height changes (Fig. 4, part c) and in modeled SEB and MB components (Fig. 9). Two examples for the transient snow line on 25 July 2010 and 24 June 2011 are shown in Figure 8. On 25 July 2010, the spatial pattern of the snow line is clearly detectable on the ortho-image and very well reproduced by the MB model. On 24 June 2011, the pattern is characterized by large snow patches that complicate the exact detection of the snow line. The respective model output corresponds to this irregular pattern. The results support that COSIMA estimates the accurate temporal evolution of the snow line to a high degree and captures also spatially differential ablation patterns.

FIGURE 8.

Examples of the direct spatial evaluation of the transient snow line for two days in 2010 and 2011; orthorectified image of Zhadang glacier in the digital elevation model (DEM) (SRTM DEM, 45 m) overlaid with mapped snow line (blue) and modeled snow line (AWS-forced model; red). For camera fields of view see Fig. 1.

Spatial patterns of the transient snow lines in most regions of the glacier tongue are clearly dominated by altitude and therefore mainly by the air temperature lapse rate. The eastern edge of the glacier tongue is an exception with consistently fast increasing snow line altitudes (Fig. 8), probably because of enhanced melt through longwave radiation and sensible heat release from the adjacent rocks. Thermal radiation emitted from valley walls or rocks was suggested to accelerate glacier melt (Paul et al., 2004). Jiskoot and Mueller (2012) calculated an increase of local melt by up to 30% at Shackleton Glacier, Canadian Rockies. However, heat advection from the valley walls was reduced with increasing katabatic wind speed. Observations at Zhadang glacier are consistent with these findings: On smaller glaciers like Zhadang, radiation effects (short- and longwave) from terrain exposure are important factors determining the spatial distribution of glacier melt (Jiskoot and Mueller, 2012). In the uppermost southeastern regions of Zhadang glacier and in the steep areas in the west, topographic shading generally results in long-lasting snow cover and slowly retreating snow lines (Fig. 8).

SEB/MB CHARACTERISTICS

Glacier-wide mean monthly SEB and MB components from the HAR-REF COSIMA run 2001–2011 are illustrated in Figure 9. In both summer (s; April–September) and winter season (w; October–March), SW_in (s: +252 W m^-2 / w: +212 W m^-2) and LW_in (s: +244 W m^-2 / w: +174 W m^-2) dominate energy input over the considered period, followed by Q_sens (s: +28 W m^-2 / w: +21 W m^-2). Energy sinks at the glacier surface are LW_out (s: -282 W m^-2 / w: -238 W m^-2), SW_out (s: -196 W m^-2 / w: -121 W m^-2), Q_lat (s: -20 W m^-2 / w: -49 W m^-2), Q_melt (s: -22 W m^-2 / w: -4 W m^-2) and Q_G (s: -7 W m^-2 / w: -0.3 W m^-2). SW_net (s: +57 W m^-2 / w: +91 W m^-2) is the most important energy source for the glacier, highly depending on α. This was also determined by Mölg et al. (2012) and Zhang et al. (2013) for Zhadang glacier and by Li et al. (2014) in the catchment containing Zhadang glacier. The lowest values of SW_net are evident in or around September when larger snowfall amounts in context with storm events increase α (Fig. 9, part a). In 2006, this structure is completely missing and the glacier SEB is characterized by high SW_net throughout the year (Fig. 9, part a) caused by low α through very little snowfall in November 2005 to March 2006 and June to August 2006 (Fig. 9, part b). The result is that the model returns the largest amounts of Q_melt and therefore surface melt from March to September (Fig. 9). This stresses the impact of the monsoon onset as reported by Mölg et al. (2012) and Kang et al. (2009) but also the importance of winter snowfall events that keep α high. In general, α decreases at the end of each year (increasing SW_net, Fig. 9, part a) because solid precipitation is minimal in that period (Fig. 9, part b). The generally dry conditions on the TP (Aizen et al., 2002) lead to continuously negative Q_lat (Mölg et al., 2012) on monthly scales with largest values in winter when air humidity is smallest and higher wind speeds enhance turbulence (Figs. 2 and 3). In absolute values, Q_lat is usually larger than Q_sens and responsible for significant mass loss through sublimation, especially in winter and spring prior to the onset of the monsoon (Fig. 9). For the total considered period, SW_net accounted for 38% of overall energy turnover. The energy turnover is calculated as the sum of energy fluxes in absolute values. LW_net accounted for 27%, followed by Q_lat (16%), Q_sens (15%), and Q_G (3.3%). These proportions are characteristic for cold and dry climates (Braun and Schneider, 2000; Mölg and Hardy, 2004). The proportions determined by Zhang et al. (2013) for 2009–2011 at AWS1 are similar to these calculations.

From measurements at AWS1, no simple altitude dependency of wind speed could be determined at Zhadang glacier. However, the observation of wind-induced redistribution of snow in some regions of the glacier implies that the wind field is influenced by topography. Snowdrift is not only an important factor for the spatial variability of snow accumulation and ablation but also affects the energy budget of a snowpack (Barral et al., 2014). Increased sublimation by blowing snow raises the moisture content of the near-surface atmosphere. Thus, SEB/MB models that do not integrate parameterizations for snow drift are probably affected by a dry bias. They also overestimate surface sublimation (Barral et al., 2014). These feedback mechanisms have to be considered when interpreting the results of COSIMA for Zhadang glacier. At larger glaciers, katabatic boundary layer processes have been observed to modify the spatial pattern of near-surface air temperature and vapor pressure (e.g., Oerlemans and Grisogono, 2002; Shea and Moore, 2010). The associated changes in the gradients between surface and ambient air influence the calculation of turbulent energy fluxes. Furthermore, the katabatic flow generates turbulence at the glacier surface and increases Q_sens and Qlat (Oerlemans and Grisogono, 2002). Thus, when applying COSIMA to larger glaciers or ice caps, neglecting katabatic winds will introduce an uncertainty in the SEB and MB results and glacier melt may in cases be underestimated.

FIGURE 9.

Glacier-wide monthly (a) surface energy balance (SEB) components (abbreviations are given in the text), and (b) MB components from October 2001 to September 2011 at Zhadang glacier. COSIMA is forced with HAR data.

The glacier-wide MB over the considered period is -1067 ± 600 kg m^-2 yr^-1. In general, sublimation (-395 kg m^-2 yr^-1) is the second largest factor of glacier-wide mass loss after surface melt (-1174 kg m^-2 yr^-1) but clearly dominates ablation in winter (Fig. 9, part b). Subsurface melt plays a smaller role -115 kg m^-2 yr^-1. Solid precipitation (+478 kg m^-2 yr^-1) and refreezing (+139 kg m^-2 yr^-1) contribute to mass gain of the glacier. In total, 11% of surface and subsurface melt refreezes. Mölg et al. (2014) determined a proportion of 18% for the total area of Zhadang glacier in 2001–2011, whereas Zhang et al. (2013) derived a proportion of 6% at AWS1 for 2009–2011. A reason for the discrepancy to Mölg et al. (2014) is probably the treatment of the snowpack as a bulk medium. Thus, the refreezing process is not resolved vertically in Mölg et al. (2014). In COSIMA, all subsurface processes are resolved vertically in a layer structure. The calculations of Zhang et al. (2013) are based on daily mean forcing variables and are restricted to the location of AWS1 in the ablation area of Zhadang glacier. For this location, a lower amount of refreezing is reasonable. At Xiao Dongkemadi glacier a proportion of 20% refreezing was determined (Fujita and Ageta, 2000; Fujita et al., 2007). A reason for the higher proportion can be found in the different glacier types according to their climate setting. After Shi and Liu (2000), Zhadang glacier is a subcontinental (polythermal) glacier type, whereas Xiao Dongkemadi glacier belongs to the extreme continental (cold) type with lower ice temperatures and therefore larger cold content. At Zhadang glacier, effective melt (surface melt + subsurface melt - refrozen water) accounts for 74% of total mass loss followed by sublimation (26%).

Mean modeled glacier-wide MB of Zhadang glacier over 2001–2011 (-1067 ± 600 kg m^-2 yr^-1) is slightly more negative than the value published by Mölg et al. (2014) (-891 ± 105 kg m^-2 yr^-1) for the same period but within the range of uncertainties. The modeled mean equilibrium line altitude (ELA) between 2001 and 2011 is at ≈5850 m a.s.l. Bolch et al. (2010) determined an ELA of 5710 m a.s.l. for Zhadang from the median elevation of the glacier based on the 2001 glacier extent. Kang et al. (2009) and Yu et al. (2013) published annual MB values of Zhadang glacier for three years based on field data that are considerably less negative than the COSIMA output.

Furthermore, we compare the results of COSIMA obtained for Zhadang with other regional MB studies for glaciers in High Asia over a similar time period using satellite data. Gardner et al. (2013) obtained a value of -270 ± 160 kg m^-2 yr^-1 for the southern and eastern TP and -800 ± 220 kg m^-2 yr^-1 for the East Himalayas (2003–2008). Results of Gardelle et al. (2013a, 2013b) range between -220 ± 120 kg m^-2 yr^-1 and -450 ± 130 kg m^-2 yr^-1 (1999–2011) for the Himalayas that correspond to the findings of Kääb et al. (2012) for the period 2003–2008. Combining satellite and aerial images, Bolch et al. (2011) obtained an MB of up to -1960 ± 530 kg m^-2 yr^-1 for single glaciers in the Everest region over 2002–2007. The few existing MB measurements over several years reveal the most negative mass balances in the Himalayas, ranging from -1100 kg m^-2 yr^-1 to -760 kg m^-2 yr^-1 between 2002 and 2010 (Yao et al., 2012). Zhadang glacier MB is thus representing the lower (more negative) end of these regional estimates, probably because of its small size and its leeward location in the western Nyainqentanglha Shan. Smaller glaciers are characterized by shorter response times to climate changes (Bahr et al., 1998; Bolch et al., 2012). Precipitation on the leeward side of the Nyainqentanglha Shan is lower than windward, resulting in fewer and smaller glaciers and higher retreat rates as a response to recent warming (Bolch et al., 2010).