We present a new open-source, collaborative “COupled Snowpack and Ice surface energy and MAss balance model” (COSIMA) that is evaluated for Zhadang glacier, Tibetan Plateau. The model is calibrated, run, and validated based on in situ measurements and atmospheric model data from the High Asia Refined analysis (HAR) over the period April 2009 to June 2012. Results for the model runs forced by both in situ measurements and HAR agree well with observations of various atmospheric, glaciological, surface, and subsurface parameters on the glacier. A time-lapse camera system next to the glacier provides a 3-year image time series of the mean transient snow line altitude and the snow cover pattern, which is used for the spatial and temporal validation of the model. The model output corresponds very well to the observed temporal and spatial snow cover variability. The model is then run for a 10-year period of October 2001 to September 2011 forced with HAR data. In general, the radiation components dominate the overall energy turnover (65%), followed by the turbulent fluxes (31%). The generally dry atmosphere on the Tibetan Plateau causes sublimation to be responsible for 26% of the total mass loss. A proportion of 11% of the surface and subsurface melt refreezes within the snowpack.
According to the overall trend of increasing air temperatures on the Tibetan Plateau and its adjacent mountain ranges (TP) for the past several decades, most glaciers on the Tibetan Plateau (Kang et al., 2010; Yao et al., 2012) and in the Himalayas (Bolch et al., 2012; Kargel et al., 2011) are retreating. Generally, the rates of area loss have been increasing in recent years (e.g., Bolch et al., 2012; Yao et al., 2012), but regional patterns are contrasting (Kääb et al., 2012). Understanding the influence of the spatial and temporal heterogeneity of climate and climate variability (Kang et al., 2010) on glacier melt requires the analysis of the surface energy balance (SEB) and mass balance (MB). In the remote high-mountain regions of the TP the collection of field data is difficult and only few detailed SEB/MB studies exist (e.g., Aizen et al., 2002; Jiang et al., 2010; Mölg et al., 2012; Yang et al., 2011; Zhang et al., 2013). In most studies, observation and/or simulation periods are short or limited to point SEB.
In this study, a newly developed “COupled Snowpack and Ice surface energy and MAss balance model” (COSIMA) is presented. We use different data sources to calibrate and validate COSIMA and to calculate glacier-wide SEB and MB of Zhadang glacier (southern central TP) for the period 2001–2011. On Zhadang glacier, detailed climate observations are conducted since 2009. The applied structures and parameterizations within COSIMA have been used in a number of SEB studies (e.g., Anderson, 1976; Mölg and Hardy, 2004; Mölg et al., 2012; Reijmer and Hock, 2008; van Pelt et al., 2012; Zhang et al., 2013). However, COSIMA was developed independently from the existing models. Compared to other commonly used SEB and MB approaches (e.g., Hock and Holmgren, 2005; Klok and Oerlemans, 2002; Pellicciotti et al., 2009) the parameterization of subsurface energy and mass fluxes within COSIMA is more sophisticated and directly coupled to the surface processes. COSIMA explicitly calculates meltwater percolation, the refreezing process within the snowpack under consideration of latent heat release and resulting subsurface melt, as well as the effects on subsurface temperature, snow density, and ground heat flux. Compared to the SEB/MB model of Mölg et al. (2012), the snowpack is not treated as a bulk medium but all subsurface processes are resolved in a vertical layer structure. The use of high-resolution atmospheric model data generated by dynamical downscaling to force COSIMA over 10 years allows a better understanding of both inter- and intra-annual SEB and MB dynamics. To support the application and further collaborative development of COSIMA, we provide the complete source code on an open access Git repository ( https://bitbucket.org/glaciermodel/cosima).
The snow line separates two glacier parts with different albedo and therefore SEB processes (Hock and Holmgren, 2005). Thus, the accurate estimation of the snow line's temporal evolution is an important factor for melt modeling, especially for the so-called “summer-accumulation type” (SAT) glaciers. Time-lapse photography has proved to be a valuable method to infer the snow accumulation distribution at glaciers on the Swiss Alps (Farinotti et al., 2010; Huss et al., 2013). For the TP, the time-lapse camera system installed at Zhadang glacier is the first of its kind. In this study, we use the inferred mean transient snow line altitude and the snow line pattern from the 3-year image time series to test its applicability for distributed MB model evaluation in a region of SAT glaciers.
The goals of this study are (1) to present and validate the newly developed “COupled Snowpack and Ice surface energy and MAss balance model” (COSIMA) by means of different data sources; the complete source code is provided open access; and (2) to test the applicability of transient snow lines derived from time-lapse camera images at SAT glaciers.
Zhadang glacier is a polythermal (Huang, 1990; Shi and Liu, 2000) northwest-exposed valley glacier (∼2 km2, as of Zhang et al., 2013) at the northern slope of the western Nyainqentanglha Shan (30°28.57′N, 90°38.71′E; Fig. 1). The glacier area ranges from 5553 m a.s.l. to 5947 m a.s.l. The region is influenced both by the Indian monsoon and mid-latitude westerlies (Kang et al., 2009; Mölg et al., 2012). Pronounced precipitation maxima in summer and minima in November and December are characteristic features of the southern central and eastern TP (Maussion et al., 2014). Glaciers in this region belong to the SAT glaciers (Ageta and Higuchi, 1984; Maussion et al., 2014) that show a strong sensitivity to air temperature variations during the monsoon season (Zhang et al., 2013). Both accumulation and ablation season coincide in summer, making the albedo feedback very important for glacier MB (Fujita, 2008a, 2008b; Fujita et al., 2007).
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 (Tair) 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, Tair and RH were linearly interpolated to the 2 m standard height to be consistent with the HAR data. Tair 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 (Rnet), longwave outgoing radiation, and incoming (SWin) and reflected shortwave radiation (SWout; 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.
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).
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 (LWin) from measurements of all other radiation fluxes.
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.
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.
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.
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.
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):Equation 1) to be positive. In this case, F is equal to Qmelt. If Ts is below the melting point, F is zero. To calculate F, COSIMA computes the atmospheric energy fluxes at the glacier surface from meteorological variables (Figs. 2 and 3) and the subsurface temperature profile inside the glacier to determine QG. F and Qlat are then converted into mass fluxes that contribute to the surface MB of the glacier together with calculated accumulation. Ts is a key variable in the model and links the surface and subsurface modules. Ts is necessary to solve Equation 1, as it is used to calculate LWout, Qsens, Qlat, and QG. Ts is calculated iteratively through Equation 1 from the energy available at the surface. The equilibrium between the SEB fluxes is required. In cases where Ts exceeds the melting point, it is reset to 273.15 K and the remaining energy flux F equals Qmelt.
For the application of the model at the location of AWS1, SWin is the only SEB component that comes from direct measurements. In the distributed model run, SWin 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 SWin at any pixel of the DEM (x,y; SWin,x,y,pot) is then corrected for cloud cover through SWin measured at AWS1 (or HAR SWin) as follows: For each time step potential SWin at the location of AWS1 (SWin,x_AWS,y_AWS,pot) is set to measured SWin (SWin,x_AWS,y_AWS) or HAR SWin (SWin,x_AWS,y_AWS,HAR) at the location of AWS1. Then, at each time step we calculate the ratio (ri) of SWin,x,y,pot to SWin,x_AWS,y_AWS,pot:
The spatially distributed cloud-corrected shortwave irradiance SWin 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, tsnow 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 SWin 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). LWin and LWout are obtained by the Stefan-Boltzmann law. For LWin, the atmospheric emissivity ϵ is calculated after Klok and Oerlemans (2002):2002). N is directly taken from HAR except during the period 1 January to 10 June 2012 that is not covered by HAR data (Fig. 2, part d). During the latter period, we determine N by comparing SWin to the theoretical values of the top of atmosphere solar irradiance (SWTOA) using the expression after Favier et al. (2004). SWTOA is calculated through the radiation model after Kumar et al. (1997) considering solar constant and geographical position.
Turbulent heat fluxes Qsens and Qlat 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, Tair, and specific humidity at 2 m; cp is specific heat capacity of air (1004.67 J kg-1 K-1); LE is latent heat of evaporation (2.514 × 106 J kg-1); LS is latent heat of sublimation (2.849 × 106 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 (Cse and Clat, respectively) are found to differ by less than 5% (Oerlemans, 2001). Therefore, we assume Cse = Clat. Cse is calculated after Braithwaite (1995):Table 2, Mölg et al., 2009).
Similar to Mölg et al. (2012), the surface roughness length z0 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 z0 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 Tair itself is absolute temperature in K.
In order to calculate QG, vertical profiles of subsurface temperature (Tsub) 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, Tsub is solved from the thermodynamic heat equation (e.g., Anderson, 1976).
QG is the sum of the conductive heat flux (QC) and the energy flux from penetrating shortwave radiation (QPS). QPS is calculated following Bintanja and van den Broeke (1995). ByTable 2), as determined by Bintanja and van den Broeke (1995).
QC 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 Tair (= Ts ) and Tb 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 Tair and Tb 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 QPS. 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):Reijmer and Hock, 2008). For K and E we take the values as determined by Herron and Langway (1980) for two stages of densification: for ρ< 550 kg m-3 we take K = 11 and E = 10,160, for 550 kg m-3 < ρ< 800 kg m-3 we take K = 575 and E = 21,400. After the densification process, the temperature increase through QPS is calculated. The available amount of surface and subsurface meltwater percolates downward and a small amount is held in each layer. We choose constant percolation velocities of 0.06 cm s-1 for unsaturated layers and 0.042 cm s-1 (70%) for saturated layers where the liquid water content (w) exceeds the irreducible water content (wi ). These values are within the boundaries for water film flow and tubular flow through homogenous snow layers (Baumgartner et al., 1990). The wi is assumed to be 5% of the total mass of the layer (Fujita and Ageta, 2000). Then w is calculated and depending on the cold content (γ) of the layer, the available water can refreeze. The γ is defined as the amount of energy that is required to raise the temperature of the respective layer to the melting point and is calculated as follows (Cogley et al., 2011):
Free parameters in COSIMA for the application at Zhadang glacier.
Therefore, the amount of refreezing is limited by w and by Tsub, 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 Tsub and ρin the respective layer.
When Tsub exceeds the melting point in the model, Tsub 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 wi, 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 (Uday ) is calculated after:Table 3). For the application of Uday to glacier-wide annual MB values, annual uncertainties are calculated from the daily means. The results are ±650 kg m-2 yr-1 for the AWS-forced model and ±600 kg m-2 yr-1 for the HAR-forced model. The estimate of uncertainty is much larger than in Mölg et al. (2012, 2014; ±105 kg m-2 yr-1), which probably results from the fact that our estimate of uncertainty explicitly incorporates the uncertainty of the fully distributed spatial model other than using the uncertainty at the location of AWS1 as a measure of uncertainty for the total glacier mass balance.
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).
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.
Results and Discussion
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 LWin, albedo, surface height, Ts, and Tsub 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 LWin and 0.096 for α are within the sensors' accuracy (Table 1) and confirm the applied parameterizations for the calculation of LWin and α within COSIMA. LWin is not directly measured at AWS1 but calculated as a residual from measured net radiation (Rnet), Ts, SWin, and SWout. The uncertainty of LWin is calculated as additive error from Rnet (20%), SWin (5%), SWout (5%), and Ts (0.3 K, corresponds to 3% for Zhadang glacier; Table 1) to ∼21%. SWin and SWout yield measured α. The considered sensors' accuracy is 7% of daily totals (Table 1). At lower values of measured LWin during winter, model values are generally too large (Fig. 4, part b). LWin is calculated from Tair, 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 LWin in winter compared to the summer months. Higher modeled LWin implies that Tair , N, and/or RH taken from HAR are higher than the respective conditions at AWS1 in winter. However, Tair and RH measured at AWS1 and HAR data agree very well. Further uncertainties in modeled LWin are introduced through the applied parameterization of the atmospheric emissivity. The uncertainty of LWin can also vary seasonally depending on, for example, changing surface properties and emissivity. Therefore, discrepancies in the seasonal cycles of measured and modeled LWin may originate either in the model or in the observations.
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 Ts 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 Tsub (Fig. 4, part e). To evaluate model performance, we selected measured Tsub 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 Tsub 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 Tsub from the subsurface temperature profile. An RMSE of only 0.15 K indicates a reasonable parameterization of QG 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 Ts, and 0.16 K for Tsub. 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.
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.
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).
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), SWin (s: +252 W m-2 / w: +212 W m-2) and LWin (s: +244 W m-2 / w: +174 W m-2) dominate energy input over the considered period, followed by Qsens (s: +28 W m-2 / w: +21 W m-2). Energy sinks at the glacier surface are LWout (s: -282 W m-2 / w: -238 W m-2), SWout (s: -196 W m-2 / w: -121 W m-2), Qlat (s: -20 W m-2 / w: -49 W m-2), Qmelt (s: -22 W m-2 / w: -4 W m-2) and QG (s: -7 W m-2 / w: -0.3 W m-2). SWnet (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 SWnet 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 SWnet 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 Qmelt 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 SWnet, 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 Qlat (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, Qlat is usually larger than Qsens 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, SWnet accounted for 38% of overall energy turnover. The energy turnover is calculated as the sum of energy fluxes in absolute values. LWnet accounted for 27%, followed by Qlat (16%), Qsens (15%), and QG (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 Qsens 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.
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).
The SEB and MB study at Zhadang glacier aims at introducing the newly developed “COupled Snowpack and Ice surface energy and MAss balance model” (COSIMA) and furthermore to verify its performance. The installed complex monitoring system including a time-lapse camera is the first of its kind on the TP. Inferred mean transient snow-line altitudes and the snow-line patterns obtained from the orthorectified time-lapse images proved to be a valuable data set for the validation of COSIMA at SAT glaciers. The method of time-lapse photography provides an excellent and inexpensive data source with a high acquisition frequency so far not exceeded by civil high-resolution satellite imaging systems. The evaluation of the AWS and the HAR forced model with various atmospheric, glaciological, surface, and subsurface data sets for the point location of AWS1 as well as in its distributed version for the ablation stake network and the transient snow lines indicates that the MB model captures all relevant SEB and MB processes. For the 10-year period 2001–2011, the energy flux at the glacier surface is dominated by radiation (65%). The generally dry atmosphere on the TP causes sublimation to be an important mass loss component for Zhadang glacier, especially in winter. However, effective melt (74%) is still the dominating mass loss, followed by sublimation (26%). Refreezing compensates for 11% of the surface and subsurface melt.
This work was supported by the German Federal Ministry of Education and Research (BMBF) Programme “Central Asia—Monsoon Dynamics and Geo-Ecosystems” (CAME) within the WET project (“Variability and Trends in Water Balance Components of Benchmark Drainage Basins on the Tibetan Plateau”) under the codes 03G0804E, 03G0804A, and 03G0804F, and by the German Research Foundation (DFG) Priority Programme 1372, “Tibetan Plateau: Formation—Climate—Ecosystems” within the DynRG-TiP (“Dynamic Response of Glaciers on the Tibetan Plateau to Climate Change”) project under the codes SCHN 680/3-1, SCHN 680/3-2, SCHN 680/3-3, SCHE 750/4-1, SCHE750/4-2, SCHE 750/4-3, BU 949/20-1, BU 949/20-2, and BU 949/20-3. Maussion acknowledges support by the Austrian Science Fund (FWF project P22443-N21). Thanks to Thomas Mölg and Peter Kuipers Munneke for their support in the MB model development. We thank Tobias Krenscher and Tino Pieczonka for their support in the development of the time-lapse camera system and the image processing. We thank the staff of the Nam Co monitoring station from the Institute of Tibetan Plateau Research, Chinese Academy of Sciences, for leading the glaciological measurements on Zhadang glacier and for providing ablation stake measurements. We also thank the local Tibetan people for their assistance during fieldwork. We thank two anonymous reviewers and the associate editor, Hester Jiskoot, University of Lethbridge, for valuable suggestions that improved the manuscript.