Interactive comment on “Observed high-altitude warming and snow cover retreat over Tibet and the Himalayas enhanced by black carbon aerosols”

role of BC increase in the tropospheric warming of the Tibetan Plateau region using a coupled AOGCM at high horizontal resolution, by observationally-based BC aerosol of High resolution coupled OAGCM with coupling between snow over land and BC deposition. this last is crucial. This analysis is also based on a work where the BC forcing is re-estimated by using satellite + ground based optical depth, with a new methodology to separate the BC contribution to solar absorption from other aerosols and its direct rad forcing


Introduction
Himalayan mountain glaciers and snowpacks have a major impact on the water systems of major rivers throughout Asia and the people living in the river basins.Recent observations suggest a continuing decline in Himalayan mountain glaciers and snow cover.Bajracharya et al. (2008) observed that the Himalayan glaciers are retreating at rates ranging from 10 to 60 m per year, and many small glaciers have disappeared.Gardner et al. (2013) also showed with satellite observations the steady reduction of western China glaciers with the most rapid decline observed in the Himalayan mountain regions.Changes in the cryosphere are accompanied by documented surface warming trends over Tibet, which reveals a strong altitude dependence of surface warming with peak warming trends of 2-2.5 • C at 5000 m from 1961 to 2006 (Liu et al., 2009).
The last few decades also witnessed rapid growth in human population and economic activities, causing intense air pollution over the Asian region.Among the many air pollutants, black carbon (BC) aerosols have been shown to have a significant impact on global and regional climate change (Ramanathan and Carmichael, 2008).Many previous studies have linked Asian aerosols (including sulfates and BC) with monsoon systems and have demonstrated the aerosol impact on the summer rainfall (Ramanathan et al., 2005;Lau et al., 2006;Lau and Kim, 2006;Meehl et al., 2008).The BC aerosols have also been shown to have an impact on warming trends over the Himalaya-Tibet region (Ramanathan et al., 2007), on the retreat of Himalayan glaciers (Menon et al., 2010;Qian et al., 2011), and on Eurasian snow cover (Flanner et al., 2009).Observationally, using ice-core samples to reconstruct historical BC content over Tibet, Xu et al. (2009) suggested that BC is a significant contributing factor in causing the glacier change.
To date global climate models forced by historical radiative forcing scenarios (such as those in the Coupled Model Y. Xu et al.: Observed high-altitude warming and snow cover retreat due to black carbon Intercomparison Project Phase 5, CMIP5) have difficulty in simulating the observed record surface warming (You et al., 2015) or its anomalously strong altitude dependence in the Tibet-Himalaya region.One possible explanation as we will investigate here is that few of these earlier studies of the Himalayan and Tibetan climate change have considered the combined effects of all the following factors: BC direct heating of the atmosphere, the heating of snowpacks and glaciers by BC darkening the snow and ice, the greenhouse effect of CO 2 , and the surface cooling effects by aerosols other than BC.
In this study, we used a state-of-the-art global climate model to conduct a suite of model experiments to understand BC's role in the cryosphere change over the Himalaya and Tibetan region.A unique feature of the present study, compared with earlier studies, is that BC radiative forcing is constrained with multiple sources of observations (satellite observed aerosol optical depths and ground network of spectral sun photometer measurements).We also used a newly developed method to separate the BC contribution to solar absorption from other aerosols (sulfates, organics, and brown carbon) and calculate its direct radiative forcing (Bahadur et al., 2012;Xu et al., 2013).Previous studies (Ramanathan et al., 2007;Lau et al., 2010;Menon et al., 2010;Qian et al., 2011) included the effects of BC on the atmosphere and the cryosphere, but the simulated BC radiative forcing by these "standard" models used in CMIP5 is strongly biased to low values (Bond et al., 2013) due to emission inventory biases and missing physical treatments (Jacobson, 2012).As shown in Bond et al. (2013), current models are underestimating BC solar absorption over South Asia by a factor of 2 to 5. In this study, we scaled the simulated BC forcing in the climate model by factors ranging from 2 to 4 to bring the simulated values closer to the observationally constrained values (Xu, 2014).Another improvement in this study is that the simulations were conducted using a fully coupled oceanatmosphere-land model at a high resolution of 1 • by 1 • , in which a new land snow module is adopted (Lawrence et al., 2011) to account for BC deposition effect on snow and ice.

