Interactive comment on “ High-resolution interactive modelling of the mountain glacier – atmosphere interface : an application over the Karakoram ” by E

The paper describes a coupled high-resolution mesoscale atmospheric and glacier mass balance modelling system and its application to the Karakoram – although the analysis is focussed on a single glacier, Baltoro. The method builds on existing work, namely, the widely-used Weather Research and Forecasting (WRF) mesoscale atmospheric model and a physically-based glacier mass balance model developed by the authors over recent years, which has principally been applied to glaciers on Kilimanjaro and in Tibet. The key advance and novelty of the approach is that the two models are coupled interactively. Not only does meteorology force the glacier model in the usual


Introduction
Spatially-distributed simulations of glacier surface and climatic mass balance (where the latter term denotes surface plus near-subsurface mass balance following; Cogley et al., 2011) require distributed meteorological forcing; however, obtaining these data is complicated both by the spatial and temporal scarcity of in situ observations and by the "scale mismatch" between the spatial scales represented in atmospheric models and those relevant for surface and climatic MB calculations (e.g.Machguth et al., 2009;Mölg and Kaser, 2011).To overcome these issues, forcing data can be obtained by extrapolation from point measurements by automated weather stations, where available, or interpolation from climate reanalyses and atmospheric model output, using surface-and free-air lapse rates.Surface lapse rates exhibit significant spatial and temporal variability, however, leading to uncertainty in temperature downscaling from altitude changes (Marshall et al., 2007;Gardner et al., 2009;Petersen and Pellicciotti, 2011).In addition, the assumption of linear lapse rates over glacier surfaces may be inappropriate (Petersen and Pellicciotti, 2011) and may under-predict nearsurface temperature over debris-covered regions (Reid et al., 2012).Finally, additional corrections are often required for the poor representation of the strength and spatial variability of processes relevant to mass balance, such as orographic precipitation, in coarse spatial-resolution atmospheric models (e.g.Paul and Kotlarski, 2010;Radić and Hock, 2011).
Dynamical downscaling has been used to address the issue of spatial resolution in the most recent studies to produce climate data at horizontal resolutions of ∼ 18 km Published by Copernicus Publications on behalf of the European Geosciences Union.(Machguth et al., 2009;Kotlarski et al., 2010a,b;Paul and Kotlarski, 2010), ∼ 11 km (Van Pelt et al., 2012) and ∼ 1-3 km grid spacings (Mölg and Kaser, 2011;Mölg et al., 2012a,b) as forcing for distributed alpine glacier surface-and climatic-mass-balance calculations.This approach provides high spatial-and high temporal-resolution atmospheric fields obtained from a physical model, and the increased resolution allows for improved representation of features such as complex topography and orographic precipitation (e.g.Maussion et al., 2011).However, most of these studies required statistical corrections to link mesoscale circulation patterns and meteorological fields simulated by regional atmospheric models to local conditions on the glacier surface.Mölg and Kaser (2011) first showed that, at sufficiently high spatial resolution (∼ 1 km), a regional atmospheric model could be used to force explicit distributed simulations of glacier CMB without statistical corrections at the glacier-atmosphere interface.This approach has since been applied successfully in multiple locations for small glaciers (Mölg and Kaser, 2011;Mölg et al., 2012a,b).
Traditional approaches to simulations of surface and climatic mass balance, including those discussed above, have been one-way, or offline, in which meteorological fields are passed to the CMB model but changing surface boundary conditions due to CMB processes are not fed back into the atmospheric model.Interactively or two-way coupled atmospheric and ice-sheet simulations with simple treatments of ablation have been performed to estimate the paleoclimate and future climate behaviours of the Laurentide and Greenland ice sheets, respectively, with significant alterations to atmospheric circulation, temperature and precipitation resulting from ice sheet evolution (Ridley et al., 2005;Pritchard et al., 2008).Although an initial effort has been made to include "interactive" alpine glaciers in a regional atmospheric model with the subgrid-scale parameterization of Kotlarski et al. (2010b), the influence of two-way coupling on the atmospheric forcing and explicitly simulated surface and climatic mass balance has yet to be assessed for alpine glaciers.
Here, we build on a new, unified and explicit approach to resolving the glacier-atmosphere interface without statistical downscaling (Mölg and Kaser, 2011), through the use of an interactively coupled high-resolution mesoscale atmospheric and physically based CMB modelling system.By allowing changes in glacier surface conditions to feed back on the atmospheric drivers, the model provides a consistent calculation of surface energy and mass fluxes.For the initial application of the coupled model, we simulate the Karakoram region of the northwestern Himalaya (Fig. 1), which is estimated to contain anywhere from ∼ 1250-4000 km 3 of ice, covering an area of ∼ 18 000 km 2 (Bolch et al., 2012).Due to its extensive glaciation, this region presents a high potential influence on atmospheric simulations resulting from the inclusion of feedbacks from alpine glaciers.In addition, Yao (2007) estimates that more than half of the glaciated area is contained in the 15 largest glaciers, thus optimizing the Karakoram for representation in a high-resolution atmospheric model, where the smallest practical grid spacing is on the order of a few kilometers.
The Karakoram is also of interest due to recent evidence of stable or positive mass balances (e.g.Hewitt, 2005;Scherler et al., 2011;Gardelle et al., 2012;Kääb et al., 2012), which contrasts with the general trend of mass loss exhibited by glaciers elsewhere in the Himalaya (Cogley, 2011).However, definitive statements about the mass balance of Karakoram glaciers have been hampered by a dearth of both in situ measurements and information on ice thickness changes.The latter limitation has been partially addressed by recent geodetic studies (Gardelle et al., 2012;Kääb et al., 2012) that support reduced mass loss or even a positive mass balance anomaly in the early 21st century but emphasize the spatial and temporal heterogeneity of recent glacier behaviour.In addition, explicit, physically based, spatially distributed numerical modelling has the potential to clarify the dynamics occurring in this region.
In this study, we first evaluate the regional and the local performance, respectively, of the coupled model against available measurements, by comparing the former with surface albedo and temperature over all glacier surfaces and the latter with meteorological fields and ablation on the Baltoro glacier.We then aim to (1) explore the importance of energy and mass exchanges between the glacier surface and boundary layer on the atmospheric forcing, and (2) assess the ultimate influence of interactive coupling on simulations of glacier mass balance.A final goal of this work is to improve the representation of alpine glaciers in mesoscale atmospheric models by introducing additional, relevant physical processes.

