1. Introduction
The High Mountain Asia (HMA) region covers an area of 3.8 million km2 (Baumann and others, Reference Baumann, He, Schmidt, Kühn and Scholten2009) and encompasses parts of Afghanistan, Pakistan, India, Nepal, China, Tibet, Bhutan, Kazakhstan, Uzbekistan, Kyrgyzstan and Tajikistan. This region is also known as the ‘third pole’ due to hosting the largest concentration of ice outside the polar regions (IPCC, 2013). Here lie the headwaters of some of the major rivers of the Asian continent, such as the Indus, Ganges, Brahmaputra, Yangtze and Yellow (Pritchard and others, Reference Pritchard, Forsythe, Fowler, O’Donnell and Li2019). Collectively, these hydrographic basins provide water to more than a billion people, constituting over 20% of the global population (GPWv3, 2005; Immerzeel and others, Reference Immerzeel, Van Beek and Bierkens2010, Reference Immerzeel2020; IPCC, 2013).
The meltwater generated by the glaciers and snow cover in HMA plays a critical role in safeguarding this vast population from severe water scarcity during summer or droughts. The additional flow integrates into the river systems, offsetting the decrease in discharge resulting from the reduced rainfall amount and frequency. Each summer, the volume of meltwater produced in HMA is ∼36 km3, increasing to 37 km3 in drought years (Pritchard and others, Reference Pritchard, Forsythe, Fowler, O’Donnell and Li2019). This buffering effect is more pronounced in the most glacierized areas of HMA (i.e., the Hindu-Kush, Karakoram and Himalayan ranges) which house extensive valley glaciers. Therefore, the contributions of glacial meltwater to the rivers are indispensable for sustaining services such as water supply, agricultural irrigation and hydropower generation.
Unfortunately, HMA cryosphere has shown degradation signals since the mid-19th century (Bolch and others, Reference Bolch2012), a situation that has become aggravated in the 21st century, as reflected in some studies conducted in the last two decades (e.g., IPCC, 2019; Zemp and others, Reference Zemp2019). This loss of critical glacier mass will inevitably lead to a long-term decrease in glacier-derived meltwater runoff. In the past century, droughts were the most devastating natural disasters in the region, resulting in ∼6 million deaths and affecting ∼1.1 billion people (National Research Council, 2012). Moreover, reduction in food production due to water deficit can lead to social instability, sudden and uncontrolled migrations and intra- and interstate conflicts (World Economic Forum, 2016).
One of the degradation factors is the emission of black carbon (BC) into the atmosphere, an aerosol with a high carbon content that can be generated by human activities through the incomplete combustion of fossil fuels or biomass burning (Bond and others, Reference Bond2013). The residence time of these BC particles is typically ∼1 week (Ramanathan and Carmichael, Reference Ramanathan and Carmichael2008; Lee and others, Reference Lee2013), facilitating their dispersion over distances that can exceed 1000 km (Rodhe and others, Reference Rodhe, Persson and Åkesson1972; Bond and others, Reference Bond2013). The HMA region is surrounded by important global BC emitting centres like East, Southeast and South Asia, with countries like China and India standing out (Lu and others, Reference Lu, Streets, Zhang and Wang2012; Bond and others, Reference Bond2013; IPCC, 2013). When BC particles are deposited on the snow or ice surface, they tend to absorb a higher proportion of shortwave solar irradiance, reducing the surface albedo, warming the particles and eventually increasing the melting rate of the snow (e.g., Flanner and others, Reference Flanner, Zender, Randerson and Rasch2007; Ménégoz and others, Reference Ménégoz2014) and ice (e.g., Xu and others, Reference Xu2009a, Reference Xu2009b; Kaspari and others, Reference Kaspari2011). Therefore, due to its geographical context, HMA is one of the most vulnerable regions in the world to the effects of BC.
This study examines the impact of BC on snowmelt during the 2019 ablation season. Samples were collected from two snow pits excavated in the seasonal snowpack on the Godwin-Austen Glacier, as well as from surface snow at K2 Camps 1 and 2. To provide additional context, measurements of water stable isotopes (δ1⁸O and δ2H) were used to evaluate the representativeness of the snowpack as a record for the 2018–19 accumulation season, enhancing interpretations of likely BC and moisture sources and the timing of BC deposition.
2. Study framework
The Godwin-Austen Glacier in the Central Karakoram is part of the Upper Indus Basin (UIB; Fig. 1a). This basin has an extensive glaciated area, covering ∼22 000 km2 (Immerzeel and others, Reference Immerzeel, Van Beek and Bierkens2010) and faces dry summers and moderately wet winters. Approximately 49% of the annual precipitation falls as snow (Pritchard and others, Reference Pritchard, Forsythe, Fowler, O’Donnell and Li2019) that covers ∼85% of the basin area in winter (Shrestha and others, Reference Shrestha2015). Various studies have reported that more than half of the annual Indus River flow comes from the snowpack (seasonal or perennial) and glacier melting (e.g., Hewitt and others, Reference Hewitt, Wake, Young and David1989; Wake, Reference Wake1989; Young and Hewitt, Reference Young and Hewitt1990; Archer and Fowler, Reference Archer and Fowler2004; Bookhagen and Burbank, Reference Bookhagen and Burbank2010). Particularly noteworthy is the contribution of seasonal snowpacks, as highlighted by the study of Yu and others (Reference Yu2013), which indicated that the meltwater of the glaciers located in the Karakoram, western Himalayas and Hindu-Kush contributes ∼18% of the total flow in the Indus basin, while the remaining 82% comes from the melting of seasonal snowpacks. In this context, the net annual volume of meltwater in the Indus River over the 4 month ablation season is 14.5 ± 3.1 km3 (Pritchard and others, Reference Pritchard, Forsythe, Fowler, O’Donnell and Li2019). The meltwater produced counteracts the reduction of river discharge caused by decreased precipitation during the thawing season or drought, making it a crucial resource. This is especially important considering that the economy of the Indus River lowlands, highly dependent on agriculture (Bhatti and others, Reference Bhatti, Suttion and Seigo2009), relies on the largest network of irrigation channels in the world, the IBIS (Indus Basin Irrigation System). The IBIS can extract up to 75% of the Indus River flow, making this basin one of the most depleted worldwide (Sharma and others, Reference Sharma2010).
Figure 1. (a) Geographical context of the Karakoram range (yellow square) and Indus river (light blue line). (b) The yellow rectangle indicates the study area within the glacier system composed of the Baltoro Glacier, Godwin-Austen Glacier and Baltoro South Glacier. (c) Position of the snow pits (P1 and P2) relative to the K2 Base Camp.
3. Methods
Between January and February 2019, we stayed at K2 Base Camp (K2BC; 5150 m a.s.l.), located on the western moraine of the Godwin-Austen Glacier (Fig. 1b and c). We excavated two snow pits on the glacier, selecting sites far from K2BC to minimize contamination from emissions. Based on the automatic weather station (AWS) data indicating prevailing west-to-east winds, we chose locations south of K2BC. Snow pit 1 (P1; 35°49.486’ N, 76°30.779’ E; 4937 m a.s.l.) was located 950 m south of K2BC and reached a depth of 120 cm, while Snow pit 2 (P2; 35°49.184’ N, 76°30.920’ E; 4927 m a.s.l.) was 500 m further south, reaching a depth of 200 cm (Fig. 1c). Both snow pits covered the full snowpack thickness (Fig. S1). Additionally, mountaineers collected surface snow samples at K2 Camp 1 (6076 m a.s.l.; 35°51.95’ N, 76°32.17’ E) and Camp 2 (6650 m a.s.l.; 35°52.15’ N; 76°31.99’ E).
3.1. In situ measurements and sampling
The snow pits were excavated in late February 2019, with P1 on 24 February 2019 and P2 on 27 February 2019. Works in P1 started at 10:56 UT (UTC + 5) and took ∼3 h, while P2 began at 11:50 UT (UTC + 5) and took around 5 h due to its greater depth. Snow temperature and density were measured at intervals of 5 and 8 cm, respectively, using a TFA® digital thermometer with a precision of ±0.5 °C, and a 125.5 cm3 box cutter based on the original design of the Institute of Low Temperature Science of the Hokkaido University. Snow samples for BC measurements and stable isotopes were collected at 6–8 cm intervals, with volumes of 1300–1800 cm3, stored in Whirlpack™ bags and melted at room temperature (∼10ºC). The melted samples were homogenized, transferred into 30 mL glass vials and kept at 1–5°C during shipping, as plastic bottles can alter isotopic signals during long-term storage (Spangenberg and Vennemann, Reference Spangenberg and Vennemann2008; Spangenberg, Reference Spangenberg2012).
3.2. Chemical analysis
The National Institute of Polar Research in Tokyo, Japan, performed the chemical analysis of the water samples. For BC measurements, a modified Single-Particle Soot Photometer (SP2; Droplet Measurement Technologies, Longmont, CO, USA) was used, employing the Laser-Induced Incandescence method to measure BC based on its refractory properties (e.g., Kaspari, Reference Kaspari2011; Sterle and others, Reference Sterle, McConnell, Dozier, Edwards and Flanner2013). The advantages of the SP2 include its lack of requirement for prior water sample filtering, low water volume requirements, the possibility of using a continuous flow system and it is not sensitive to other materials, avoiding potential sources of uncertainty associated with filtration methods. In this particular case, we used a Wide-Range Single-Particle Soot Photometer (Mori and others, Reference Mori2016), a modified version of the standard SP2 with a high-efficiency nebulizer (Marin-5, CETAC Technologies Inc., Omaha, NE, USA). This combination allows to measure BC particles with diameters between 70 and 4000 nm, whereas a standard SP2 would only allow us to measure BC particles with diameters between ∼70 and 500 nm.
Complementing the BC analysis, water stable isotopes (1⁸O and 2H, or deuterium) were measured using an Isotopic Ratio Mass Spectrometer, with a precision of ±0.5‰ to ±0.6‰ for δD and ±0.01‰ to ±0.11‰ for δ1⁸O (Uemura and others, Reference Uemura, Matsui, Motoyama and Yoshida2007). Isotopic data are particularly useful, as the δ1⁸O content in precipitation reflects atmospheric condensation temperatures, with colder conditions leading to lower δ1⁸O values, which can be preserved in the snowpack (Dansgaard, Reference Dansgaard1964; Friedman and others, Reference Friedman, Smith, Gleason, Warden and Harris1992; Mayer and others, Reference Mayer2014). Additionally, δ2H and δ1⁸O follow a global linear correlation known as the Global Meteoric Water Line (GMWL), with the y-axis intercept (10‰) referred to as deuterium excess or d-excess (Epstein and Mayeda, Reference Epstein and Mayeda1953; Friedman, Reference Friedman1953; Craig, Reference Craig1961; Dansgaard, Reference Dansgaard1964). Variations in atmospheric humidity during evaporation modify the GMWL and, consequently, the d-excess. Higher humidity at the source results in lower d-excess values, and vice versa, helping to identify moisture sources and trace the origin of water vapour (Jouzel, Reference Jouzel, Merlivat and Lorius1982; Rindsberger and others, Reference Rindsberger, Magaritz, Carmi and Gilad1983; Pfahl and Sodemann, Reference Pfahl and Sodemann2014; Galewsky and others, Reference Galewsky, Steen-Larsen, Field, Worden, Risi and Schneider2016).
3.3. Atmospheric reanalyses
For this study, it was essential to have meteorological records (e.g., temperature, precipitation, solar irradiance, cloud cover) from the Godwin-Austen Glacier during the winter of 2018–19 and the 2019 ablation season. To obtain that information, an AWS was set up at K2BC during the expedition (January 2018–March 2019), providing highly accurate data on temperature and wind speed and direction. Unfortunately, due to technical issues, the data were only available from 3 February onwards. Moreover, meteorological variables like temperature and precipitation vary significantly with elevation, and data tend to be biased because stations are typically located in valley bottoms (Wake, Reference Wake1989). For instance, on the K2 southern slope, precipitation ranges from 150 to 200 mm a−1 in the elevation range 2000–3200 m a.s.l., with values of 1600 at 6100 m a.s.l., and reaching 2500 mm a−1 at 8000 m a.s.l. (Decheng, Reference Decheng1978; Mayer and others, Reference Mayer, Lambrecht, Belò, Smiraglia and Diolaiuti2006). Consequently, we resorted to atmospheric reanalysis to obtain continuous meteorological data, specifically ERA5 (Hersbach, Reference Hersbach2016), as validated for the UIB by studies from Dahri and others (Reference Dahri2021) and Baudouin and others (Reference Baudouin, Herzog and Petrie2020). For regional atmospheric BC concentrations, we used NASA’s MERRA-2 data from the Global Modeling and Assimilation Office.
3.4. Back-trajectory analysis
To analyse the snowpack isotopic data and BC content from a global and regional atmospheric circulation perspective, we employed the Hybrid Single-Particle Lagrangian Integrated Trajectory model (HYSPLIT; Stein and others, Reference Stein, Draxler, Rolph, Stunder, Cohen and Ngan2015; Rolph and others, Reference Rolph, Stein and Stunder2017). This model calculates daily back-trajectories of air parcels before reaching the Godwin-Austen Glacier. For meteorological data, we used the Global Data Assimilation System, developed by NOAA in conjunction with the National Centers for Environmental Prediction (NCEP), with a 0.5° × 0.5° latitude/longitude resolution. From the start of the accumulation season in October 2018 (Archer and Fowler, Reference Archer and Fowler2004; Hasson and others, Reference Hasson2016) to the snow pit excavation on 24–27 February, we computed 180 back-trajectories at two altitudes (500 and 1500 m above ground level, a.g.l.) to capture variations in air mass origins across different heights, using the snow pits as endpoints. To simplify the analysis, we grouped the trajectories by month using Total Spatial Variance and the K-means clustering algorithm (Fig. S2; Hartigan and Wong, Reference Hartigan and Wong1979; NOAA, 2023a; NOAA, 2023b). Notably, different algorithms may yield variations in results, even when applied to the same dataset (Cui and others, Reference Cui, Song and Zhong2021). These clusters help identify probable regional BC emission sources, the origins of predominant air masses and potential moisture sources contributing to snowpack formation (Sinclair and Marshall, Reference Sinclair and Marshall2009).
Back-trajectories were classified considering the last major water body that air parcels passed over before reaching the study site. Although moisture gains and losses could occur along their path, the isotopic signature of the precipitation is likely predominantly influenced by the moisture from these large water bodies. Using this criterion, we established six domains: one circular (Domain 1) and five wedge-shaped (or radial; Domain 2–6), using snow pits coordinates as reference points (see Figure 4). Each radial domain consists of two regions based on the distance to the study site. This division is similar to the one proposed by Juhlke and others (Reference Juhlke2019) for the Pamir Mountains, based on the relative humidity data provided by HYSPLIT. Domain 1 is centred on the snow pits and covers a radius of ∼500 km, whereas the wedge-shaped domains are North Central Asia and Siberia (Domain 2), the Caspian Sea and Black Sea (Domain 3), the Middle East and Mediterranean Sea (Domain 4), the Persian Gulf and Red Sea (Domain 5) and Northern India and the Arabian Sea (Domain 6).
3.5. Albedo reduction, radiative forcing and snowmelt
To assess the impact of BC, we needed to determine the span of the 2019 ablation season in the Godwin-Austen glacier region. The removal of the K2BC weather station in early March 2019 left us without direct temperature data, so we relied on previous studies. On the one hand, our field observations in late February showed no signs of liquid water or refreezing layers, indicating the snowpack was still ‘dry’ at that time. On the other hand, studies by Lund and others (Reference Lund2020), Ahmad and others (Reference Ahmad2020), Mukhopadhyay and others (Reference Mukhopadhyay, Khan and Gautam2015) and Yi and others (Reference Yi, Liu, Zhu, Wu, Xie and Saifullah2021) indicate that snowmelt in the Karakoram generally spans from May to September. Consequently, this is an extrapolation from previous research rather than direct in situ data.
For estimating the effects of BC on snow albedo reduction, radiative forcing and melt rates, we used parameterizations by Dang and others (Reference Dang, Brandt and Warren2015, Reference Dang, Warren, Fu, Doherty, Sturm and Su2017), which account for BC particle concentration, shortwave solar irradiance, average snowpack ice crystal size (Flanner and others, Reference Flanner, Zender, Randerson and Rasch2007) and cloud fraction (CF; Dang and others, Reference Dang, Warren, Fu, Doherty, Sturm and Su2017). This method can only be applied to ‘semi-infinite’ snowpacks, thick enough to hide the underlying surface for all wavelengths (Dang and others, Reference Dang, Brandt and Warren2015), and has been already used in the Andes and Antarctic Peninsula (Cordero and others, Reference Cordero2022a, Reference Cordero2022b).
We calculated the albedo reduction (ΔA) under two atmospheric conditions: clear sky (ΔA clear) and cloudy sky (ΔA cloudy). The parameterizations consider the combined effect of BC particles and snow crystal size variation on ΔA. A given concentration (C) of BC in a snowpack with an average snow crystal radius (r) can produce an ΔA similar to those in a snowpack with lower BC concentrations and larger average radius. For instance, snow with a BC concentration of 20 ng g−1 and an average snow crystal radius of 100 μm reduces albedo by 0.5%, whereas if the radius is 1000 μm (and equal BC concentration), the albedo is reduced by 1.6% (Dang and others, Reference Dang, Brandt and Warren2015). In the infrared spectrum, snow albedo is primarily controlled by the size and shape of ice crystals (Wiscombe and Warren, Reference Wiscombe and Warren1980; Warren, Reference Warren1982). On the other hand, in the visible spectrum, albedo is more influenced by particles with high carbon content (e.g., BC), which increase photon absorption, thereby decreasing the reflective capacity of the snow (i.e., lower albedo; Chýlek and others, Reference Chýlek, Videen, Ngo, Pinnick and Klett1995).
As proposed by Dang and others (Reference Dang, Brandt and Warren2015), the r and C variables can be combined to calculate ΔA through the definition of the following parameter:
\begin{equation*}H = \left( {\frac{C}{{{C_0}}}} \right){\left( {\frac{r}{{{r_0}}}} \right)^s}\,{\text{,}}\end{equation*}where C 0 is a BC reference concentration (10−6 ng g−1) and r 0 is a reference crystal radius (100 μm). The exponent s takes the value 0.73 (Table 3 in Dang and others, Reference Dang, Brandt and Warren2015). The albedo reduction (ΔA) is related to H by a quadratic function:
\begin{equation*}{\text{lo}}{{\text{g}}_{10}}\left( {\Delta A} \right) = {m_1} \cdot {p^2} + {m_2} \cdot p + {m_3},\end{equation*}where
\begin{equation*}p = {\text{lo}}{{\text{g}}_{10}}H.\end{equation*}The coefficients m 1, m 2 and m 3 are defined in Table 3 of Dang and others (Reference Dang, Brandt and Warren2015) for both clear sky and cloudy conditions; therefore, we can compute the parameters ΔA clear and ΔA cloudy. It is important to note that these calculations are only valid for crystal radii ranging from 50 to 2500 μm. According to Dang and others (Reference Dang, Warren, Fu, Doherty, Sturm and Su2017), the total albedo reduction (ΔA) can be estimated by
\begin{equation*}\Delta A = CF \cdot \Delta {A_{{\text{cloudy}}}} + \left( {1 - CF} \right) \cdot \Delta {A_{{\text{clear}}}},\end{equation*}where CF represents the average cloud fraction over the melting season, ranging between 0 and 1, with 0 indicating the complete absence of clouds and 1 denoting full coverage. The radiative forcing of BC (RF BC) is computed by multiplying the total albedo reduction (ΔA) by the average incoming solar irradiance or shortwave irradiance (I):
\begin{equation*}R{F_{{\text{BC}}}} = \Delta A \cdot I.\end{equation*}The values of CF and I are derived from the ERA5 atmospheric reanalysis (Hersbach, Reference Hersbach2016). The RF BC is a product of an increased absorption of solar irradiance by BC particles, resulting in snowpack surface warming. This extra energy absorbed per unit surface of the snowpack (Ex) can be calculated by multiplying RF BC by the duration of the melting season (d):
\begin{equation*}{E_x} = R{F_{{\text{BC}}}}d.\end{equation*}Finally, we calculate W, which is the snow that melts sooner due to RF BC:
\begin{equation*}W = \frac{{\Delta A \cdot I \cdot d}}{E},\end{equation*}where E represents the water latent heat of fusion (Cohen, Reference Cohen1994).
4. Results
4.1. Temperature and mass density
Both snow pits (P1 and P2) showed similar temperature variations despite differences in depth (121 cm in P1 and 200 cm in P2). Near the surface (top 10 cm), temperature differences were minimal with −13.7ºC in P1 and −13.5ºC in P2, whereas the bottom showed a temperature of −9.5ºC in P1 and −9.1ºC in P2. A pronounced temperature decrease occurred in the upper layers, with a gradient of −0.32ºC cm−1 in P1 and −0.26ºC cm−1 in P2, reaching lows of −16.7ºC at 10 cm in P1 and −18.8ºC at 20 cm in P2. Below these points, temperatures decreased gradually at a rate of 0.065ºC cm−1 in P1 and 0.053ºC cm−1 in P2 down to the snowpack base (Fig. S3).
Average mass density was 292 kg m−3 in P1 and 333 kg m−3 in P2, with lower densities in the surface layers (211 kg m−3), some layers in the central section (249 kg m−3) and at the snow–ice interface (289 kg m−3). The low-density surface layer had crystal sizes ranging from <1.0 to 1.5 mm and exhibited irregular and rounded shapes. In contrast, the low-density layers in the central and basal sections had larger crystals (2.0–3.5 mm) characterized by well-defined faces and edges. Although the P2 basal low-density layer was thicker compared to P1, the crystals in P1 were slightly larger and faceted. Sampling the intermediate and basal low-density layers was quite challenging, as crystals tended to disintegrate due to low cohesion. In the remaining snowpack, the ice crystals ranged between 1.5 and 2.0 mm (Figs. 2 and 3).
Figure 2. Measurements along P1. (Left) Isotopic variations with δ18O (top) and d-excess (bottom). (Centre) variation in black carbon (BC) concentrations (top) with differentiation between zones Z1 (orange) and Z2 (blue), along with mass density fluctuations (bottom). (Right) Evolution of grain size range with depth and grain shape classification according to Colbeck (Reference Colbeck1990) and Fierz and others (Reference Fierz2009).
Figure 3. Measurements along P2. (Left) Isotopic variations with δ18O (top) and d-excess (bottom). (Centre) variation in black carbon (BC) concentrations (top) with differentiation between zones Z1 (orange) and Z2 (blue), along with mass density fluctuations (bottom). (Right) Evolution of grain size range with depth and grain shape classification according to Colbeck (Reference Colbeck1990) and Fierz and others (Reference Fierz2009).
4.2. Isotopic variations
Samples were collected at 6 cm intervals along P1 and 8 cm intervals along P2, with the average isotopic values assigned to the midpoints of the sampling intervals. Despite some data gaps, the resolution was adequate to establish correlations between the snow pits.
The δ18O had a similar trend in both snow pits. At the surface, the values were −22.2‰ for P1 and −21.7‰ for P2. As depth increased, we observed a gradual depletion of isotope values, reaching the minimum at 33 cm in P1 (−30.4‰) and at 75 cm in P2 (−30.5‰). From this point, δ18O in P1 increased until reaching the maximum enrichment (−16.3‰) at 108 cm and then descended to −16.5‰ at the snowpack base. The second snow pit (P2) reached its less depleted value (−17.8‰) at 192 cm before reaching the ice contact with a value of −18.5‰. Regarding d-excess, the values at the surface were also comparable: 10.2‰ in P1 and 10.5‰ in P2. As depth increased, d-excess decreases until reaching its minimum values of 6.1‰ in P1 at 27 cm depth and 4.4‰ in P2 at 68 cm depth. With further depth increase, values progressively rose, peaking at 18.6‰ at 63 cm in P1 and 18‰ at 140 cm in P2. Towards the snowpack bottom, the d-excess profiles of the two snow pits differed somewhat. The snow pit 2 (P2) simply undergone a progressive decrease to 13.8‰. In contrast, the P1 values first decreased to a minimum of 12.5‰ at 105 cm, and then experienced an increase to 18.3‰ at the base (Figs. 2 and 3).
4.3. Back-trajectory analysis
Following the proposed regional classification (Fig. 4) and considering the different altitudes at which the back-trajectories were evaluated, significant variability was observed in September 2018 at 500 m a.g.l. Local back-trajectories predominated (47%), followed by those from North Central Asia and the Caspian Sea (13% each), while the Mediterranean Sea and Northern India contributed 10% each, and the Arabian Sea accounted for the remaining 7%. At 1500 m a.g.l., local trajectories even more dominant (70%), with the Mediterranean Sea (13%) and Northern India (17%) being the only other notable contributors, and no trajectories from North Central Asia, the Caspian Sea or the Arabian Sea. Although September 2018 is not part of the accumulation season, we included this month to observe the preceding back-trajectory trends. In October 2018, the pattern at 500 m a.g.l. showed an increase in local trajectories (52%) and those from the Caspian Sea (35%), with North Central Asia steady at 13%. At 1500 m a.g.l., compared to 500 m a.g.l., local trajectories increased to 65%, while contributions from the Caspian Sea and North Central Asia decreased to 32% and 3%, respectively. By November 2018, 500 m a.g.l. trajectories were more diverse, with the Middle East (28%) becoming the dominant source, and the Caspian Sea, Mediterranean Sea and Persian Gulf each contributing 20%, alongside a minor influence from Siberia (12%). In contrast, at 1500 m a.g.l., the Mediterranean Sea dominated (67%), followed by the Middle East (33%), with no contributions from the Caspian Sea, Persian Gulf or Siberia. December 2018 at 500 m a.g.l. saw dominance from the Black/Caspian Sea regions, but at 1500 m a.g.l., the Mediterranean Sea was the sole contributor, a trend that continued into January 2019. In January, the influence of the Black/Caspian Sea decreased to 55% at 500 m a.g.l., with the Mediterranean Sea (13%) and Persian Gulf (32%) reappearing, though at 1500 m a.g.l., neither the Persian Gulf nor Caspian Sea were present. Finally, in February 2019, at 500 m a.g.l., the Caspian Sea was absent, replaced by the Middle East (46%). Mediterranean Sea comprised by 43%, with smaller contributions from the Persian Gulf (7%) and Siberia (4%). At 1500 m a.g.l., the Mediterranean Sea (40%) and Middle East (25%) origins persisted, with a reappearance of the Caspian Sea (29%) and no contribution from the Persian Gulf, while Siberia remained steady at 4% (see bars in Fig. 6).
Figure 4. (Top) Geographic divisions used for back-trajectory classification (domains; D1 to D6). (Bottom) Clusters generated from 180 daily back-trajectories (September 2018 to February 2019) at 500 m a.g.l. (black lines) and 1500 m a.g.l. (red lines).
4.4. BC concentrations
For BC quantification, the snow pit profiles (P1 and P2) showed data gaps similar to those in the isotopic record; additionally, we collected surface snow samples from K2 high-altitude Camps 1 and 2. In P1, the BC concentration ranged from 5 to 44 ng g−1, with an average value of 14 ng g−1. Meanwhile, in P2, the concentration ranged from 2 to 28 ng g−1, with an average of 10 ng g−1. If we combine the BC concentrations of both snow pits, the average value is 12 ng g−1. Based on concentration levels and fluctuation patterns, we divided the BC profiles into two main zones: Zone 1 (Z1) and Zone 2 (Z2). In P1, the average BC concentration in Z1 was 8 ng g−1, with this zone extending to a depth of 77 cm. Conversely, in P2, the Z1 average concentration was comparable to P1 (8 ng g−1), but it extended deeper, reaching 160 cm. Zone 2 (Z2) expanded from the limit with Z1 to the snowpack–glacier ice interface and exhibited a higher variability and values. In P1, Z2 average concentration was 27 ng g−1, with maximum values reaching 44 ng g−1, while in P2, the average concentration was 18 ng g−1, with maximum values of 29 ng g−1. Finally, regarding the high-altitude camps, the average concentration for Camp 1 was 7 ng g−1, and 26 ng g−1 for Camp 2 (Figs. 2 and 3).
Observations in the snow pits revealed darker layers, especially in P1, where a basal-dark layer aligns with higher BC concentrations (Fig. S1). Another layer at ∼80 cm depth lacked chemical data, leaving its cause – BC or mineral dust – undetermined. Although beyond the scope of this article, mineral dust, common in HMA due to the nearby Taklamakan Desert, also contributes to snow darkening, increasing albedo and accelerating melt (Ge and others, Reference Ge, Huang, Xu, Qi and Liu2014; Jia and others, Reference Jia, Liu, Chen, Zhang and Huang2015; Chen and others, Reference Chen2017; Yuan and others, Reference Yuan2019; Xing and others, Reference Xing2024).
4.5. BC effects on the snowpack
According to the daily data provided by ERA5 (Hersbach, Reference Hersbach2016), during the 2019 melting season, the shortwave solar irradiance (I) varied between 242.24 and 0.0035 W m−2, with an average value of 106.68 W m−2. Another necessary parameter obtained from ERA5 was the average CF, which was 0.6. We only considered the data corresponding to daylight hours for I and CF; therefore, values during night-time (I = 0 W m−2) do not affect the average. For the calculation of the mean snow crystal radius (r), the grains that made up the intermediate and basal low-density layers, with radii ranging between 2 and 3.5 mm, were excluded, as these crystals result from post-depositional processes, a fact discussed in the following section. Therefore, we considered the average ice crystal radius was 1.5 mm. It is crucial to mention that we made only approximations to obtain a general overview of the potential consequences of the BC concentrations found in the snowpack. The results provide an initial idea, but it is essential to interpret them with caution since they are preliminary estimates rather than a comprehensive and precise assessment.
To perform the quantitative analysis, we took the average shortwave solar irradiance as a reference value (106.68 W m−2). Concerning the snowpack BC concentration, our focus is not on variations along each snow pit, as we did with the isotopic profiles, but rather on the overall bulk concentration of both combined. The BC load in the Godwin-Austen Glacier snowpack ranged from 2 to 44 ng g−1, with an average of 12 ng g−1. Given the high capacity of BC particles to absorb shortwave solar irradiance, they induce the snowpack to absorb additional energy ranging from 7.8 to 50.35 MJ m−2, with an average of 22.1 MJ m−2. In turn, this increase in energy absorption would reduce the snowpack surface albedo by between 0.0068% and 0.044%, with an average of 0.019%. The extra energy absorption would lead to snowpack warming and, consequently, a positive radiative forcing ranging from +0.73 to +4.73 W m−2, with an average of +2.07 W m−2. These data allow us to estimate the volume of snow that melted sooner per unit area, ranging from 23.3 to 150.7 mm w.e. (or kg m−2), with an average of 66 mm w.e. Regarding the BC load in the high-altitude K2 camps, it could have produced a volume of melted snow of 45.6 mm w.e. in Camp 1 and 109.9 mm w.e. in Camp 2 (Table 1).
Table 1. Black carbon (BC) concentration in the snow pits (P1 and P2) and K2 high-altitude camps (camp 1 and camp 2), and its effects on snow considering the average shortwave irradiance (106.68 w m−2)
ΔA: albedo reduction, RF: radiative forcing, Ex: extra energy absorbed by the snowpack, W: snow melting sooner.
5. Discussion
5.1. Snowpack temporal scale
The correct interpretation of BC concentration data requires determining the time interval of the snowpack accumulation, as the implications of BC levels in a multi-year snowpack differ from those in a seasonal one. To address this, we focused on water stable isotopes, examining the correlation between δ1⁸O values and temperature, as well as the relationship between moisture sources and d-excess levels. However, it should be considered that the original isotopic composition of the snowpack may be altered from its state at the time of deposition.
The low-density surface layers are likely the result of recent snowfall (Colbeck and others, Reference Colbeck1990; Fierz and others, Reference Fierz2009). The crystals with well-defined faces and edges in the middle and basal low-density layers are characteristic of depth hoar (DH; Seligman, Reference Seligman1936). The DH crystals typically form under high-temperature gradients (Akitaya, Reference Akitaya1974; Colbeck, Reference Colbeck1982; Figs 2 and 3), as observed in the Godwin-Austen snowpack (∼6ºC m−1). The development of DH crystals was more pronounced in P1, possibly due to a greater temperature gradient relative to the snowpack thickness, while higher gradients in the upper 60–80 cm (up to ∼18ºC m−1; Flin and Brzoska, Reference Flin and Brzoska2008) may explain the middle DH layers (Fig. S3). High-temperature gradients, such as those associated with DH, stimulate molecular diffusion within the pore space, which can potentially alter the isotopic signal in the snowpack. Similarly, diffusion at the snowpack surface may occur due to interactions with the atmosphere (Ekaykin and others, Reference Ekaykin, Lipenkov, Barkov, Petit and Masson-Delmotte2002; Konishchev and others, Reference Konishchev, Golubev and Sokratov2003). In this vein, no δ1⁸O variations were detected in relation to the DH layers. However, the higher d-excess values in the DH basal layer of P1 likely suggest intense diffusion, presumably related to the greater development of DH crystals. In the case of surface diffusion, the high precipitation during the 2018–19 winter in the Karakoram region likely preserved the original isotopic signal by reducing atmospheric exchange (IMD, 2019; National Weather Forecasting Centre, 2019). This was further supported by ERA5 data, which identified 2018–19 as one of the three wettest seasons of the past decade, with an average daily precipitation of 2.20 mm w.e. (Fig. S4).
Given the minimal alteration of isotopic values, ERA5 reanalysis (Hersbach, Reference Hersbach2016) was employed to estimate temperatures for the entire accumulation season (October 2018–April 2019) and correlate them with δ1⁸O levels in the snowpack. Despite its lower accuracy and resolution compared to the AWS data (available for only 24 days starting on 3 February 2019), ERA5 effectively replicates the temperature patterns recorded by the AWS, making it a reliable tool for assessing long-term trends, even with an underestimation of ∼5ºC. For the correlation, data from P2 (200 cm) was prioritized over P1 (120 cm) due to its location, which provides greater relief, allowing for more vertical accumulation and better protection from erosion. Consequently, the snowpack thickness in P2 was equated to the accumulation season, with the snow–ice interface corresponding to approximately October 2018 and the snowpack surface marking late February 2019, when the snow pits were excavated. The δ18O variations in P2 closely followed the ERA5 atmospheric temperature patterns (Fig. 5). Temperatures peaked at −14ºC at the start of the season (October 2018), dropped to −32.7ºC in late December 2018, and rose again to −25ºC by late February 2019. Similarly, δ18O peaked at −17.8‰ in October 2018, reached a minimum of −30.5‰ in early January 2019 and increased to −21.7‰ by the end of February 2019.
Figure 5. Comparison of δ18O variation along P2 (black line) with the daily mean temperature values at 2 m above the ground obtained from ERA5 (solid grey line). the temperature values acquired in situ by the Base Camp automatic weather station (AWS) are depicted by the dashed grey line.
To confirm the relative scale of the snowpack suggested by the δ1⁸O-temperature correspondence, similar correlations can be established between the d-excess levels along P2 and the origin of the predominant trajectories during the accumulation period (Fig. 6). During October and November 2018, the high average d-excess value (16‰) coincided with a strong influence from westward trajectories (Mediterranean and Middle East). Local moisture recycling processes, such as re-evaporation or evapotranspiration, may have also contributed to these elevated d-excess values (Froehlich and others, Reference Froehlich, Kralik, Papesch, Rank, Scheifinger and Stichler2008; Aemisegger and others, Reference Aemisegger, Pfahl, Sodemann, Lehner, Seneviratne and Wernli2014). The values decreased to ∼5‰ in late December 2018 and early January 2019, likely due to the dominance of Black Sea and Caspian Sea trajectories at 500 m a.g.l., typically associated with low d-excess values in winter (Juhlke and others, Reference Juhlke2019). By February 2019, d-excess rose again to 10‰, probably due to the return of westward trajectories. These variations can be divided into two groups: higher d-excess values in October–November 2018 and late February 2019 (10.2‰), and lower values from December 2018 to early February 2019 (7.5‰). The corresponding Local Meteoric Water Lines for these groups reveal different intersections with the y-axis (d-excess) compared to the GMWL, confirming the distinct origins of the air masses (Fig. S5). Therefore, this correspondence further reinforces the seasonal nature of the snowpack.
Figure 6. Frequency diagram illustrating back-trajectories contributions to each domain at both altitudeS: 500 and 1500 m a.g.l. isotopic variations of both snow pits (P1 and P2) are depicted behind the bar chart. domain 1 (L: local), domain 2 (NCA: North Central Asia and S: siberia), domain 3 (CS: Caspian Sea and BS: Black Sea), domain 4 (ME: Middle East and MS: Mediterranean Sea), domain 5 (PG: Persian Gulf and RS: Red Sea) and domain 6 (NI: Northern India and AS: Arabian Sea).
Finally, to establish a relationship between P1 and P2, d-excess variations were compared, as their patterns were more distinct, although δ1⁸O variations were also considered. Four comparable sections (S1–S4) were identified between both profiles (Fig. 7). The basal section (S1) showed the highest δ1⁸O values, with similar d-excess patterns, except for a peak in P1, likely due to intense molecular diffusion. In S2, maximum d-excess values were comparable. In S3, δ1⁸O reached a minimum, with corresponding lower d-excess values in both profiles. Finally, the last section (S4) displayed similar isotopic values in both records. Additionally, three hiatuses (H1–H3) were identified in P1 between the sections, possibly correlating with periods of lower precipitation, which may have increased exposure to erosive factors.
Figure 7. (Top) potential correlations between P1 and P2 based on their Isotopic variations: δ18O (solid lines) and d-excess (dashed lines). (1) P2 δ18O max. (−17.8‰), (2) P1 δ18O max. (−16.3‰), (3) P2 δ18O min. (−30.5‰), (4) P1 δ18O min. (−30.4‰), (5) P1 d-excess max. (19‰), (6) P2 d-excess max. (18‰), (7) P1 d-excess 2nd max. (18.6‰), (8) P1 d-excess min. (6.1‰), (9) P2 d-excess min. (4.4‰). correlation sections are depicted as grey bands with the letter S, while sections with missing records in P1 appear as white bands with the letter h (hiatus). (Bottom) seasonal precipitation data (dotted line) from ERA5. Hiatuses (H) may correspond to periods characterised by less frequent and/or less intense total precipitation events.
Previous studies have found high winter d-excess values in Central Asia due to moisture from the Eastern Mediterranean (Hussain and others, Reference Hussain2015; Liu and others, Reference Liu2015; Juhlke and others, Reference Juhlke2019). However, the Godwin-Austen snowpack showed low d-excess values during the 2018–19 winter, likely influenced by Black and Caspian Sea trajectories. This highlights the complexity of atmospheric dynamics in the western HMA, particularly during the positive precipitation anomaly of that particular season (Fig. S6; Dimri and others, Reference Dimri2015; Hunt and others, Reference Hunt, Curio, Turner and Schiemann2018; Attada and others Reference Attada2022; Nischal and others, Reference Nischal and Hunt2022).
5.2. Implications of BC-enhanced snowmelt
The clear difference in the BC concentration between Zones 1 (Z1) and 2 (Z2) in both snow pits is consistent with previous studies that reported higher concentrations in the bottom of the snowpack (e.g., Thakur and others, Reference Thakur2021). This distribution can be explained by the partial melting of snow layers formed during snowfalls at the beginning of the accumulation season (autumn), which can potentially concentrate BC particles. This occurs due to the still high solar radiation and the fact that the Godwin-Austen Glacier is debris-covered, that reduces the albedo of its surface, leading to greater absorption of solar energy, increased surface heating and ultimately accelerating the melting of the initial snow layers. As a result, the BC average concentration in Z2 (27 ng g−1 in P1 and 18 ng g−1 in P2) could be attributed to these reconcentration processes. Conversely, the Z1 average concentration in both snow pits (8 ng g−1), with lower variability compared to Z2, is likely less influenced by reconcentration processes, providing a more accurate depiction of the local average snowpack BC concentration during the 2018–19 winter season and its related effects. Consequently, the BC concentration in Z1 would produce 53 mm w.e. of meltwater.
To give a sense of the scale of meltwater generated by the total average BC during the 2019 melting season (66 mm w.e.), we can reference the average precipitation during the study season (from October to February) over the previous 10 years (2009–19), which was 283 mm w.e. In this context, the meltwater generated could potentially represent ∼23%, while considering only Z1, this percentage would reduce to ∼18%. However, caution is necessary when interpreting these values due to inherent uncertainties in the data, particularly regarding the accuracy of precipitation estimates, which rely on interpolations from the ERA5 atmospheric reanalysis. Additionally, in remote mountainous regions with limited instrumental data, such as the Karakoram, ERA5 often struggles to capture local precipitation dynamics accurately. Furthermore, the average BC concentration is based on measurements from only two isolated points on the glacier (P1 and P2), which limits the representativeness of the data. Regarding the BC load in surface snow at high-altitude camps, we cannot quantify the effects in the same way as we did for the snowpack, since these are even more localized measurements, as they are limited to surface snow samples. However, their analysis provides insights into how the presence of this pollutant impacts regions at high altitudes and how the activities conducted at these camps contribute to environmental contamination. In this context, the presence of this pollutant could produce 45.6 mm w.e. of meltwater at Camp 1, increasing to 109.9 mm w.e. at Camp 2.
5.3. Potential local and regional BC sources
The Baltoro and Godwin-Austen Glaciers serve as the main route to the base camps of K2 and Broad Peak, involving a seven-day trek over ice and debris (see Fig. 1b). Expeditions set up temporary camps near permanent army bases, both reliant on kerosene for cooking and heating. During our expedition, the AWS at K2BC recorded predominant strong winds (>10 m s−1) from the west (Saboy Glacier), with some from the northwest and fewer from the south, suggesting limited pollutant dispersion from K2BC to the snow pit sites. The minimal impact from nearby camps like Concordia, located southward, can be attributed to the infrequency of southern winds (Fig. S7). This indicates that the BC average concentration in the snow pits (12 ng g−1) may be representative of the overall concentration in the Godwin-Austen snowpack for the 2018–19 season, with contributions from distant and local sources. Although no wind data are available prior to 3 February 2019, similar patterns were observed starting from early January 2019.
To evaluate regional BC sources, we used monthly average atmospheric concentrations from the MERRA-2 atmospheric reanalysis (Fig. S8). Between September 2018 and December 2019, BC concentrations increased across South Asia, followed by a decrease from January to February 2019. This trend was particularly notable in the Indo-Gangetic Plain (IGP), where concentrations reached up to 8 μg m−3 in December 2018. This increase was influenced by high population density, industrialization, higher heating demand due to lower temperatures and a lack of precipitation, which reduced BC removal through wet deposition. Additionally, the surrounding mountain ranges, such as the Himalayas to the north and the Safēd Kōh Range to the west, acted as barriers, limiting the dispersion of pollutants. Furthermore, from September 2018 to February 2019, other regions of Central Asia and the Middle East exhibited similar patterns of BC atmospheric concentration, though to a lesser extent.
By combining BC atmospheric concentration with back-trajectories, we can estimate the contribution of long-distance BC sources to the snowpack during the accumulation period (Fig. S8). In September and October 2018, the main sources at 500 and 1500 m a.g.l. were the northern Indus River basin lowlands, with secondary contributions of South Asia, eastern IGP, the Middle East, Central Asia and Western China. In November 2018, the Middle East was the major contributor at 500 m a.g.l., followed by Central Asia and the Near East. In December 2018, the westerly jet stream increased the contribution of Central Asia, Eastern Europe, the Middle East and the Near East, while the northern Indus basin lost influence. The trajectories at 1500 m a.g.l. from November 2018 to February 2019 show a similar contribution to that of December 2018 at 500 m a.g.l. Continuing at 500 m a.g.l., in January 2019, the westerlies persisted, but some clusters may have carried BC from the central Gangetic Plain. In February 2019, contributions mainly came from the West, with influence from the Persian Gulf region, Middle East, and, to a lesser extent, from Central Asia. The connection between BC distant sources and back-trajectories is an estimation with inherent limitations: BC concentration in the snowpack does not reflect accurately the predominant long-distance sources due to factors like dry or wet deposition (outputs) and inputs from regions along the transport path.
The trajectories during the accumulation season in K2 Camps 1 and 2 matched those for the snow pits, likely due to the nearby sampling points and the resolution of meteorological data used in HYSPLIT. Even with a higher resolution, the general patterns would not have important variations, so we consider this result valid. High-altitude snow samples were taken at the end of February 2019, coinciding with snow pit sampling. The BC values at K2 Camp 1 (6076 m a.s.l.; 7 ng g−1) were similar to surface snow in the snow pits, constituting a reliable representation of BC concentration in the surface snow during late February 2019. In contrast, Camp 2 (6650 m a.s.l.) had much higher concentrations (26 ng g−1), likely due to activities from successive expeditions, indicating their impact on the high-mountain environment.
5.4. Godwin-Austen Glacier BC in HMA: A comparison
The HMA region consists of several sub-regions, each influenced by distinct BC sources and atmospheric patterns (Table 2; Fig. 8). In comparing results across studies, it is important to consider that differences in sampling methodologies and post-depositional processes can alter BC concentrations and their effects. Snow pit excavation studies generally follow one of two approaches: either calculating the average BC concentration for the entire snowpack and treating it as a representative value, or conducting a more detailed layer-by-layer analysis, which requires higher sampling resolution. Additionally, surface snow sample analyses show that BC concentrations can fluctuate based on the time between deposition and sampling. Some studies also distinguish between precipitation snow and freshly fallen snow, as BC concentrations can change even over short periods. Surface snow modification mainly results from phase changes, which are primarily influenced by temperature variations linked to altitude and seasonality. As a result, in lower-elevation areas and warmer seasons, surface snow undergoes more significant modification, while at higher altitudes and during colder seasons, these changes are less pronounced. Selecting the sampling point appropriately is crucial, taking into account local BC sources, their proximity, and potential influences (e.g., industrial areas, urban areas or roads), as these factors can lead to an overestimation of background BC levels.
Figure 8. Locations and distribution of studies focused on BC in snow across HMA. Each code (number) corresponds to a study listed in Table 2.
Table 2. Summary of BC studies in snow across the HMA region
Studies are organized by subregion and listed in ascending order of average BC concentration. Geographic locations are indicated by corresponding codes in Fig. 8. The ‘Study area’ indicates the type of location: glacier (Gl.), region (Reg.), locality (Loc.) or snow station (S.St.). Seasons are abbreviated as Sum. (summer), Spr. (spring), Aut. (autumn), Win. (winter) and M.S. (melting season). Sample types include Ss (surface snow), SsA (aged surface snow, > 24 h post-deposition), SsF (fresh surface snow, <24 h post-deposition), Swp (snow pit) and P (precipitation). The letter ‘n’ indicates the number of samples, and ‘depth’ refers to the maximum snow pit depth (n/a for surface snow samples). The BC concentrations are given in ng g-1, with minimum, maximum and average values provided. A dash (-) indicates unavailable or unspecified data. The ‘Ref.’ column contains the reference of each study.
Due to the heterogeneity of sampling conditions, in terms of seasonality, altitude and methods, establishing a direct comparison between our data and previous studies would be inappropriate, even if some values are similar. To address this issue, it is crucial to define clear standard criteria, such as the maximum depth to consider a sample as surface snow. The lack of well-defined criteria in the literature can lead to ambiguous interpretations and hinder data comparability. Besides, following and reporting sampling protocols is essential to maintain data integrity, including using sterilized materials and personal protective equipment, especially in isolated areas with low concentrations, such as our study site. Nevertheless, conducting a comparative analysis of existing BC studies in the HMA region would be a valuable exercise to situate our study within a broader context (Table 2).
Compared to previous snow pit excavation studies, the work conducted on the Godwin-Austen Glacier is notable as the only study performed during the winter season, since this type of fieldwork is often carried out in the summer. Two factors may explain why the average BC concentration in our snow pits (12 ng g−1) is significantly lower than the average concentration reported in previous snow pit-based studies (127 ng g−1). First, the concentration in summer studies is likely influenced by BC reconcentration processes resulting from snowmelt. Second, the unusually high snow accumulation during the winter of 2018–19 may have diluted BC levels. These factors also account for the fact that the snow pits in this study recorded both the lowest minimum BC concentration (2 ng g−1) and one of the lowest maximum concentrations (44 ng g−1) in the dataset. Our snow pits reached a depth of 200 cm, only comparable to the study of Kaspari and others (Reference Kaspari, Painter, Gysel, Skiles and Schwikowski2014) on the Mera Glacier (Himalayas), and included the highest number of samples (43). In terms of altitude, the snow pits were excavated at a relatively lower average elevation (∼4930 m a.s.l.) compared to the average elevation of ∼5330 m a.s.l. observed in similar studies. This difference is likely due to the more favourable summer conditions for snow pit excavation. Notably, this investigation is one of the few conducted in the Karakoram region, alongside the studies by Thakur and others (Reference Thakur2021) on the Kardung and Phuche glaciers.
Surface snow studies are common throughout the HMA, as they are more time-efficient and less labour-intensive compared to snow pit excavations. These samples are typically collected at an average altitude of ∼4740 m a.s.l., with the samples from Camp 1 (6076 m a.s.l.) and Camp 2 (6650 m a.s.l.) standing out due to their significantly high altitudes. They also exhibit greater seasonal heterogeneity, with samples collected during spring, autumn and winter. The average BC concentration in the dataset of 193 ng g−1 is significantly higher than the concentrations found in our samples, which range from 7 to 26 ng g−1. The BC concentrations in the winter surface snow samples are expected to be low due to higher accumulation rates and reduced melting, but the exceptionally high concentrations reported by Zhang and others (Reference Zhang2017) and Gul and others (Reference Gul2018) skew the winter average concentrations of the dataset to 1112 ng g−1. Excluding these outliers, the average drops significantly to 38 ng g−1. However, it is important to note that our study has one of the lowest levels of sample representativity for surface snow compared to previous studies in the region.
6. Conclusions
This work presents the first attempt to quantify the BC concentration in the seasonal snowpack over the Godwin-Austen Glacier, Central Karakoram Range, Pakistan. To achieve this, we used snow samples extracted from two snow pits (P1 and P2) excavated south of the K2BC, along with surface snow samples from K2 Camps 1 and 2.
The alignment between the δ1⁸O profile and temperature records from ERA5 reanalysis and the K2BC AWS, along with the correlation of d-excess variations with potential source trajectories during the accumulation season, suggests that the studied snowpack has a seasonal character, having accumulated at the beginning of the accumulation season around October 2018 and disappeared during the 2019 ablation season. Determining the timing of snowpack formation is essential for accurately interpreting BC data.
The back-trajectory analysis reveals shifts in dominant regional BC sources during the accumulation season, with contributions from the northern Indus basin in early autumn, followed by increasing influence from the Middle East, Central Asia and Eastern Europe through winter. Local sources, including army camps and base camps, such as the K2BC, also contribute significantly to BC levels in the area. Together, these sources account for the BC concentrations observed in the Godwin-Austen Glacier snowpack. Based on snow pit samples, the average BC concentration is 12 ng g−1. The similar surface concentrations and close sampling dates of the snow pit and surface snow samples at K2 Camp 1 (∼7 ng g−1) suggest that this concentration may represent the prevailing BC level in the area in late February 2019. The elevated BC average concentration in the lower section (Z2) of the snowpack (ranging from 18 to 27 ng g−1) could result from early-season snowmelt and subsequent particle reconcentration, indicating that the average concentration measured in Z1 (8 ng g−1) may be more representative. Additionally, the high concentration observed at Camp 2 (26 ng g−1) likely reflects the impact of mountaineering activities and their associated environmental effects.
Considering the average solar irradiance (106.68 W m−2) and mean BC concentration, BC-induced albedo reduction caused 66 mm w.e. of snow to melt sooner during the 2019 melting season. To contextualize these findings, this would represent ∼23% of the regional average precipitation volume measured over the past 10 years for the same reference period, from October to February. This percentage decreases to 18% when considering the lower concentration in Z1. However, caution is necessary due to data limitations and reliance on a small number of sampling points.
When comparing the results with previous studies, it is important to consider that variability in BC concentrations in the snow across the HMA region is influenced by diverse sources and atmospheric patterns. Besides, differences in sampling methodologies and post-depositional processes complicate direct comparisons of results. Consequently, to ensure reliable data comparability, standardized sampling criteria and protocols, along with consideration of local BC sources, are essential.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/jog.2024.111.
Acknowledgements
This research was supported by the Spanish Government through the AEI-funded MicroFAME project (Ref. PID2022-140975NB-I00), as well as through Ref. MDM-415 2017-0714 (incl. BALELUR-Sherpa project) and Ref. CEX2021-001201-M funded by MICIN/AEI/10.13039/501100011033. This study has been also partly supported by the Arctic Challenge for Sustainability (ArCS) Project (Program Grant Number: JPMXD130000000), the Arctic Challenge for Sustainability II (ArCS II) Project (Program Grant Number: JPMXD1420318865), and the Environment Research and Technology Development Funds (JPMEERF20172003, JPMEERF20202003 and JPMEERF20232001) of the Environmental Restoration and Conservation Agency of Japan. We are indebted to Alex Txikon (https://www.alextxikon.com) and his team of mountaineers for their logistic support and willingness to collect snow samples from the K2 camps. Finally, we extend our gratitude to Deniz Bozkurt from the University of Valparaiso for his valuable contributions.
Competing interests
The authors declare no competing interests for this paper.
Open access