The global climate model
CESM1 (Community Earth System Model 1) is a coupled ocean-atmosphere-land-sea-ice model.CESM1 climate simulations have been documented extensively (Meehl et al., 2013).The CESM1 (CAM5) used in this study is a version with a finite-volume nominal 1 • horizontal resolution (0.9 • by 1.25 • ) and 30-level vertical resolution.The highest model level is about 36 km (4 hPa) in the stratosphere, and lower levels close to surface (boundary layers) have vertical resolutions of about 100-200 m.CESM1 (CAM5) includes forcings from greenhouse gases (GHGs) as well as concentrations of tropospheric ozone and stratospheric ozone (Lamarque et al., 2010).The concentrations of various gases were calculated offline and prescribed into model simulations, unlike the aerosol loading calculated online from the emissions.The three-mode modal aerosol scheme (MAM3) has been implemented (Liu et al., 2012) and provides internally mixed representations of number concentrations and mass for Aitken, accumulation, and coarse modes of various aerosol species (sulfates (SO 4 ), BC, organic carbons, dust, sea salt).The new cloud microphysics scheme (Morrison and Gettelman, 2008) allows the number concentration of cloud drops and ice crystals to be affected by aerosol concentrations and therefore accounts for the "indirect radiative forcing" of aerosols.
The land model (Community Land Model, CLM4) also includes major updates, making it more versatile in simulating snowpacks (Lawrence et al., 2011).The sub-grid processes including melting, metamorphism, deposition, and redistribution are considered in a snow cover fraction parameterization (Niu and Yang, 2007).Other parameterizations include snow compaction (Lawrence and Slater, 2010) and the albedo calculations for snow on or around vegetation (Wang and Zeng, 2009).Compared to the previous model versions, the albedo contrast between snow-covered and non-snowcovered area is more consistent with observations.

BC effects on surface albedo
The deposition of BC particles, due to gravity or rainfall removal, is a mechanism to remove aerosols from the atmosphere, and therefore a sink term for the atmospheric BC mass balance.BC particles deposited onto the surface of high-albedo snow or ice would reduce surface albedo.The snow model of CLM4 is significantly modified via the incorporation of SNICAR (snow and ice aerosol radiation) module, which represents the effect of aerosol deposition (BC, organic carbon and dust) on albedo, introduces a grain-sizedependent snow-aging parameterization, and permits vertically resolved snowpack heating (Flanner et al., 2007).This new module considers the albedo change by counting the surface concentration of BC and it calculates the surface radiative energy flux at multiple wavelengths.The surface albedo change will consequently alter the energy balance at the surface and in the atmosphere.

BC atmospheric radiative forcing
The present-day BC emission is adjusted from the standard model emission inventory (Lamarque et al., 2010) to account for the potential model underestimation of BC forcing.Emissions over East Asia regions are increased by a factor of 2 and South Asia regions by 4. The emissions are adjusted by the same ratio in all economic sectors (energy, industrial, etc.) and all seasons by the same ratio.Our previous analysis showed that such a correction would improve modelsimulated radiative forcing compared with direct observations (Xu et al., 2013;Xu, 2014).Without the observationally constrained values, the modeled forcing and simulated temperature change would be lower by about a factor of 2 to 4.

Model experiments
To isolate the climate impact of individual forcing agents, we contrasted the perturbed model simulations with present-day forcing (Sect.2.3.2) to the long-term preindustrial control simulations (Sect.2.3.1).The approach is similar to the classical instantaneous CO 2 doubling experiment (Manabe and Wetherald, 1975).Additionally we conducted a fixed-SST (sea surface temperature) experiment for radiative forcing diagnostics (Sect.2.3.3) and the 20th century transient runs to better attribute the observational changes (Sect.2.3.4).The details of these simulations are given below.

Control simulation for preindustrial climate
We have a 319-year-long preindustrial control run, and extended it with an additional 75-year run to test whether there was any discernible drift in the mean climate state.The Northern Hemisphere temperature does not show any statistically significant drift.Therefore, we lay the foundation for our analysis by employing the original 319-year run and the extended 75-year run (394 years in total) as a control case.Natural variability of the climate system can be examined from the unforced 394-year preindustrial simulations.