Methodology
The coupled modelling system (hereafter "WRF-CMB") consists of two components: the advanced research version of the nonhydrostatic and fully compressible Weather Research and Forecasting (WRF) mesoscale atmospheric model version 3.4 (Skamarock and Klemp, 2008, Sect. 2.1) and the process-based surface-energy and CMB model of Mölg et al. (2008, 2009, 2012a, Sect. 2.2).The CMB model has been incorporated into the WRF source code as an additional physics option, and, thus, the user may select via runtime ("namelist") options whether the CMB simulation is offline (conventional one-way forcing, with feedbacks only from WRF's land surface model) or interactive (feedback from CMB model to WRF over glaciated grid cells; Sect.2.3).We performed two simulations, one interactive (INT) and one offline (OFF), for the months of June-August 2004, to coincide with a limited number of glaciological and meteorological measurements from the Baltoro glacier available for evaluation (Sect.2.4), with the period of 1-25 June discarded as model spin-up time.Here we use the term interactive to denote surface-atmosphere exchanges through heat, moisture and momentum fluxes only and not through topographic feedbacks, as glacier geometry is held constant over our brief simulation.As a first approximation, we focused on the meteorologically driven fluctuations of mass balance and neglected the influence of debris cover.

Mesoscale atmospheric model
For these simulations, WRF was configured with three nested domains of 33, 11 and 2.2 km spatial resolution, centered over the northwestern Himalaya (D1-3; Fig. 1).By increasing the spatial resolution over the region of interest, the use of multiple grid nesting improves the representation of complex terrain and associated processes such as orographic precipitation, and has been found to increase the simulation skill of WRF for mountain summit conditions (Mölg and Kaser, 2011).Model physics and other settings were selected following the recommendations of the National Center for Atmospheric Research (NCAR) for regional climate simulations with WRF (Table 1; outlined in WRF ARW user's guide).Note that no cumulus parameterization was employed in the highest-resolution, convection-permitting model domain, WRF D3 (e.g.Molinari and Dudek, 1992;Weisman et al., 1997).The range of terrain elevation represented in this domain at 2.2 km resolution is 916 to 7442 m a.s.l., which encompasses the most heavily glaciated altitudes in the Karakoram (∼ 2700-7200 m, as shown in Fig. S2 of Bolch et al. (2012), with the mean basin-wide glacier elevation located at 5326 m).
In this study, WRF was coupled with the Noah land surface model (LSM; Chen and Dudhia, 2001).The land-ice mask was updated using glacier outlines for the Karakoram region based on the glacier inventory of China (Shi et al., 2009) as well as inventories generated by ICIMOD (2007) and GlobGlacier (Frey et al., 2012).Other modifications made to glaciated grid cells included assigning (1) zero vegetation cover, (2) maximum and minimum albedo values consistent with the parameterization in the CMB model (Sect.2.2), and (3) a soil moisture availability of 1.0 (from an original value of 0.95).The conventional bulk computation of the latent heat flux in the WRF surface module is multiplied by the last parameter; therefore, this change was made for consistency with the CMB model.
The atmospheric model was forced with ERA-Interim data at 0.75 • × 0.75 • spatial-resolution and 6-hourly temporalresolution, as provided by the European Centre for Medium-Range Weather Forecasts (ECMWF; Dee et al., 2011).In the ERA-Interim reanalysis, snow depth is arbitrarily set to 10 m in the analysis for grid cells with greater than 50 % glacier coverage (Paul Berrisford, personal communication, 2012), which results in unphysical snow depths over the Karakoram.We therefore obtained the initial snow condition from the microwave-derived Global EASE-Grid 8-day Blended  (Brodzik et al., 2007), assuming a snow density of 300 kg m −2 and assigning an initial depth of 2 m over large glaciers where these data are missing (less than 0.1 % (8 in total) of data points in the region spanned by WRF D1).WRF employs a terrain-following hydrostatic-pressure coordinate in the vertical, defined as eta (η) levels (Skamarock and Klemp, 2008).For these simulations, the lowest atmospheric model level was specified at η = 0.997585 (∼ 20 m) to maintain the validity of the constant-flux assumption in the bulk computation of the turbulent heat fluxes, as the surface mid-layer height (less than 10 m) is used in the calculation following the approach of the Noah LSM.We selected the recently revised Monin-Obukhov surface layer (Jiménez et al., 2012), which was found to improve the simulation of the diurnal amplitudes of near-surface meteorological fields over complex terrain with a horizontal spatial resolution of 2 km.We also used positive-definite explicit 6th order diffusion (Knievel et al., 2007), in order to dampen grid-scale noise in the atmospheric fields and because Mölg and Kaser (2011) found this option improved the simulated magnitude of precipitation at high elevations on Kilimanjaro.For the simulations presented here, we selected the default value of the diffusion coefficient (0.12) for all model domains except D3, for which we used a value of 0.36.The choice of the diffusion parameter value is uncertain; sensitivity runs revealed that increasing the strength increased simulated precipitation at high elevations, which may be attributable to increased diffusive transport, with the best agreement with the Urdukas AWS data found for the selected value.

Surface energy and mass balance model
The CMB model is described fully by Mölg et al. (2008Mölg et al. ( , 2009) ) with the most recent updates in Mölg et al. (2012a), but we will review some important features here.The model computes the column specific mass balance from solid precipitation, surface deposition and sublimation, surface and subsurface melt, and refreeze of both meltwater and liquid precipitation.To determine the mass fluxes, the model first solves the surface energy balance equation: in which the terms correspond to, from left to right: incoming short-wave radiation, broadband albedo, incoming and outgoing long-wave radiation, turbulent fluxes of sensible and latent heat, ground heat flux and heat flux from precipitation.The ground heat flux, QG, consists of a conductive component (QC) as well as a component due to subsurface penetration of short-wave radiation (QPS).The net flux, F NET , represents the energy available for melt, QM, provided the surface temperature is at the melting point, T M = 273.15K.
The model treats both surface and subsurface processes, including surface albedo and roughness evolution based on snow depth and age; snowpack compaction and densification by refreeze; and the influence of penetrating solar radiation, refreeze and conduction on the englacial temperature distribution.The CMB model is forced by air temperature, humidity, wind speed, and air pressure, all of which were taken from the lowest model level (z = 20 m).Note that the diagnostically updated 2 and 10 m meteorological fields were not used as forcing so as to (1) be consistent with the approach of the Noah LSM (Chen and Dudhia, 2001), and (2) prevent decoupling of the atmosphere and land surface, wherein the lower atmosphere is no longer influenced by surface conditions.The CMB model also takes as input: total precipitation and its frozen fraction; incoming short-and long-wave radiation; and time between snowfall events.Some model parameter values are provided in Table 2.The initial subsurface temperature was specified through linear interpolation of the input data to the Noah LSM, available at 0.1, 0.4, 1.0, and 2.0 m depths, and assigning a constant value of 268.6 K below this level.The lower boundary is specified at 268.6 K during the simulation, based on measurements taken from a Tibetan glacier (Mölg et al., 2012a).We address uncertainties in the subsurface temperature initialization by including a long (25 day) model spin-up period.

Coupling architecture
For both offline and interactive simulations, the CMB model calculates glacier surface-energy and surface-and nearsurface-mass fluxes among other variables at every time step (e.g.every 4 s in D3) over glaciated grid cells, while the Noah LSM values are retained over non-glaciated grid cells.The integrated modelling approach permits some advantages in the CMB model's forcing strategy.For example, topographic shading, incoming short-and long-wave radiation, and the fraction of frozen precipitation are now obtained from the atmospheric model's surface, radiation, and microphysics modules, respectively.Another important advantage is that WRF provides high-resolution, dynamically derived, and spatially distributed forcing data without the need for traditional statistical methods, such as those mentioned in Sect. 1.The incorporation of the CMB model for all glacier grid points in the coupled model adds negligible computational expense to WRF simulations.For interactive simulations, the CMB model updates over glaciated areas in WRF, at every time step: (1) surface heat and moisture fluxes, (2) surface and subsurface (including deep soil) temperature, (3) snow depth, water equivalent and fractional cover, (4) surface albedo and roughness, and ( 5) surface specific humidity.The inclusion of feedbacks represents a more consistent approach, as it permits the nearsurface forcing variables to be modified by exchanges of mass, momentum and moisture between the glacier and the atmospheric surface layer.In this study, the CMB model output accumulated energy and mass fluxes every hour that were then converted into hourly averages for analysis; these data will be referred to as "hourly".
As indicated at the beginning of Sect.2, it is not absolutely correct to label the two forcing approaches as "offline" and "interactive" because the atmospheric model currently receives surface feedbacks through the Noah LSM.There have been recent efforts to improve the simulation of snow processes in WRF, such as with the introduction of the Noah-MP land surface parameterization (Niu et al., 2011), which, for example, introduces separate vegetation canopy and surface layers and the possibility of multiple vertical layers in the snowpack.However, the simplified treatment of glacier grid cells in the Noah LSM is retained.Thus, by incorporating the CMB model, we are able to simulate more physical processes relevant for glaciers, such as refreezing of meltwater in the snowpack, englacial melt, and formation of superimposed ice.Other improvements to the treatment of snow and ice physics, compared with the Noah LSM, include introducing multiple layers in the snowpack, increasing the column depth from 2 to 9 m, consideration of snow porosity, and allowing for full snowpack ablation to expose bare ice.The latter point is especially critical, as the Noah LSM imposes minimum snow depth and water equivalent values over land-ice grid cells.

Measurements for model evaluation
In Sect.3.2, we compare the coupled model results with a limited number of available ablation stake measurements as well as automated weather station (AWS) data that were acquired in summer 2004 on the Baltoro glacier (35 Mihalcea et al., 2006).The glacier is approximately 62 km long, with an average (maximum) width of 2.1 (3.1) km (Mayer et al., 2006).Therefore, in WRF-CMB the Baltoro glacier is represented by at least one grid point in the along-glacier direction and we resolve longitudinal rather than transverse gradients in surface conditions.We use data from 6 sections (SF, U, G1, C, BN, and BS) as well as from longitudinal transects along the glacier (L1, L2 and L3), comprising 53 stakes in total that provide sufficient spatial coverage (cf.Figs.1b or 4b) to evaluate both the spatial pattern and the magnitude of ablation in the coupled model applied to the Baltoro glacier.The ablation measurements were taken at different intervals between 1-15 July 2004; a brief summary of the location and other details of the stake measurements are given in Table 3 (a more detailed description of the data can be found in Mihalcea et al., 2006).While the data represent a brief period, they provide the only available direct ablation measurements in the Karakoram.For the comparison, total simulated surface lowering was interpolated to the mean location of the stake section or transect using inverse distance weighting.
The AWS was situated adjacent to the glacier on a moraine ridge at an elevation of 4022 m (35 • 43.684 N, 76 • 17.164 E) and provided hourly mean data after 18 June 2004 (Mihalcea et al., 2006).We compared these data with WRF data from the nearest model grid point (located at an elevation of 4322 m), which was also nonglaciated and therefore was more consistent in the land surface type.However, the data therefore do not include direct feedbacks from the CMB model.Note that the assumptions discussed in Sect.2.1 for snow initialization were not applied over the stake sites on the main glacier area (cf.Fig. 1b).
To supplement these field measurements, we also evaluate in Sect.3.1 the basin-scale performance of the coupled model E. Collier et al.: Interactive glacier-atmosphere climatic mass balance modelling using MODIS/Aqua (1) MYD10A1 daily snow albedo available at 500 m resolution, and (2) MYD11A2 eight-day land surface temperature available at 1 km resolution, with daily data obtained by averaging day-and night-time temperatures where both fields were available and were assigned the highest quality assurance flag for MODIS products.Due to the prevalence of missing data in the snow albedo dataset, we considered only the grid cells with at least 25 % valid observations during the 67 day period for comparison with WRF.Both MODIS datasets were re-projected to the WRF D3 grid before completing the analysis.

Results and discussion
We first compare our simulated results with remote sensing data (Sect.3.1) and with meteorological and glaciological measurements from the Baltoro glacier (Sect.3.2).The role of interactive coupling on the atmospheric forcing data and on simulations of CMB will then be discussed.Results are presented from the finest-resolution atmospheric model domain only, since it provides the most realistic terrain representation.

Remote sensing data
Figure 2 presents a comparison between WRF-CMB and the MODIS/Aqua datasets discussed in Sect.2.4.The elevational profile of land surface temperature (LST) averaged over the simulation period produced by the CMB model is in good agreement with the MODIS data above ∼ 5200 m and is an improvement on the Noah LSM values at all resolved elevations (Fig. 2a).The strong divergence of modelled and observed LST below 5200 m likely results from neglecting debris cover, since its presence allows the glacier surface to be warmed by solar radiation above the melting point.Supraglacial debris extent has also been found to increase with distance down glacier in remote-sensing case studies of the central Karakoram (e.g.Scherler et al., 2011).Specific to the Baltoro glacier and its tributaries, Mayer et al. (2006) found that debris coverage increased to 70-90 % of the glacier area below 5000 m, with 100 % coverage found below the U site (Fig. 1b).A time-series analysis of LST, averaged only over elevations greater than 5200 m is presented in Fig. 2c.The CMB model gives an improved performance over the LSM alone, although LST is generally under-predicted, with mean biases of −1.0, −1.3, and −6.1 K in the INT, OFF, and Noah LSM simulations, respectively.Conversely, snow albedo in WRF-CMB is in good agreement with MODIS below ∼ 5200 m (Fig. 2b), although simulated values are constrained by the lower bound of α ice = 0.3 as snow depth goes to zero and, thus, slightly overestimate the observational values at the lowest elevations.Above 5200 m, WRF-CMB over-predicts snow albedo compared with MODIS.However, it produces values and an altitudinal gradient that are in much better agreement with observations than the Noah LSM.
The strong discrepancy between Noah LSM and MODIS data is in part related to the treatment of grid cells defined as glacial ice: the LSM in WRF v3.4 imposes minimum values of snow depth and water equivalent of 0.5 and 0.1 m, respectively, thus preventing the exposure of bare ice or debris and the associated lowering of surface albedo.In addition, the Noah LSM employs a time-decaying snow albedo formulation (based on the scheme of Livneh et al., 2009) and determines surface albedo using fractional snow cover to correct a background snow-free albedo.Although snow albedo is likewise an exponential function of age in the CMB model (following Oerlemans and Knap, 1998), the actual surface albedo also depends on snow depth to account for surface darkening when the snowpack is thin.It is clear from Fig. 2 that this formulation, in combination with permitting snowfree conditions, gives more realistic values.
The evaluation of modelled albedo is sensitive to the simulated timing of snowfall events for both models, due to the nature of the parameterization schemes, and is limited by the large number of missing data in MODIS.In addition, comparison is hindered by the fact that the MODIS daily albedo product is not a daily-averaged quantity (as the simulated data are) but rather a collection of pixels with the highest quality from that day for which acquisition time is not trivial to retrieve (Stroeve et al., 2006).Finally, the small positive biases in snow albedo below ∼ 4300 m and the negative biases in LST are both physically consistent with neglecting the influence of debris cover in WRF-CMB, as will be discussed further in Sect.

Baltoro glacier
Figure 3 presents a time series of modelled and observed near-surface meteorological data from the Urdukas AWS that is situated adjacent to the Baltoro glacier.WRF-CMB is skillful in simulating air temperature at 2 m, and its evolution over the study period, including capturing periods of reduced diurnal variability at the beginning and between 30 July and 6 August.However, the good agreement in near-surface temperature despite a difference in real and modelled elevation of ∼ 300 m (4022 vs. 4322 m a.s.l., respectively) suggests that there is a positive temperature bias in WRF-CMB at this grid point.The greatest contributing factor is higher incoming short-wave radiation: averaged over the simulation period, the surface in INT receives an additional 112 W m −2 more radiation than measured by the AWS (not shown).The discrepancy is most likely due to insufficient simulated cloud cover and humidity (Fig. 3b), with a potential contribution also from the computation of topographic shading at 2. cycle is smaller in WRF-CMB, which may be attributable to differing thermal properties of the real and modelled land surface or to the fact that the AWS sensor was not aspirated.
The magnitude of the near-surface wind speed is also in agreement with the AWS data (Fig. 3c).However, an important discrepancy is the underestimation of precipitation at this particular grid cell in both INT and OFF simulations: the AWS records a total of 122.8 mm of precipitation between 25 June and 31 August, while INT and OFF simulate 46.9 and 48.4 mm, respectively (Fig. 3e).Missing precipitation events are also reflected as discrepancies in the time series of relative humidity (cf.Fig. 3b, e) and are consistent with an overestimation of incoming short-wave radiation as a result of too little cloud cover.The disagreement in measured and simulated humidity and precipitation may reflect several sources of error, such as in the forcing data at the lateral boundaries.In addition, the spatial resolution of WRF D3 may be insufficiently fine to fully resolve orographic uplift or microscale complex flow features that affect precipitation at the AWS.Furthermore, we do not use a cumulus parameterization in the finest model domain and therefore assume that convection is explicitly resolved.However, previous studies indicate that a grid spacing on the order of 100 m (Bryan et al., 2003;Petch, 2006) or even 10 m (Craig and Dörnback, 2008) is needed to capture the dominant length scales of moist cumulus convection.A final potential error source is the difference in the land surface type adjacent to the AWS and model grid cell: the Baltoro glacier is debriscovered at the Urdukas site, while WRF-CMB has a clean snow/ice surface.The differing thermal properties of the adjacent surface area, specifically the limiting of temperature at the melting point in WRF-CMB, may also contribute to differences in localized convection.
WRF-CMB produces 40 to 80 cm of ablation along the main body of the Baltoro glacier between 1-15 July 2004 (Fig. 4a).Spatial comparison of the two simulations reveals only small differences, generally on the order of a few centimetres, consistent with the short nature of the study period (Fig. 4b).There are slightly positive anomalies at lower elevations, corresponding to less ablation in INT; conversely, there are negative anomalies at higher elevations, corresponding to more ablation in INT.Total simulated ablation is in order-of-magnitude agreement with measurements (Fig. 4c); however, the model overestimates ablation at all sites, in part because it does not capture four all-phase precipitation events, amounting to 17.6 mm, during the measurement period (cf.Fig. 3e).In comparing daily simulated/measured ablation rates and mean debris thickness (Table 4), the rates tend to be in better agreement for sites with thinner mean debris cover (SF, L2, L3, C) and more strongly overestimated by WRF-CMB where supraglacial debris is thicker (U, L1, G1, BN, BS).Although differences between INT and OFF are small, INT is in closer agreement with observations at all but one site (BN, as a result of less simulated refreeze than in OFF), with the strongest improvement at BS.The improvement at this site stems from faster complete snow cover removal (∼ 1 day earlier in INT), which reduces subsurface penetration of short-wave radiation and, thus, subsurface melt production.Finally, the overestimation of ablation by WRF-CMB tends to diminish as the observation period increases (Fig. 4d), which then suggests that the coupled model as configured in this study may be best suited for "climatological" simulations of glacier mass balance due to its sensitivity to the timing of precipitation.Mihalcea et al. (2008) performed distributed surfaceenergy sub-debris melt modelling, using the Urdukas AWS data as forcing for the same study period.The authors determined debris extent, thickness and thermal properties from satellite imagery, and considered only the elevation range of 3650-5400 m a.s.l., which gives a corresponding glacier area of 124 km 2 .Mihalcea et al. (2008) computed 0.058 km 3 w.e. of ablation, or a mean surface lowering of 0.47 m between 1-15 July, which was found to be a slight underestimation (on average, −0.016 m) of the observed ablation rates at the SF site.For comparison, INT and OFF simulate 0.069 and 0.070 km 3 w.e. of surface melt, respectively, over an area of 126 km 2 that produces an average thickness change of approximately −0.55 m.The actual CMB calculation in the model also includes additional processes, such as snowpack ablation and surface vapour fluxes, that bring the total simulated mass loss to 0.078 km 3 w.e. between 1-15 July.We employ the same glacier outline, that of Mayer et al. (2006); however, discrepancies in our estimates may arise from its projection to the WRF D3 grid, differences in removal of tributary glaciers, and the coarser representation of the Baltoro glacier at 2.2 km spatial resolution (vs.90 m in Mihalcea et al., 2008).Despite these and other sources of disagreement, comparing the two estimates gives an approximate measure of the effect of neglecting debris, which is thought to cover 38 % of the Baltoro glacier (Mayer et al., 2006) and 73 % of the altitude range of the main glacier tongue considered in Mihalcea et al. (2008), in our simulations.

Influence of interactive coupling
Figure 5 presents a time series of daily means of the nearsurface WRF meteorological data used as forcing for the CMB model and provides the context for the fluctuations of surface energy and mass fluxes discussed in this section.Near-surface air temperatures in INT are higher by 0.3 • C on average than in OFF (Fig. 5a).The difference arises primarily from a reduced amplitude of the diurnal cycle, with higher nocturnal temperatures (Fig. 5a subpanel).INT simulates higher surface temperatures (T sfc ), as well as higher subsurface temperatures in the top 0.5-1 m (peak differences are ∼ 0.7 • C, not shown), as a result of stronger downward long-wave radiative forcing (see Fig. 5f for daily average curves).The increase in L↓ is expressed between evening and early morning and is a direct result of higher mixing ratios at 2 m in INT (not shown).The change in radiative forcing in INT translates into less heat extraction from the surface layer, through a reduced nocturnal QS, and, in turn, into the near-surface temperature difference.Note that the near-surface air temperature evolution simulated in INT may represent an improvement, as Mölg et al. (2012b) found that WRF + Noah LSM can produce an excessively large diurnal cycle as a result of a nighttime cold bias at 2 m compared with AWS measurements on Kilimanjaro.Interactive coupling also results in a reduction of mean incoming short-wave radiation (−9.0 W m −2 ) and, as previously mentioned, a mean increase in incoming long-wave radiation (2.4), changes that arise from alterations to atmospheric clouds and moisture (see Fig. 6).Basin-scale daily-mean differences in the other forcing variables for the CMB model are negligible.
The atmospheric changes induced by including feedbacks from the CMB model are generally small in magnitude and limited in vertical extent, but still appreciable.Air temperature and mixing ratio anomalies are generally confined to the lowest 10 model levels, which correspond to the layer between a mean surface pressure of 543 hPa and the level of 450 hPa (Fig. 6a).Vertical changes in the mean cloud cover fraction are variable, with the greatest differences present near the levels of 375 (η = 12-14) and 125 hPa (η = 26-29; Fig. 6d).However, interactive coupling has a strong warming influence on the subsurface temperature distribution (on average, +2.6 K; Fig. 6c), as a result of (1) the inclusion of the energy flux from penetrating solar radiation, and (2) the method for updating deep soil temperature (T ds ), which is defined at a depth of 3 m.With regard to the latter point, T ds in INT is taken from the CMB model subsurface scheme, which resolves the column to a depth of 9 m but is constrained by a lower boundary temperature of 268.6 K in this study.In contrast, the Noah LSM updates T ds using a weighted combination of the annual mean T sfc of the previous year and of the last 150 days as the data become available, with no lower threshold imposed.The resulting minimum values for T sfc in the CMB model and Noah LSM are ∼ 245 and 224 K, respectively.
The non-negligible influence of interactive coupling on the near-surface meteorological forcing data translates primarily into reduced ablation of snow and ice in INT (Fig. 7a and b).Area-averaged modelled surface height lowering is smaller and total mass balance is less negative in INT, with a mean reduction in ablation over the Karakoram basin of 0.1 m w.e.(−6.0 %), to a cumulative value of −1.5 m w.e. by 31 August.The difference in the total mass balance arises despite higher T sfc in INT (Fig. 7c).The inclusion of additional processes, such as the refreezing of meltwater, and the different method of subsurface temperature calculation both contribute to higher T sfc in both INT and OFF compared with the Noah LSM.distribution is characterized by a shallowing of the VBP above ∼ 5000 m, associated with (1) an increase in the positive vertical gradient of the fraction of solid precipitation that contributes positively to CMB (Fig. 8c), and (2) cooling of mean surface temperature with height to below the melting point (cf.Fig. 2a).Above 5875 m, the VBP profile again steepens as a result of large increases in accumulated, solid precipitation.Below 5875 m, INT produces less ablation (on average, 117.7 mm w.e.), while above this level it simulates a mean increase in accumulation (13.4) in part due to small increases in both accumulated precipitation and its frozen fraction (cf.Fig. 8b and c).Averaged over the whole period, the equilibrium line altitudes are 5469 and 5536 m in INT and OFF, respectively, which exceed the annual and generalized estimate of 4500 m by Hewitt (2005) and of 4200-4800 m by Young and Hewitt (1993), because we only simulate the ablation season.
Figure 9 presents the surface fluxes of energy and mass from the interactive simulation.On average, the main energy sources are incoming radiation, S↓ (374.3W m −2 ) and L↓ (220.4), with smaller contributions from QS (9.5) and QC (8.1; Fig. 9a).The main energy sinks are outgoing L↑ (−306.1)and reflected S↑ (−186.4),followed by QM (−74.0),QPS (−32.1), and QL (−13.9).Mass gains are, in general, dominated by refreeze (1.0 kg m −2 ) and solid precipitation (0.9; Fig. 9b), while mass loss is primarily through surface (−19.1) and subsurface (−3.2) melt.the snowfall event that occurs at the beginning of August (cf.solid precipitation bars in Fig. 9b) is clearly associated with (1) a reduction in S↓ and an increase in L↓, (2) a spike in both surface albedo and thus S↑, and (3) a reduction in absorbed short-wave radiation that translates into reduced energy for surface and subsurface melt.Furthermore, changes in glacier surface conditions have a noticeable feedback on the atmosphere during and after the snowfall event (cf.e.g.S↓ in Fig. 5e or cloud fraction changes in Fig. 6d).
Interactive coupling has the strongest influence on the net short-wave and ground heat fluxes in the atmospheric model (Fig. 10a).The average QG in INT (−23.7 W m −2 ) greatly exceeds that simulated by the Noah LSM alone (−0.5), due to (1) the inclusion of penetrating short-wave radiation, which always represents an energy sink at the surface, and (2) higher surface temperatures, which result in a stronger (more negative) flux downward to the subsurface.Mean absorbed short-wave radiation is also much larger in INT (184.8 vs. 107.9),as a result of lower average surface albedo (0.49 vs. 0.71; see Fig. 2c) over the simulation period.Smaller changes in the turbulent heat fluxes reflect in part different treatments of surface roughness, which is a spatially and temporally varying parameter in WRF-CMB that ranges between 0.8 and 2.6 mm as a function of snow age and generally exceeds the constant value of 1 mm specified by the Noah LSM for snow/ice surfaces.
In the CMB model, interactive coupling induces the largest magnitude change in the net short-wave (−3.3 %,) and longwave (+1.8 %,) radiative fluxes, as a result of the changes to S↓ and L↓ in the atmospheric model discussed previously.The CMB model fluxes of QL and QPS decrease (become less negative) on average in INT by 8.9 % and 3.5 %, respectively (Fig. 10b).In addition, QC and QS are reduced by 9.0 % and 5.0 %, respectively.The reduction in QPS is associated with less snow-cover-free glacier area in INT, which reduces the subsurface absorption of short-wave radiation  (Mölg et al., 2008).The reduction in the turbulent heat fluxes appears to stem from reduced 10 m wind speeds (−0.2 to 4.3 m s −1 ), due to the surface roughness changes discussed above.The changes also occur despite (1) a weaker correction for atmospheric stability on average, according to a modified version of the Monin-Obukhov stability function (Eq.12 in Braithwaite, 1995); and (2) stronger mean gradients in vapor pressure ( VP; −1.2 in INT vs. −0.7 hPa in OFF) and temperature ( T ; 1.8 vs. 1.2 K) between the nearsurface and lower boundary.
The net result is a decrease of 4.8 % in the average energy available for surface melt, QM, and less negative mass balance over the INT simulation (cf.Fig. 7b).The difference is primarily reflected in a reduction in surface melt (−5.8 %; Fig. 10c), and is compounded by increases in both refreezing in the snowpack and the formation of superimposed ice (8.5 % combined) as well as by greater solid precipitation (4.0 %).In general, mass exchanges between the glacier surface and the overlying boundary layer are smaller in INT, with the weakening in QL resulting in less sublimation (particularly at night) and deposition (at all times).

Remarks and perspectives for future research
The explicit approach to modelling alpine glacier climatic mass balance using WRF first demonstrated by Mölg and Kaser (2011) has been applied in offline mode for simulations of small glaciers (Mölg and Kaser, 2011;Mölg et al., 2012a,b) and has yielded important insights into the physical processes and the atmospheric forcing underlying mass fluctuations.However, here we demonstrate that feedbacks due to CMB processes are evident for the heavily-glaciated Karakoram region.The basin-scale influence of interactive coupling on the atmospheric forcing data, while moderate, acts to reduce the energy available for surface melt and, in concert with both reduced mass exchanges between the surface and boundary layer and increased refreezing, reduces modelled ablation during the summer of 2004.Furthermore, we demonstrate that the inclusion of additional real processes such as CMB feedbacks renders WRF-CMB capable of simulating observed magnitudes of CMB.
To the best of our knowledge, only one previous study, that of Kotlarski et al. (2010b), has performed interactive and distributed simulations of alpine glacier mass balance, achieved by introducing a subgrid-scale parameterization for glaciers and their areal changes into a regional climate model configured with a relatively coarse spatial resolution of ∼ 18 km.However, their implicit treatment "pools" the glaciers located in a grid cell into one ice mass at a fixed altitude and with a uniform snow depth.Furthermore, they quantify the effect of two-way coupling by comparing their interactive simulation with a control run that contains no glaciers (i.e.snow and ice surfaces are compared with bare soil or vegetationcovered surfaces), thus obscuring the exact role of feedbacks.Therefore, this paper presents the first assessment of the importance and strength of interactions between alpine glaciers and the atmosphere on explicit simulations of CMB.
We have shown results for only a short study period, of one ablation season, to evaluate the novel approach and its performance against the available measurements.The Karakoram, and the Himalayan region in general, are very data-sparse (e.g.Bolch et al., 2012), due to the expense and logistics of field surveys in this remote region.The model captures the magnitude of the 53 stake measurements of Mihalcea et al. (2006) the greatest test for some of the simplifying assumptions employed, such as zero debris cover.However, a longer application is needed to assess the year-long and inter-annual influence of interactive coupling, as well as the long-term performance of the atmospheric model under the climatesimulation forcing strategy we employ (i.e. with no nudging or model re-initialization; e.g.Maussion et al., 2011).Simulations of glacier mass balance are also inherently sensitive to the modeled solid precipitation (Mölg and Kaser, 2011), which is influenced in our study by the choice of microphysics scheme.Furthermore, the optimal choice of diffusion scheme, its strength, and its influence on simulated precipitation and therefore glacier CMB are beyond the scope of this paper and have not been investigated fully for our area of interest and model configuration.The simulation of near-surface meteorological fields by WRF over glacier surfaces has been found to be relatively insensitive to the choice of physical parameterizations (Claremar et al., 2012); however, the extent to which modelled CMB is dependent on the model physics, the choice of numerics, and the spatial resolution of the finest domain represents an important uncertainty that will be explored in a future study.
The mean proportion of debris covered-area on Karakoram glaciers is estimated to be 18-22 % (Scherler et al., 2011;Hewitt, 2011), which is higher than the pan-Himalayan average of ∼ 10 % (Bolch et al., 2012).Specific to the Baltoro glacier, Mayer et al. (2006) estimate that ∼ 38 % of the total glacier area is debris covered.The presence of debris above a threshold, or "critical thickness", of ∼ 2 cm has been shown, empirically and through surface energy balance modelling, to reduce glacier ablation as a result of its insulating effect (e.g.Østrem, 1959;Kayastha et al., 2000;Nicholson and Benn, 2006;Reid et al., 2012).The range of mean debris thicknesses at the stake sites is 2.0-18.0cm (Table 4), suggesting that on the whole insulation effects should dominate over the lowering of surface albedo except at the sites L2 and L3 where debris thickness is approximately equal to the critical value.Indeed modelled ablation closely matches the measured rate at these two sites and elsewhere is overestimated by WRF-CMB (Fig. 4c), physically consistent with the exclusion of debris in this study.This interpretation is supported by the first distributed ablation modelling study for debris-covered ice, that of Reid et al. (2012), which found reduced sub-debris ablation when depth exceeded 2 cm.However, it is noteworthy that geodetic estimates of early 21st century elevation changes in the Karakoram (Gardelle et al., 2012;Kääb et al., 2012) do not show a difference between clean and debris-covered ice.
Given the similarity of the underlying surface types in INT (snow/ice) and OFF (snow) influencing the atmospheric forcing data, the difference in simulated CMB for the clean glacier simulations is relatively small.From the results presented here, it could be expected that the inclusion of feedbacks is not essential for small glaciers or less glaciated basins.However, we would expect the interactive inclusion of the CMB model to have a larger influence for glaciers with significant debris cover, as its presence alters surface temperature and moisture properties and thus turbulent exchanges with the surface boundary layer (e.g.Takeuchi et al., 2000).To assess the role of feedbacks for debris-covered glaciers and to allow the WRF-CMB modelling system to provide long-term, accurate simulations in the Karakoram, including the effects of debris cover on surface conditions and glacier ablation represents important future work.Treating debris cover in distributed mass balance modelling is also becoming more important in light of observations of increasing debris-covered area in many regions (e.g.Stokes et al., 2007;Bhambri et al., 2011).Another process that is thought to be important for Karakoram glaciers is accumulation via snow and ice avalanching (e.g.Hewitt, 2011), which may be useful to parameterize.Finally, dynamical ice flow changes have been shown to be important when quantifying the response of Himalayan glaciers to climate fluctuations on multiannual timescales (e.g.Scherler et al., 2011;Gardelle et al., 2012;Kääb et al., 2012;Azam et al., 2012).

Conclusion
CMB feedbacks have been introduced into a new, multi-scale modelling approach for explicitly resolving the surface and climatic mass balance processes of alpine glaciers, and this technique has been extended to the regional scale.Although validation data is sparse, the model captures the magnitude of available in situ measurements, with improvements arising from including feedbacks from the CMB model to WRF.Furthermore, discrepancies between observed and simulated ablation can be attributed to physical processes neglected as simplifying assumptions, particularly debris cover effects.
Both components of WRF-CMB are based on physical principles, with no statistical downscaling at their interface.The direct linkage increases the applicability of this approach for the simulation of the past-and future-climate response of glaciers, since the modelling system produces a physicallyconsistent response to changes in external forcing.Incorporation of the CMB model also increases the number of physical processes important for glaciers represented in the atmospheric model, and provides a consistent calculation of surface energy and mass fluxes, since changes in glacier surface conditions are permitted to influence the atmospheric drivers.Perhaps the most important advantage, however, is that WRF-CMB permits direct causal attribution of glacier mass changes to both physical processes and the main atmospheric drivers.With further development, the model has the potential to bridge the data gap in the Karakoram and shed light on the role of climate forcing in the anomalous behaviour of glaciers in this region.
The offline approach to simulations of CMB, as well as the simplified representation of glaciers in regional atmospheric models, essentially treats either the atmosphere or alpine glaciers simply as a boundary condition.We suggest that this unified, explicit approach should be increasingly adopted in future studies, particularly for heavily glaciated regions.

Fig. 1 .
Fig. 1.(a) WRF atmospheric model domains, configured with horizontal spatial resolutions of 33, 11, and 2.2 km.Terrain elevation from the GTOPO30 dataset is shaded in units of meters.(b) Outline of Baltoro glacier and its tributaries, which are included in WRF D3, with the mean stake locations labeled and denoted by stars.

Fig. 2 .
Fig. 2. Comparison between WRF-CMB D3, Noah LSM, and MODIS for land surface temperature, (a) averaged from 25 June-28 August 2004, and in 50 m elevation bins; (c) the mean time series above 5200 m elevation; and for (b) snow albedo, averaged between 25 June-31 August 2004, over glacier grid cells where at least 25 % of the daily times are available.

Fig. 3 .
Fig. 3. Hourly (a) air temperature at 2 m, (b) relative humidity at 2 m, (c) wind speed at 10 m, and (d) surface air pressure, as well as (e) daily total precipitation.Solid black (blue) curves display data from the interactive (offline) simulations while the dashed grey curve is the Urdukas AWS station data.Note the difference in elevation of the AWS (4022 m a.s.l.) and the terrain height in the closest WRF grid cell (4322 m).
2 km resolution.Finally, the amplitude of the diurnal temperature www.the-cryosphere.net/7Ablationrates (cm day −1 ) and debris thickness on the Bal- Fig. 4. (a) Total and (b) INT-OFF surface height change between 1-15 July 2004, in the vicinity of the main tongue of the Baltoro glacier.Stake site and transect locations shown in (b), with additional information provided in Table 3. White grid cells correspond to non-glaciated area.(c) Measured mean ablation (triangles) at stake locations, with range of observed values denoted by bars.Simulated INT (OFF) ablation shown by black squares (blue circles).

Fig. 5 .
Fig. 5. Daily mean (a) air temperature, (b) relative humidity, (c) wind speed, (d) total precipitation, incoming (e) normal short-wave and (f) downward long-wave radiation at ground surface, (g) air pressure, and (h) frozen fraction of precipitation, area-averaged over all glaciated grid cells.Data for (a-c), and (g) are taken from the lowest model level (z = 20 m).The subpanel in (a) presents the average diurnal temperature cycle over the simulation period.Black (blue) curves display data from the interactive (offline) simulation.

Fig. 6 .
Fig. 6.Vertical and subsurface distribution of the influence of interactive coupling over glaciated areas illustrated by hourly time series of the change (INT-OFF) in area-averaged (a) air temperature,(b) specific humidity, (c) subsurface temperature, and (d) cloud fraction.

Fig. 7 .
Fig. 7. Daily basin-scale averages of (a) accumulated surface height change, (b) accumulated total mass balance, and (c) surface temperature.Black (blue) curves display data from the interactive (offline) simulation.For reference, surface temperature simulated by the Noah LSM is the dashed grey curve in (c).
Fig. 8. (a) The vertical mass balance profile of the Karakoram basin at the end of the simulation.The altitudinal dependence of (b) total accumulated precipitation, and (c) mean frozen fraction, averaged over the simulation.Data are area-averaged in 50 m elevation bins.

Fig. 9 .
Fig. 9. From the interactive WRF-CMB simulation: daily (a) mean surface energy balance components (left y-axis; see Eq. 1 for explanation of symbols) and albedo values (grey right y-axis), and (b) sums of mass fluxes.The radiation variables are shown in (a) as solid (directed downward) and dashed (upward) lines, albedo as grey dots, and the other surface energy fluxes as bars.The heat flux from precipitation (QPRC) is negligible and not shown.Values are averaged over glaciated grid cells only.

Fig. 10 .
Fig. 10.Area-averaged mean difference (INT − OFF) over the simulation period and over glaciated grid cells of (a) the main components of the WRF surface energy budget, (b) the MB model energy fluxes, and (c) the MB model mass fluxes.Symbols in (a) represent, from left to right, net short-and long-wave radiation, ground heat flux, and turbulent fluxes of sensible and latent heat.Note that the sign convention for the turbulent fluxes in (a) is opposite to (b).Symbols in (b) are discussed in Sect.2.2 , and simulating the ablation season likely represents www.the-cryosphere.net/7

Table 3 .
Summary of available ablation stake measurements from the Baltoro glacier.