Four sets of perturbed simulations with
instantaneously imposed present-day forcing: BC, SO 4 , CO 2 , and all three forcings combined The forcings were imposed by instantaneously increasing the emissions of BC, or the emission of SO 4 's precursor, sulfur dioxide, or by increasing CO 2 concentration to presentday level (400 ppm).Except for the adjusted BC present-day emission as detailed in Sect.2.2.2, all other emissions are from the standard inventory adopted by CMIP5 models, as described in Lamarque et al. (2010).We run the perturbed simulations in fully coupled mode for 75 years, starting from the end of the 319th year of the control simulation.The difference between the last 60 years (allowing the first 15 years for model spin-up) and the long-term control simulation provides the response signal due to the imposed forcing.With the concern that BC signal is potentially small compared with natural variability, five ensemble members of BC forced simulations are conducted to increase signal-to-noise ratio.Each model year costs about 2000 processor hours in a highperformance computing system.

Three sets of perturbed simulations but with fixed SST
These are also forced by the instantaneous increase in BC, SO 4 , and CO 2 , separately, but the model runs in atmosphere and land only mode with SST fixed at preindustrial level.These simulations are used only for diagnosing the radiative forcing due to various species.

The 20th century transient single-forcing simulations
The simulations as part of CMIP5 experiments were conducted using the same model configuration as above, except with time-evolving transient forcing of individual species (all forcing, GHGs, aerosols, and BC).Three ensemble members are available for each single-forcing run.In addition to the standard BC runs, we also conducted a new BC singleforcing simulation with adjusted BC emission factor as described in Sect.2.2.2.

Observations
The key model output in this high-altitude region to be compared with observations is surface temperature and snow cover.For temperature trend, we adopted both in situ data recorded at meteorological sites as reported in previous studies (Sect.2.4.1) and a high-resolution temperature reanalysis data set (Sect.2.4.2).For snow cover, we adopted a long-term data set (Sect.2.4.3) as well as the direct satellite measurement but only dated back to 2000s (Sect.2.4.4).The details of these observational data sets are below.

Ground-based temperature record
Monthly mean daily-minimum temperatures from 116 weather stations in the eastern Tibetan Plateau and its vicinity (with elevations ranging from 300 to 5000 m) during 1961-2006 are reported in Liu et al. (2009).Liu et al. (2009) only analyzed daily-minimum temperature because a wellrecognized feature associated with climatic warming is less warming observed in maximum temperatures and substantially more warming in minimum temperatures (Easterling et al., 1997).Previous studies also show such asymmetric changes in maximum and minimum temperatures are particularly true for Tibet (Liu et al., 2006) and the Alps (Weber et al., 1997).

Surface temperature reanalysis data set
Global Historical Climatology Network (GHCN) data set provides high-resolution (0.5 • by 0.5

NOAA climate data record of snow cover extent (Robinson et al., 2012)
Prior to 1999 the NH snow cover extent is based on satellitederived maps produced weekly by trained NOAA meteorologists.After 1999 NOAA NH snow cover extent maps were replaced by output from the Interactive Multi-sensor Snow and Ice Mapping System (IMS) processed at Rutgers University.Data are downloaded from http://climate.rutgers.edu/snowcover/docs.php?target=datareq.

Moderate Resolution Imaging Spectroradiometer (MODIS) snow cover observations
The MOD10CM product is a climate modeling grid product at a 0.05 • resolution with global coverage and monthly availability.Pixel values depict the percentage of snow cover (Hall et al., 2006).

Observed snow cover reduction linked with BC
The observations (Robinson et al., 2012) show that the snow cover extent over the Himalayan mountain range has declined at a rate of more than 10 % per decade since the 1960s (Fig. 1).The snow cover retreat along the Himalayan mountain range is greater than in Eurasia during the same period.In situ studies on regional glaciers and snowpack also reported strong declining trends.For example, Ma and Qin (2012) used 754 stations in China to document statistically significant declining trends of spring snow for the Qinghai-Tibet Plateau for the 1951-2009 period.Consistently, permafrost degradation has been reported on the Tibetan Plateau (Cheng and Wu, 2007;Li et al., 2008).Several satellite observations since the year of 2001 provided additional record in snow cover extent.The observed trend over this shorter period (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012) is less significant (Fig. S1a in the Supplement) and negative trends are only found along some portion of the Himalayan range.Consistently, the 5 km MODIS data set (Fig. S1b) also shows that the snow cover extent averaged over the entire Tibet region only has a slight decrease.But as other studies have pointed out (Pu et al., 2007;Pu and Xu, 2009), the highest altitudes of 5750-6000 m exhibit larger negative trends (−6 % decade −1 ).We note that at a shorter timescale (10 years), the snow cover trend is heavily influenced by natural variability and less significant.For example, during 1980-1991 (Fig. S1c) or 1990-2001 as shown in Fig. 5 of Menon et al. (2010), the declining trends are much larger.Nevertheless, the declining trend in the 40-50 year timescale (Fig. 1) is more robust and warrants further investigation of its causes, which is the main objective of this study.
To understand the causes of the observed trends of snow reduction over the multi-decadal timescale, we conducted global climate model simulations, in which BC emissions, CO 2 concentration or SO 4 emissions are increased instantaneously from preindustrial to present-day levels.Figure 2 (left column) shows the simulated change in snow fraction due to the increase in BC, CO 2 , and SO 4 aerosols.The pattern of snow cover decline in the BC model simulation captures the broad features of the observed decline (Fig. 1), with the largest snow reduction along the mountain range.The Tibetan Plateau on average showed a reduction in snow fraction of 1.9 % due to BC.The snow fraction shrinks by 2.9 % due to present-day CO 2 .Along the Himalayan mountain range, where the near-permanent snow cover exists, the reduction of snow fraction exceeds 10 % in both BC and CO 2 cases.Menon et al. (2010) attempted to simulate the snow reduction trends during the 1990s but the spatial distribution of the observed trend was not well captured mainly due to the coarse resolution of the model.Qian et al. (2011) also acknowledged their model's limitation in representing the snow cover climatology and therefore may have biases in estimating BC impact on snow.It is well known that global models tend to overestimate the snow cover of the Tibetan Plateau, and one potential reason is that the blocking effect for the moisture transport crossing the Himalayas is too small due to the coarse resolution of the global models and too much snowfall is simulated (Ménégoz et al., 2013).This limitation can partly be overcome with models using higher spatial resolutions.The modeling work presented here is a major step forward in terms of spatial resolution (about 1  Middle: same as left, but for surface albedo (%).Right: same as middle, but for snow depth (water equivalent, cm).The regionally averaged statistics are shown in Table 2. 2009, andQian et al., 2011;4 • by 5 • in Menon et al., 2010), which helps to better resolve the complex topography in this region.As a result of increased spatial resolution and also the improved land scheme, the biases in snow cover simulation is significantly reduced from its earlier model versions (Lawrence et al., 2011; also contrast Fig. S2c with Fig. 2 of Qian et al., 2011).However, we note that the precipitation over the Tibetan Plateau is still overestimated (Fig. S2b), and future studies, especially using regional climate models with even higher resolutions, are needed to improve the fidelity of model simulations of snowpack and glaciers over this topography-complicated region.
One consequence of the snow cover reduction is the decrease in surface albedo, which provides a positive feedback mechanism to localized warming.Such a surface albedo change in response to sea-ice loss has been observationally detected (Kay and L'Ecuyer, 2013;Pistone et al., 2014) and is important in explaining amplified Arctic warming.Flanner et al. ( 2011) also used observations during recent decades to calculate the surface albedo feedback in Northern Hemisphere large-scale snow-covered regions.Our simulations show that surface albedo over the Tibet region decreased by over 2 % (Fig. 2, right column) in response to BC.The maximum reduction occurs right along the Himalayan mountain range and part of the Tibet-Sichuan mountain regions.
The surface albedo decrease due to CO 2 shares a similar spatial pattern with BC (Fig. 2, right column) but with a smaller magnitude (Table 2b).Moreover, the snow depth reduction in response to CO 2 is only 30 % of that due to BC (Table 2c and Fig. 2), and this highlights the larger effect of BC in causing the regional cryospheric change over the Himalayas and Tibet.Not surprisingly, in the simulations the snow cover and surface albedo are increasing in response to cooling aerosols like SO 4 (Fig. 2), but in a smaller magnitude than that of the decreases due to BC and CO 2 .

Warming at high altitudes enhanced by BC
The Tibet region has witnessed increasing surface temperature by 0.3 • C per decade -more than twice the global average (Wang et al., 2008).One feature of the surface-warming trend over Tibet is that the warming magnitude increases significantly with altitude (Liu et al., 2009).To understand this anomalous feature, we show in Fig. 3 the tropospheric temperature responses (as a function of altitude and latitude) to BC as well as CO 2 and SO 4 .BC-induced heating rate (Fig. 3a) is more concentrated over the Northern Hemisphere (NH) due to larger emissions there from industrial activities, consistent with radiative forcing distribution (Table 1).The notable feature of BC response is the elevated warming at altitudes of 4000 to 8000 m and 30 to 60 • N (Fig. 3b), in particular over the Tibetan Plateau.The CO 2 warming pattern (Fig. 3b) features an amplified warming at the surface of the Arctic and in the upper tropical troposphere.CO 2 -induced warming in the upper atmosphere (500 hPa) over Tibet is 1 • C, larger than BC-induced warming of 0.5 • C, but the vertical gradient is much smaller (Fig. S3).SO 4 cooling fea-    tures an even stronger north-south asymmetry (north cooling; south slightly warming) but is more confined to the surface (Fig. S3).The temperature response in the troposphere is associated with strong meridional circulation change.The mechanisms behind the free-atmosphere circulation change, especially for the SO 4 case that does not have strong atmospheric forcing, are discussed in detail in Xu and Xie (2015).
Ground-based observations have shown that the last three decades were subject to a factor of 2 greater warming in the high-altitude interior of the Tibetan Plateau than at the edge of the plateau and at lower altitudes.The observations in Liu et al. (2009) were made between 1965 and 2006 from ground meteorological stations on the Tibetan Plateau region, and they revealed clear altitude dependence in the daily-minimum surface temperature (purple line in Fig. 4).The vertical profile of temperature change based on daily-average measurement from another reanalysis data set (Fig. 5) also reveals similar altitude dependence.What is driving the larger warming at high-altitude regions?
Figure 4 shows the model-simulated change of the dailyminimum surface temperature as a function of elevation due to three different forcing agents (CO 2 , SO 4 and BC).The surface temperature responses are calculated from all of the model grid cells over the Tibetan Plateau and the surrounding region (20-50 • N, 70-110 • E) to capture the altitude variation in this region.As shown in Fig. 4, the altitude dependence of the surface warming is mostly determined by the response to BC forcing (red dots).At altitudes below 1000 m the warming is minimal, but with increasing altitudes the magnitude of the warming increases up to 2 • C at 5000 m.The dependence of the surface warming on altitude is much smaller in the CO 2 case, which only increased from 1.2 • C warming at low altitudes to 1.6 • C warming at higher altitudes (yellow dots).The combined temperature response (black triangles in Fig. 4) by adding the individual trends due to BC, CO 2 , and SO 4 is largely consistent with the observed trend.To test the additivity of the temperature response, we conducted another simulation in which all of the three forcings were imposed simultaneously.The warming profile simulated by the combined anthropogenic forcing experiment largely agrees with the sum of the individual responses within 30 % (Fig. 5).Some non-linearity is expected as discussed in other modeling studies (Ming and Ramaswamy, 2009).The agreement of the simulated and the observed warming profiles provides a qualitative estimate of the relative contributions of BC, CO 2 , and SO 4 .Over the entire Tibetan Plateau, CO 2 -induced surface warming is 1.3 • C, compared to the BC-induced warming of 0.84 • C (Table 1b).Almost half of the surface warming at the highest altitudes (around 5000 m) is due to BC.
A potential complexity arises due to the internal variability of the climate systems, which has been shown to be im-portant in determining decadal trends over individual regions (Deser et al., 2012).To examine the role of natural variability, we calculated the temperature trend from all 350 consecutive 45-year periods out of 394 years of simulation (i.e.year 1-45, year 2-46,..., year 350-394).The 80 % (10th to 90th percentile) probability range of temperature change is shown in Fig. 5 (light-blue shading).The magnitude of the warming rarely exceeds 0.5 • C in any 45-year period in the long-term preindustrial control simulations without any external forcing.Therefore we infer that the vertical gradient of the temperature trend found in the simulations is very unlikely due to natural variability.One further concern is that the internal variability deduced from the long-term preindustrial control simulation can be model dependent.However, comparing a 30-member ensemble of CESM1 simulation (the same spatial resolution and configuration as in our control simulation) with the 38-member CMIP5, Kay et al. (2015)  in CESM1 ensemble is statistically the same with the spread in trends within CMIP5.Note that the climate simulations shown in Figs. 4 and 5 are driven by the instantaneous increase in present-day forcing.Since the real forcing trends were time dependent, we further analyzed a set of 20th century transient simulation output from the same model (Fig. 6), as part of CMIP5.The relative contributions of CO 2 and SO 4 to the simulated warming are consistent between the two sets of simulations (the instantaneous forcing and transient forcing).But the trends estimated from the 20th century transient forcing simulations are smaller than the quasi-equilibrium response to instantaneous forcing (parentheses in Table 1b).The reason is that only 70 % of SO 4 forcing and about 60 % of CO 2 forcing in the transient simulation were applied after 1960.The standard BC forcing (red solid line in Fig. 6) only lead to a weak warming, not exceeding the range of natural variability.As a result, the combined all-forcing responses (black triangles) did not capture the altitude dependence in observations very well.
Only when we adjusted historical BC forcing using the same scaling factors constrained by present-day observations, transient BC forcing induced a robust warming and amplification over high altitudes (red dots in Fig. 6), similar to what is shown in the instantaneous forcing experiment.However, note that the historical time dependence of the BC forcing is more uncertain, and we were also only able to produce one ensemble of adjusted BC simulation that was more subject to the influence of decadal variability.Therefore, the response to the instantaneous present-day BC forcing seems a more reliable indicator of the BC effects.While the absolute values of the warming profile need more model tests, our inference regarding the relative role of BC and CO 2 in the observed decrease in snow cover as well as the major role of BC on the altitude dependence of the warming trends is robust.

Physical mechanisms of elevated warming due to BC
Both CO 2 and BC contribute to the elevated warming at 5 km as shown in Fig. 4.However, BC is mostly responsible for the vertical gradient of the simulated warming trend.Most of the BC aerosols in the region are emitted over India and China and subsequently transported to the Tibetan Plateau and the Himalayan mountain range.The physical mechanisms for the amplified warming at higher altitude due to BC are at least threefold.

Direct heating in the atmosphere
BC absorbs a significant amount of solar radiation, as much as 25 % in typical pollution events as directly measured by multiple unmanned aircrafts over the northern Indian Ocean (Ramanathan et al., 2007).The BC layer placed at higher altitude is even more efficient in absorbing solar radiation than at sea level, due to stronger solar radiation and the brighter underlying cloud surface.In our model simulation, the BC atmospheric heating rate is concentrated in the Northern Hemisphere (maximum at 30 • N), coincident with the location of the maximum temperature change (Fig. 3).The elevated BC layer, due to the topography of the Tibetan Plateau and the Himalayan mountain range, contributes to the elevated heat- ing which is more than 0.1 • C day −1 and about 0.03 • C day −1 at 4 km (Fig. S3a).Such anomalous heating in the atmosphere over the elevated regions will contribute to the loss of ice and snow in two ways: (a) it will increase melting of the glaciers and snowpack, and (b) more of the precipitation will fall as rain instead of snow.CO 2 increase also induces longwave heating of the atmosphere (Fig. 3c), but it is well known that warming enhancement at the upper tropical troposphere is mostly due to moist convection processes (Manabe and Wetherald, 1975), and the warming enhancement at high altitudes is not showing a sharp gradient as in the BC case (Figs.S3b and 4).

Surface darkening by BC deposition
Snow and ice have a high surface albedo and reflect as much as 50 to 90 % of incoming solar radiation.Transported BC aerosols over the Himalayas and the Tibetan Plateau are removed from the atmosphere due to precipitation.When BC aerosols are deposited over the snow and ice, they increase the absorption of solar radiation and cause surface warming (Wiscombe and Warren, 1980;Chýlek et al., 1983).Recent studies have also suggested the influence of BC aerosols over regions like the Alps (Painter et al., 2013) and Eurasian land (Flanner et al., 2009).Menon et al. (2010) found that when the model includes snow albedo change due to BC the snow cover reduction is twice as large as the simulation with BC atmospheric heating effect only.Flanner et al. (2009) also suggested that BC surface albedo darkening effects are important in causing Eurasian springtime snow cover decline and are comparable to that of CO 2 .
The surface radiative forcing due to BC deposition over Tibet in this model is estimated to be 4.6 W m −2 based on the 5-year fixed-SST simulation (Fig. S4b).Because of this strong positive surface forcing associated with surface darkening, the shortwave forcing due to BC at the surface increased from −1.5 W m −2 (initially due to BC dimming effect) to a positive value of 3.1 W m −2 .This positive forcing imposed directly at the surface is even larger than the adjusted atmospheric heating due to BC over Tibet (1.6 W m −2 ).A recent modeling study by Ménégoz et al. (2014) also examined the role of BC deposition over snow in this region (with smaller forcing estimates of 1 to 3 W m −2 ), but their study did not estimate the atmospheric heating effect of BC.
However, we note that model estimates of radiative forcing due to BC deposition on snow have large uncertainty.Using the same atmospheric and land model (but driven by realistic meteorological field in the year of 2000), Zhang et al. (2015) showed that simulated BC concentration in snow is biased high with respect to in situ sampling.Although the large spatial variations in BC deposition can affect the representativeness of BC-in-snow measurements for the model evaluation purposes, this potential model bias should be kept in mind.In addition to the uncertainty in BC loading, the forcing magnitude is also sensitive to model parameterization (Yasunari et al., 2013), as well as the simulated background snow cover because wrongly simulated melting dates of the snowpack can lead to an incorrect radiative forcing (Jacobi et al., 2015).Therefore, both in situ (Ming et al., 2012;Wang et al., 2013;Zhao et al., 2014) and laboratory measurements (Hadley and Kirchstetter, 2012) are needed to constrain model representation of BC in snow.

Snow albedo feedback
The melting snow in response to the two initial heating mechanisms discussed above will further decrease surface albedo and increase solar absorption at the surface.The results based on the 60-year coupled model simulation suggest that the surface albedo will further decrease by 1.4 % and effectively impose an additional 3.2 W m −2 shortwave forcing at the surface.In summary, the elevated heating and surface darkening due to BC are simultaneously causing local warming and snow melting.The snow cover reduction further reduces surface albedo and then provides a positive feedback.A look at the seasonality of snow depth change suggests the early spring melting is important for this feedback.The net result of such a positive loop is an amplification factor of 4 for BCinduced Tibet warming from the global average values and significant snow and ice retreat.
Beyond the three main factors as we have discussed above, the changes in water vapor and clouds are also possible mechanisms contributing to the elevation-dependent warming in the mountain regions.As shown in the schematic of a recent review paper (Mountain Research Initiative EDW Working Group, 2015), in a warmer and moister atmosphere, the latent heat release at the cloud condensation level may induce larger warming at high altitudes (cloud feedback) and the downward longwave radiation increase particularly fast in a higher and drier atmosphere (water vapor feedback).It is difficult to identify or separate the contribution of these individual feedbacks from our current experiment setup.However, we note that these feedback mechanisms are operating regardless of forcing agents and therefore cannot explain the particularly large elevated warming in response to BC. Lastly, it is also worth commenting on the role of other snow impurities.In this study we used BC, a strong solar radiation absorber, to understand the climate response and the mechanisms due to absorbing aerosols that also include dust (Painter, et al., 2007;Di Mauro et al., 2015;Gabbi et al., 2015) and organic aerosols (Qian et al., 2015).Similarly, we used SO 4 to characterize all other reflecting aerosols.Any changes in dust and organics may induce changes to the snow cover, as their atmospheric heating and surface deposition are readily captured by this model, although the magnitude of response might be smaller since they are partially reflecting as well.

Conclusions
The observed surface warming over the Tibetan and Himalayan region of about 0.5 • C at sea level to about 2-2.5 • C at 5000 m (from 1961 to 2006) has been an outstanding feature of climate trends.The more than 2 • C warming is close to the peak warming trend observed anywhere on the planet.For comparison, the Arctic warming associated with large sea-ice retreat during this period is 1.2 • C.
The high-resolution coupled ocean-atmosphere model in this study was able to attribute the observed warming trends and their high-altitude enhancement to imposed increases in CO 2 , BC, and SO 4 aerosols.The simulated changes with all forcing imposed were consistent with the observations.The key to the success is that we obtained the BC forcing from the reconstruction of ground-based and satellite-based observations.The imposed BC forcing was about 2 to 4 times (depending on the regions) larger than that simulated by the models using bottom-up emission inventories.The analysis of model simulations highlights that the high-altitude warming due to BC is as large as CO 2 warming over the Tibetan Plateau and the elevated warming profile is unique in BC responses.
The observed record warming is accompanied by retreat of glaciers and snow cover as well as thinning of the snowpacks.In response to the preindustrial to the present-day increase in BC emissions, the annual averaged snow fraction over the Tibetan Plateau is reduced by more than 6 % (relatively), and the snow depth by approximately 19 %.The surface albedo decreases by more than 5 % along the Himalayan mountain range and 1.4 % over the entire Tibet region, providing a positive local feedback to the enhanced local warming.In stark contrast, despite having 5 times larger effect in global mean temperature than BC, over Tibet CO 2 impact is only 1.5 times stronger in snow cover decrease, and only one-third in snow depth decrease.We conclude that BC is instrumental in causing snow retreat and its effects are manifested simultaneously through a threefold process: (i) direct atmospheric heating, (ii) darkening of the snow surface, and (iii) the snow albedo feedback.It is important to note that, without the scaling factor we applied to bring the model BC forcing to the observationally constrained values, the impact of BC on the observed temperature trends would have been marginal.This perhaps explains why the models used in IPCC assessments have not simulated the role of BC in the large warming trend over the Himalayas.
The Supplement related to this article is available online at doi:10.5194/acp-16-1303-2016-supplement.
For the period March 2000 to December 2006, the algorithm version 4 is used and after that version 5 of the algorithm is used.Snow cover products derived from MODIS are based on a spectral ratioing of MODIS band 4 (green) (0.545-0.565 µm) and band 6 (near infrared) (1.628-1.652µm).Data are downloaded from http://nsidc.org/data/MOD10CM.

Figure 1 .
Figure 1.Observed snow cover extent change (% per decade) from 1967 to 2012.The trend is calculated based on snow cover extent data in the entire period.Insignificant changes (confidence interval < 75 % calculated from a Student t test) are not shown.

Figure 2 .
Figure 2. Left: (a) simulated change in snow fraction (%) due to present-day BC (versus its preindustrial level), (b) CO 2 , and (c) SO 4 .Middle: same as left, but for surface albedo (%).Right: same as middle, but for snow depth (water equivalent, cm).The regionally averaged statistics are shown in Table2.
• N and 80 to 100 • E. (b) Surface temperature change ( • C) in response to different forcings in (a).Surface temperature change is calculated by averaging the last 60 years of a 75-year coupled model simulation.The values in parentheses are temperature change in the 20th century time-dependent forcing simulations.The linear trend ( • C decade −1 ) is first calculated and then multiplied by 4.5 to obtain the change with 45-year time frame.BC responses include the range of using "standard" and adjusted emissions.

Figure 3 .
Figure 3. Left: globally zonal averaged radiative heating rate ( • C day −1 ) as a function of altitude and latitude due to (a) BC, (b) CO 2 , and (c) SO 4 , calculated from the 5-year fixed-SST simulations using the instantaneous radiative diagnostic procedure.Shortwave fluxes are shown for BC and SO 4 , and longwave flux for CO 2 .Right: the temperature response ( • C) due to (a) BC, (b) CO 2 , and (c) SO 4 , calculated as the difference of the last 60 years of 75-year perturbed simulation and the 319-year long-term control.Figure S3 shows the normalized heating rate and temperature profile averaged just over the Tibetan Plateau.

Figure 4 .
Figure 4.The change of daily-minimum temperature ( • C) as a function of elevation (km).The observation from 1961 to 2006 are from of Liu et al. (2006).The simulated temperature responses due to instantaneous increase in forcings (CO 2 , SO 4 and BC) are calculated from model grid cells over the Tibetan Plateau and its vanity region (20-50 • N, 70-110 • E) including low-lying regions and high-altitude regions.The standard deviation due to spatial variation of temperature response is shown as error bars.The sum of CO 2 , SO 4 , and BC responses are shown in black triangles.

Figure 5 .
Figure 5. Similar to Fig. 4, but with the following differences: (1) daily-mean surface temperature, not daily-minimum temperature, are shown; (2) the spread of five-ensemble member of BC simulations are shown in red lines; (3) the observations are from GHCN data set (1961 to 2006); (4) the range of temperature change found in unforced preindustrial control simulations is shown in blue shading; and (5) the all-forcing simulation (black line) is shown in comparison with the sum of individual responses (black triangles).

Figure 6 .
Figure 6.Similar to Fig. 5 but showing results from the 20th century transient simulations (three ensemble members for each singleforcing run).Note that the "standard" BC single-forcing simulation used smaller BC emissions as in other CMIP5 models (red solid line).An additional simulation with adjusted larger BC emissions is shown (red dotted line, one ensemble member only).

2016 1306 Y. Xu et al.: Observed high-altitude warming and snow cover retreat due to black carbon itoring
• ) global land surface temperature records from 1948 to near present (Fan and van den Dool, 2008).The data set uses a combination of two large individual data sets of station observations collected from the GHCN version 2 and the Climate Anomaly Monwww.atmos-chem-phys.net/16/1303/2016/Atmos.Chem.Phys., 16, 1303-1315, System.Data are downloaded from http://www.esrl.noaa.gov/psd/data/gridded/data.ghcncams.html.

Table 2
. (a) Snow fraction (%), (b) surface albedo (%), and (c) snow depth over land (water equivalent, cm) change in response to different forcings.The relative change as a percentage is shown in parentheses next to the absolute change.