High altitude accumulation and preserved climate information in the western Pamir, observations from the Fedchenko Glacier accumulation basin

Abstract

The accumulation region of Fedchenko Glacier represents an extensive snow reservoir in the Pamir Mountains feeding the longest glacier in Central Asia. Observed elevation changes indicate a continuous ice loss in the ablation region of Fedchenko Glacier since 1928, while the mass balance of the accumulation region is largely unknown. In this study, we show that accumulation varies considerably in the main accumulation basin, with accumulation rates up to 2400 mm w.e. a⁻¹ in the West, decreasing to <1000 mm w.e. a⁻¹ in the center, although the elevation difference is <200 m. The combination of snow/firn samples and ground-penetrating radar profiles suggests that this accumulation pattern is persistent during the recent past. The recent accumulation history is reconstructed from internal radar reflectors using a firn densification model and shows strong interannual variations, but near constant mean values since 2002. Modeling of trajectories, based on accumulation and glacier geometry, results in an estimate of the depth/age relation close to the main divide. This region provides one of the most suitable locations for retrieving climate information with temporal high resolution for the last millennium, with a potential to cover most of the Holocene in less detail.

1. Introduction

The glaciers in High Mountain Asia (HMA) play an important role for the regional cryosphere and hydrology. They contain ~9600 km³ of ice compared to, for example, 117 km³ stored in the glaciers in the Alps (Huss and Farinotti, 2012). The glaciers of Central Asia alone represent more than 5000 km³ of ice and many of these glaciers are considered to be almost in balance (e.g. Gardelle and others, 2013; Brun and others, 2017). Due to their size and extensive glacier tongues, their meltwater contribution to the rivers is important and reaches ~13% for all glaciers in the Amu Darya watershed (Luo and others, 2018). However, ablation conditions are strongly influenced by debris cover and connected surface features such as supra-glacial ponds and ice cliffs (e.g. Juen and others, 2014; Thompson and others, 2016; Irvine-Fynn and others, 2017). The large glaciers in this region also play an important role for the regional climate, because of high albedo across extensive accumulation basins and thus strong reflection of incoming radiation. Apart from the Karakoram (e.g. Biafo-Hispar Glaciers, Baltoro Glacier, Siachen Glacier) and the Tian Shan (e.g. Inylchek Glaciers), one of the largest glacier systems in HMA is found in the Pamir Mountains, with Fedchenko Glacier in its center. Investigations of accumulation conditions and glacier geometry are the key to understand the mass-balance situation of this dominating glacier in the Pamir Mountains. In addition, these basins might also represent important climate archives.

Quantification of ice ablation improved with advances in the understanding of the energy transfer through the supra-glacial debris cover (Evatt and others, 2015) and the development of more sophisticated mass-balance models (Groos and others, 2017). In contrast, high elevation precipitation and thus accumulation usually shows very low correlation to data from nearby valley stations (Immerzeel and others, 2015). Thus, the least known parameter for establishing mass budget calculations of the high Asian glaciers is the high elevation accumulation (e.g. Cannon and others, 2017). This parameter cannot be retrieved simply from remote-sensing information, as dry snow is transparent for active microwave sensors, while elevation changes from sensors in the visible spectrum do not work well on surfaces without contrasting features (Paul and others, 2017). Even though recent progress in snow depth retrieval from active microwave sensors indicates a strong potential of such sensors (Lievens and 16 others, 2019), the low spatial resolution and the crude universal calibration prohibit a useful application in high mountain ranges with rough topography. Also, downscaling from meteorological models does not work well in such regions, due to complicated local conditions (Cannon and others, 2017). Climate models themselves are similarly unable to provide reliable accumulation rates across regions with large elevation differences in the absence of ground-based calibration information (Reggiani and others, 2016). Ground truth information provides important input to improve such model results. In this paper, we present direct measurements of snow accumulation from snow pits and radar which we use to constrain the recent accumulation history.

At high altitudes where other paleoclimate records are scarce, ice cores from non-temperate mountain glaciers can provide important records of past environmental and climate conditions (Bohleber, 2019). Some of the accumulation basins in HMA might represent highly valuable climate archives, because they extend to high elevations and the ice thickness is probably rather large. These conditions prevent the occurrence of layer-disrupting surface melt and enable a long and continuous history of ice deposition. An ice core from Inylchek Glacier, for example, provided recent climate and trace element records for most of the 20th century (Grigholm and others, 2016). Some of these cores even cover a few hundred years (e.g. Aizen and others, 2004) and there are indications for the existence of much older ice (Thompson and others, 1997). Shallower cores in different regions of High Asia provided climatic and atmospheric information for the last decades (e.g. Aizen and others, 2006; Duan and others, 2015). However, only some short firn/ice cores, covering the last centuries are available in the wider region of the Pamir mountains (e.g. Wang and others, 2015; Yao and others, 2017), while a long-term climate record from an ice core is still missing. The geometry and elevation of the accumulation basin of Fedchenko Glacier suggests it could be one of the best sites in the Pamir Mountains to recover a climate record of the entire Holocene and possible beyond. Investigations by Aizen and others (2009) already showed that snow and firn layering is preserved in the uppermost accumulation area of this glacier. Accordingly, a regional long-term paleoclimate reconstruction from Fedchenko ice cores promises an important complement to other proxy records, for example, from tree-rings (Esper and others, 2002).

This is only one of the reasons for investigating these high elevation basins with respect to accumulation conditions, ice distribution and mass flux. Also, the accumulation conditions in the main accumulation basin may be the key to understand the close to balanced glacier mass budgets in this region (Gardelle and others, 2013; Brun and others, 2017). As a basis for these investigations, we collected information about ice thickness, surface topography, accumulation distribution and meteorology. We then calculated the basal geometry, the recent accumulation history, and evaluated the potential for retrieving an ice core with a maximum quality and length of paleoclimate record.

2. Region of investigation

The Pamir Mountains contain more than 11 000 glaciers (Mölg and others, 2018). The most extensive glacier system in the Pamir Mountains is Fedchenko Glacier with its tributaries (e.g. Nalifkin Glacier) and neighboring glaciers (e.g. Grum Grijmailo glacier). A major tributary in the lower ablation region is the surge-type Bivachniy Glacier (Wendt and others, 2017). Fedchenko Glacier covers an area of 745 km² (Lambrecht and others, 2014) and extends over a total length of 72 km (Fig. 1). The accumulation area ratio is ~0.64 (Lambrecht and others, 2014) based on the observations since 2000 and the accumulation area is distributed across several extensive basins. The largest of these basins (total glacier area of 68 km²) is the origin of the main glacier trunk and extends from ~5100 m elevation up to Jasgulem Pass at 5300 m elevation and into some higher sub-basins. We concentrated our investigations on this basin because it contains the least disturbed layer sequence, both due to its size and smooth surface topography. However, some minor areas are influenced by avalanching from the main peaks in the vicinity of Peak Istiqlol (formerly Pik Revolution).

Fig. 1. The location of Fedchenko Glacier in Central Asia is displayed in panel (a). Panel (b) shows the glacier boundaries and the location of the main accumulation basin (red rectangle). Details of the main accumulation basin with ground-penetrating radar profiles (purple: high frequency, blue: low frequency) and the location of the snow pits are presented in panel (c). The colored, thick segments represent the GPR profiles presented in Figure 4. The background image of (b) and (c) is based on a Landsat 8 scene from 2013.

3. Data and methods

We collected data during three summer seasons in 2009, 2015 and 2016, concentrating on the recent accumulation distribution, near surface stratigraphy, topography and ice thickness. The information on accumulation rates, snow and firn densities from snow pits and shallow cores, combined with ground-penetrating radar (GPR) information allows us to extend the information retrieval to the deepest detectable firn pack reflectors. This information is also spatially interpolated to gain insight into the variability of accumulation across the upper basin. Low-frequency GPR data provide information about ice thickness in the basin. The potential temporal archive in the entire vertical glacier column is then estimated by analyzing a simple trajectory model, forced with the above input parameters.

3.1 Snow and firn investigations

We excavated three snow pits of ~2 m depth at elevations of 5200 m (S1) in the lower center of the basin, 5300 m (S2) at Jasgulem pass and 5360 m (S3) on the saddle to the North of the pass in 2015, in order to obtain information about recent accumulation rates at selected sites across the basin (Fig. 1). The combination of depth and sample densities results in a vertical accumulation distribution. Conductivity measurements and stable isotope ratios are used to determine the seasonal distribution and also annual accumulation rates. At two sites (S1, S2), we cored from the bottom of the pits to extend the record to 5 m and more than 6 m, respectively. The lowest snow pit (S1) was repeated in 2016, down to the summer horizon of 2015. A summary of the site elevations and the reached depths is given in Table 1.

Table 1. Field sites of direct snow and firn investigations

Presented are the site elevations, the pit depths and the additional core lengths. The total depth in 2015 is 6.66 m for S1 and 5.38 m for S2. Snow pit probing and drilling was done on August 12 (S1), August 16 (S2) and August 17 (S3) in 2015, and August 15 (S1) in 2016, respectively.

Density was measured in situ at snow and firn samples from the pit and the core by weighing defined sample volumes, with an uncertainty of ~5%. The uncertainty results from the 5 g-accuracy of the used scale and the estimated uncertainty of the probe volume, which is never filled exactly to 100%. The occurrence of ice layers was documented by visual/manual stratigraphy, where ice layers usually are clearly discernible from snow and firn. The vertical snow and firn sequence of each pit was sampled at 10 cm intervals, where the interval resulted as a compromise between high-resolution sampling and the restricted amount of sample material which could be transported. The samples were filled in previously cleaned vials, melted, and transported to the Institute of Environmental Physics, Heidelberg (IUP-HD), for further analysis. Meltwater conductivity was measured by manually injecting a WTW LF500 conductivity probe into each vial (measurement uncertainty ~0.5%). Following the methods of a simultaneous study (Bohleber and others, 2018), co-isotopic measurements of stable water isotopes (δ ¹⁸O and Deuterium) were performed at IUP-HD using a Picarro cavity ring down spectrometer. Typical measurement uncertainties range within 0.1 and 0.4 permil for δ ¹⁸O and Deuterium, respectively.

3.2 Ground-penetrating radar measurements

We used GPR for two different purposes: collecting information about the total ice thickness and investigating the layering of the firn pack. Ice thickness was determined from almost 30 km of low-frequency radar data in the upper accumulation basin. A small pulse transmitter (Narod, Icefield Instruments Inc., Whitehorse, YT Canada) combined with resistively loaded dipole antennas (center frequency 5 MHz) acted as energy source, while the reflected electromagnetic waves were collected by identical dipole antennas and a digital oscilloscope with a sampling rate of 4 ns. Pulling the system across the snow surface resulted in a trace sampling frequency of about one trace every 2–4 m. Combined errors from the radar measurements and geolocation were calculated according to the method described in Lapazaran and others (2016) and result to an error range of 8–18 m, depending on the total ice thickness. A first set of ice thickness profiles was collected in 2009 covering ~4.8 km in the upper accumulation area (Lambrecht and others, 2014). A more extensive survey was carried out in 2016 resulting in an additional 24.8 km of profiles in the same region.

The firn layering was determined by a commercial GPR system using a dual-frequency antenna (IDS K2 Fast Wave, 200 and 600 MHz). About 36 km of radar profiles were collected in July and August 2015. The penetration depth of both antennas is considerably less than the ice thickness and thus only internal structures of the upper 30–50 m of the glacier are revealed with a vertical resolution of 0.08–0.12 m.

All radar data were processed with a commercial software package (ReflexW, Sandmeier geophysical research). Processing steps include start time shift, dewow, depth-dependent gain function and frequency filtering for enhancing the signals from the internal reflectors and allow the correct travel time detection. Reflectors were picked semi-manually, using the rising flank of the reflection signal and correcting the phase follower at obvious deviations from the reflector. The bedrock reflection in the low-frequency data is clearly identifiable, because the ice above the bed is basically free from reflections. All snow pits in the lower accumulation basin since 2009 have had ice lenses and accumulated dust concentrations from snow melt and sublimation associated with the end of summer layer. Those layers produce strong reflectors in the high-frequency data.

Continuous position mapping with GNSS receivers provides the location information for the radar profiles, based on precision point positioning post-processing, with a horizontal position error of <1 m.

3.3 Vertical accumulation distribution

In order to convert the firn pack stratigraphy from the radar data into accumulation information, it is necessary to convert the layer thicknesses to water equivalent units using a depth-dependent density distribution. This enables the extension of the accumulation history to the end-of-season layers identified in the radar measurements. Densification can be estimated by the ambient air temperature and accumulation rates for dry snow conditions (Herron and Langway, 1980). In the presence of meltwater, densification occurs faster and refreezing alters the vertical density structure. By separating an annual layer into a firn fraction and an ice fraction, Reeh and others (2005) and Reeh (2008) accounted for surface melting and refreezing in the snow/firn pack and thus the presence of ice lenses. The density of the firn fraction of an annual layer ρ _s then can be determined by

(1)$$\displaystyle{{\rho _{\rm i}-\rho _{\rm s}} \over {\rho _{\rm i}-\rho _{{\rm s0}}}}\left( {\displaystyle{{\rho_{\rm s}} \over {\rho_{{\rm s0}}}}} \right)^{({\rm SIR}/b)-1} = exp\left( {-\displaystyle{c \over b}d} \right).$$

here, ρ _i is the ice density, ρ _s0 the snow/firn density of the uppermost layer, SIR the fraction of refrozen surface melt, b the specific net balance, d the depth below the surface and c a parameterization of the pressure-determined densification process which is dependent on air temperature. For further details, especially the choice of c, we refer to Reeh and others (2008). The optimal fit between layer density and specific net balance can only be determined by an iterative procedure, where accumulation rates are adjusted until the resulting layer depths fit the radar stratigraphy.

This empirical model assumes that the meltwater from one season will refreeze within the corresponding annual layer. This is probably only the case for regions with small melt amounts (e.g. upper percolation zone), such as those observed in the upper part of Fedchenko Glacier, where ice layers are at maximum a few centimeters thick in the snow pits and firn cores. We use this simple model, because only very limited information is available about the meteorological conditions in the upper Fedchenko basin, which prohibits the application of more sophisticated approaches. The model requires mean monthly air temperatures, mean annual accumulation rate and the surface snow density as input for firn densification with depth. The amount of melt is estimated by a degree-day approach, which, again, is sufficient to provide a good approximation of the ice fraction in the annual layer (Reeh, 1991; Reeh and others, 2005). The required temperature input for calculating annual melt is derived from the High Asia Reanalysis simulations (HAR, Maussion and others, 2014), by applying a mean vertical lapse rate for the elevation difference between the level of the model data and the glacier surface. For this purpose, we used data from an automatic weather station, which was running on Fedchenko Glacier at 4950 m elevation from August 2009 until September 2012. This station was installed by colleagues from Nagoya University, Japan with the aid of some of the authors in the framework of the CADIP program (http://www2.umaine.edu/climatechange/Research/projects/CADIP.html). No other correction has been applied to the HAR temperature data, apart from this lapse rate calibration.

In 2015, the density structure was measured in the snow pits for the first 5–6 m. Because we do not know a priori the exact vertical radar wave velocity distribution, it is necessary to iteratively adapt the wave velocity/depth function and firn densification parameters until the accumulation amounts fit the measured vertical radar stratigraphy. We apply different approaches to simulate the vertical accumulation distribution, which are described in the Results section.

3.4 Flowpath modeling for investigating the age structure in the ice column

Apart from investigating the recent accumulation conditions, the second aim of the study was the assessment of the Fedchenko accumulation basin as a potential climate archive. We reconstruct ice trajectories from the Jasgulem Pass into the accumulation basin with a simple model approach, in order to estimate the age/depth distribution and to assess the potential to retrieve an intact climate record from the glacier. For this purpose, we use the simplified approach of ice deformation along a flowline with parallel ice flow (Konrad and others, 2013). This model can be driven with the data collected at Fedchenko Glacier: glacier geometry and spatial accumulation distribution. The result will be a distributed age–depth relationship along the model domain, which is used for a comparison of potential drill sites.

A 2-D velocity field is calculated, based on ice thickness distribution and surface mass balance. The mean horizontal velocity is derived from mass conservation assuming no significant change in surface elevations, which is likely for the upper Fedchenko Glacier for at least the last century (Lambrecht and others, 2014). The depth-dependent horizontal velocity can then be derived from an analytical model (Lliboutry, 1981) and the vertical component results from mass continuity. Basal melt is neglected in the model, while divergence/convergence can be included by a simple parameterization if necessary. A kinematic correction is introduced to ensure a zero vertical velocity at the glacier bed. Ice trajectories and the accompanying age can be derived from the 2-D velocity field. The required vertically integrated mean horizontal velocity is calculated from the continuity equation by dividing the integrated accumulation along the upstream flowline by the local ice thickness, assuming close to balanced conditions. In addition, ice thickness and surface elevation information is required along the flowline. Ice thicknesses are interpolated from nearby profiles for the lower part of the flowline, while surface elevation is derived from the available GNSS measurements. Firn densification is not considered in this model as a process, but the firn layer is reduced to a respective ice thickness by applying a mean density of 650 kg m⁻³. This results in a correction to the total ice thickness of the order of 8.5 m. The influence of this simplification on the trajectory through ice thicknesses of several hundred meters is negligible. Trajectories are then calculated stepwise along the flow path. For details about the mathematical formulation and the discretization, as well as the applicability of the model, we refer to Konrad and others (2013).

4. Results

4.1 Recent accumulation conditions at the field sites

Information from the snow pits and the shallow cores was used to determine the accumulation rates for up to 3 years prior to the field measurements. Even though ice layers of 1–2 cm thickness were detected at all sites, the profiles of the stable oxygen isotope ratios show a clear seasonality (Fig. 2).

Fig. 2. Stable water isotope profiles in the form of δ ¹⁸O measurements for the three sites are shown in (a) (Deuterium profiles look identical and have been omitted). S1 and S2 show a clear seasonality, while the snow pit of S3 did not cover a full annual cycle. Black dots in (a) indicate summer seasons in the isotope profiles. Uncertainies between 1 and 4 permil are considerably smaller as the seasonal signal. (b) Meltwater conductivity profiles. Note that for S2, the peak at ~380 cm has been identified to correspond with a continuous GPR reflector.

This indicates that meltwater percolation is minimal and the conditions correspond to the so-called cold infiltration zone (Shumskii, 1964). The potential disturbance by meltwater percolation can be expected to be clearly restricted to less than one annual snow layer. Accordingly, the layer structure remains well preserved at these sites. Calculating ordinary linear regression of Deuterium against δ ¹⁸O yields slopes of 7.8, 7.9 and 8.0 for S1, S2 and S3, respectively. These values are similar to the values previously reported for Fedchenko snow pits (Aizen and others, 2009), and lie close to the global meteoric water-line. Hence, they do not show significant signs of isotopic change (e.g. by melting and refreezing). In contrast, seasonality is not recorded in the density profiles (Fig. 3), which show a slight increase with depth, but are dominated by density variations from local weak zones (e.g. depth hoar, recrystallization events) and numerous of ice layers.

Fig. 3. In situ density measurements (in dark blue) for the three snow pits marked in Figure 1. Red lines denote the corresponding water equivalent depth. The light blue bars represent ice layers observed in the snow and firn samples. Note, the thickness of the bars represents relative ice layer thickness and is not scaled to the depth axis. The uncertainty of the density measurements is ~5%.

The maximum period covered by direct sampling at the field sites is 3 years. We use a combination of the stable isotope ratios, conductivity and the location of the ice layers to determine the total annual layer thickness: The end of summer level is defined where the δ ¹⁸O values are high, thicker ice layers or aggregations of ice layers occur, and conductivity is high (e.g. a depth of 3.8 m at S2). This combination of parameters is based on the assumptions that temperatures are still high at the end of the ablation season, most meltwater during the ablation period refreezes close to the surface due to cold night temperatures and impurities accumulate in the uppermost part of the snow pack. These assumptions are supported by our numerous observations of expressed end-of-season ice layers with clear dust accumulation across the lower Fedchenko Glacier accumulation region since 2009. Layers with a high conductivity can also occur during the winter season due to, for example, dust deposition. An example is the conspicuous conductivity peak at 2.1 m in S2 and corresponds to a conductivity maximum at 1.2 m depth in S1. For both peaks, the δ ¹⁸O profiles show a local maximum within a period of low values, which corresponds to the winter season. The annual accumulation rates are calculated for the period September until August, which represents the annual layering during the field investigations and matches the end of the summer period with a higher probability of thick melt layers. For S1 (the longest record), this results in a mean accumulation rate of 1025 mm w.e. a⁻¹ (Table 2) between 2012 and 2015 and 1094 mm w.e. a⁻¹ for the period 2012 until 2016.

Table 2. Annual layer thickness, density and accumulation rate for S1

Please note that the mean density is scaled with the layer thickness and the snow pit data in 2016 (first line) are not included in the mean values.

Samples from S2 cover somewhat more than 1 year, according to our seasonality criteria. The layer thickness for the annual period 2014–2015 is 378 cm with a mean density of 536 kg m⁻³ and a mean accumulation rate of 2070 mm w.e. a⁻¹. The pit at S3 does not cover a full year. Therefore, the pit depth of 195 cm with a mean density of 516 kg m⁻³ does not allow to infer an annual accumulation rate. The overall mean δ ¹⁸O ratios are −16.06, −15.49 and −15.9 permil for S1, S2 and S3, respectively. Notably, these values appear to be significantly higher than average values ~−17.2 permil reported by Aizen and others (2009) for firn cores drilled at comparable altitudes in the same basin, albeit comprising the period of 2002–2005. The study of Aizen and others (2009) found a mean accumulation rate of 1380 mm w.e. a⁻¹ for 2002 until 2005 at a location close to S1 (at 5206 m elevation) and a mean accumulation rate of 2090 mm w.e. a⁻¹ close to S3 (at 5365 m elevation).

Data from the HAR dataset (Maussion and others, 2014) for the gridpoint corresponding to the central accumulation basin (model elevation 4964 m) indicate a mean accumulation rate of 1947 mm a⁻¹ for the entire period (2001–2014) and a minor negative trend (−20 mm a⁻¹) with rather high year-to-year variability. The mean accumulation rate results to 1020 mm w.e. a⁻¹ for the period 2012–2015, which is rather similar to our observations of 1094 mm w.e. a⁻¹ at S1 (Table 2). However, our investigations show increasing accumulation since 2013. The difference in observed accumulation rates between S1 and S2 is almost 800 mm w.e. a⁻¹ or 175%. This is in the same range as the difference of 152% found by Aizen and others (2009) between their two sample sites.

4.2 Areal accumulation distribution

We used the 600 MHz GPR data for this study, because the 200 MHz data provide no additional information, while the 600 MHz data have a better vertical resolution. The firn–ice transition was detected at 25–35 m depth in both data sets, with no visible internal layering below. The continuous reflection horizons identified in the high-frequency radar data can be attributed as seasonal markers according to the conductivity investigations in the snow pits and firn cores (Fig. 2) and by assuming the following seasonality criteria: accumulation variation is <100% (layer distance is less than double the mean layer distance) and layer reflectivity indicates a melt event with accumulated impurities (see also Section 4.1). It was possible to determine the near-surface radar wave velocity at S1 and S2 based on the identification of reflection horizons as well-established end-of-summer ice layers found at the field sites. An example is S1, where the summer layer was found at 2.1 m depth in the pit, which corresponds to a rise in conductivity in Figure 2b and the depth of the reference reflector at the end of the radar profile in Figure 4a (the location of S1). The same agreement was found for S2 (at 3.8 m depth, corresponding to the beginning of the radar profile in Fig. 4b).

Fig. 4. Ground-penetrating radar profiles Line A at the base camp (a, 5191–5215 m), Line C at Jasgulem Pass (b, 5295-5316 m) and Line Bin the eastern part (c and d, 5203-5273 m), for positions see Figure 1. For Line B also the terrain corrected geometry is given (d), in order to show the true geometry of the water table in the firn pack. The reference horizon for determining the near-surface wave velocity is shown as red line. The black arrows indicate the additional identified annual layers (dashed line in (a): indications of the reflector for 2011).

Several prominent reflectors were identified in the majority of the radar profiles by comparing their two-way travel times at consecutive profile sections or profile crossing points (Fig. 4). While the layer architecture appeared largely homogeneous, considerable differences in accumulation amounts were detected. Only minor indications of surface melt (e.g. low humidity of the surface snow) were observed during the field visits, especially in the higher regions of the basin. In contrast, some meltwater migrates and accumulates above ice layers in local depressions of the central basin (e.g. Line B at 5200 m a.s.l., Figs 4c, d). The reflections are generally stronger and continuous in the higher regions (Fig. 4b), indicating less surface melt and considerably less water migration. The better reflection conditions in these regions also allow the identification of more horizons. The layer thickness distribution is spatially consistent along all profiles, which indicates a largely stable accumulation pattern over time across the accumulation basin.

The annual layer thickness can be converted into accumulation amounts for the uppermost layers, because we can calibrate the radar profiles at the sampling sites with the measured density profiles. For the spatial interpolation across the basin, it was assumed that the density is elevation dependent (compaction in the surface layers is a function of near surface air temperature). The elevation gradient was then inferred from the elevation difference between S1 and S2.

A multivariable regression analysis indicates that the distance from the western basin margin clearly dominates the accumulation distribution compared to the surface elevation. The normalized β coefficient for the easting component is about eight times larger as for elevation. A similar weak elevation dependence of accumulation is also observed at Abramov Glacier, ~100 km to the North (Barandun and others, 2015). The distributed accumulation field (Fig. 5) confirms the strong gradient of accumulation between the western rim of the accumulation basin and its center. This situation is also found in the results of Aizen and others (2009). The spatial pattern suggests that the major precipitation events arrive from the West, depositing large amounts of snow just at and behind the top of the steep precipice west of the basin which rises almost 1500 m from the Abdu Kagor Valley. Potential wind-driven snow redistribution probably does not reach far into the basin, considering the strong gradient close to the western basin boundary. Some 3 km behind the crest and ~150 m lower, the accumulation rates are already reduced by ~900 mm w.e. (~44% of the maximum values).

Fig. 5. Distributed accumulation conditions for the uppermost annual layer, representing the accumulation period August 2014 to August 2015. This layer was identified as the first reflector in the radar data and calibrated with the respective summer horizon at the sampling sites S1 and S2. Interpolation is done by the Topotoraster tool of ArcGIS based on the radar profiles shown in the figure. The background image is based on a Landsat 8 scene from 2013.

4.3 Temporal accumulation evolution

The temporal evolution of the accumulation rates can be derived from dated layers in the radar profiles, given that the density of the individual layers is known. The high-frequency radar profiles show a characteristic vertical structure with a clear layering in the uppermost 200–300 ns two-way travel time, which corresponds to ~19–28 m with the mean bulk wave velocity in firn. The lower parts of the radar data, in contrast, show almost no reflections. This vertical structure indicates that the firn to ice transition occurs in this depth range and accumulation rates can be retrieved for at least the last 7 years before the time of our field measurements (2008–2015). It is possible to track most of the prominent layers across the entire basin by comparing cross points and adjoining profile sections. Unfortunately, it was not possible to identify a precise layer depth at S1 for 2011 in the radar data, because of an unclear reflection pattern, although the existence of this layer is obvious in the data (Fig. 4).

The firn stratigraphy is estimated from the identified radar reflectors, the air temperature and a mean accumulation rate using the extended firn densification model of Reeh (2008), which includes the possibility of incorporating refreezing of surface melt. Based on the stratigraphy, we derived distributed accumulation rates similar to Barandun and others (2015). The comparison of radar reflectors with the snow pit and firn core stratigraphy allows the determination of a characteristic annual layer thickness, because the end of summer layers are connected to thick ice lenses with a high dust concentration for all our field observations since 2009, while no strong change in stratigraphy occurs during winter. This leads to the identification of annual radar reflectors and thus the stratigraphy in the radar domain.

Air temperatures at a modeled pixel elevation of 4964 m from the HAR dataset are compared to air temperatures measured at an automatic weather station at 4915 m between October 2009 and September 2012 (data provided by K. Fujita, pers. Comm. Nov. 2018). Monthly mean temperatures of the two data sets show a correlation of R ² = 0.77, which supports the use of HAR air temperatures in the firn densification model. Also the degree-day sums of 445 °C for the weather station and 424 °C for the HAR data match very well after applying a mean lapse rate of 0.007 °C m⁻¹ for the overlapping period between station data and HAR data.

An estimate of the multi-year mean precipitation rate is required as an initial guess for the densification model. We used the mean precipitation rates from the corresponding HAR pixel for September to August periods scaled by our direct measurements as initial mean accumulation rate. The HAR precipitation rate for September 2008 until August 2014 is 1931 mm a⁻¹, while the mean HAR precipitation rate for the period overlapping with our field data (September 2012 to August 2014: 947 mm) is 1717 mm a⁻¹. Therefore, the scaled initial mean accumulation rate is 1065 mm a⁻¹ for S1 and 1849 mm a⁻¹ for S2. The initial surface layer density is taken from the field measurements. The layer thicknesses derived from the radar data are the basis for a first density/depth distribution calculated by the model. Assuming a bulk water content of 3.4% for S1 and 2% for S2 results in a perfect fit of the radar reflectors for the layer depths at these field sites, when applying the snow pack mixing model of Frolov and Macheret (1999). For this analysis, we used the mean snow densities in the snow pits and firn cores and calibrated the wave velocities by the measured reflector depths and conductivity profiles. The resulting bulk velocities are 0.185 m ns⁻ at S1 and 0.191 m ns⁻ at S2.

The vertical density profile derived by the above strategy is then used to improve the layer thicknesses in an iterative procedure until an optimal match between layer thicknesses and measured accumulation rates of the upper layers is reached. The resulting accumulation rates are displayed in Figures 6 and 7. We also used different velocity models (Kovacs and others, 1995; Looyenga, 1965; Frolov and Macheret, 1999) which cover the spread of potential radar wave velocities in the firn pack to investigate the robustness of our layer depth inversion. The resulting variability of accumulation rates is a measure of the error range of our reconstruction of accumulation history.

Fig. 6. Results of the accumulation reconstruction for S1 from radar layers (purple triangles) and the firn densification model (blue diamonds). The results are also compared to the HAR precipitation scaled to the accumulation at S1 (dark green x) and the field results (red diamonds). The error bars represent the spread of the accumulation rates, based on different reasonable radar velocity models. Note that the layer approach is not possible for the year 2011 (no suitable, clearly identifiable radar reflector). Therefore, the sums for the mass-balance periods 2010 and 2011 are displayed for the year 2010 (blue solid diamond). The sum of the HAR precipitation for 2010 and 2011 is displayed as yellow cross for comparison.

Fig. 7. Results of the accumulation reconstructions for S2, with the same approach as in Figure 6.

The scaled HAR reanalysis precipitation data follow the accumulation rates derived from the radar data quite well (high and low annual amounts), even though deviations can be several 100 mm. The interannual variability of accumulation is up to 1000 mm. The resulting mean accumulation rates, based on the ensemble means of the reconstructions, are 1223 mm for S1 and 2058 mm for S2, respectively, between August 2008 and August 2015.

The depth of the deepest identifiable annual layer was at ~22 m. Applying the scaled HAR precipitation to the densification model results in a firn density of 850 km m⁻³ at 30 m depth, which corresponds to our observations of losing radar reflectors in this depth range. The densification model estimates the amount of refrozen meltwater at S1 to a maximum of 16% of the accumulation rate in 2015, while the mean value is ~6%. At S2, the maximum amount of refrozen water is estimated to be 2%, while the mean value is <1%, considerably less than at the lower site of S1. This is also in accordance to the radar observations, where the reflectors in the upper regions of the accumulation basin are clearer even for deeper parts of the firn pack. In the lower regions, most energy is already reflected from the upper and more massif ice lenses, blurring the stratigraphy below.

4.4 Ice thickness in the accumulation basin

The knowledge of ice thickness is required for investigations on ice flux and ice volumes, and is also an important factor with respect to ice core site selection. We use this information also to evaluate the age–depth structure of the uppermost Fedchenko Glacier. Bed reflection was well detectable along all profiles after processing and could be determined semi-automatically (Fig. 8). The results of the automatic picking tool were manually checked and corrected at obvious deviations of the automatic detection from the bed reflection. The ice thickness reaches maximum values of more than 700 m at the central lower gateway of the accumulation basin (Fig. 9). At the center of Jasgulem Pass, the highest region of the main glacier, the ice thickness is still more than 400 m. In contrast, the small tributary basin toward the northwest shows ice thicknesses of ~200 m only. The ice thickness data prove that a well-expressed trough exists along the main glacier, extending up to Jasgulem Pass.

Fig. 8. Example of a low-frequency GPR profile for ice thickness determination in the lower part of the accumulation basin. The strong reflection between 5000 and 7500 ns two-way travel time (left axis) represents the bedrock, while some internal reflections can be identified above 3000 ns.

Fig. 9. Color-coded ice thickness along the ground-penetrating radar profiles (left) and the corresponding total error (right). The flowline for the trajectory modeling, starting at the ice divide, is indicated as a thin blue line in the left panel. The background image is based on a Landsat 8 scene from 2013.

4.5 Age–depth structure along a flowline

The calculation of the age with depth requires the knowledge of the ice deformation along the advection path. Because the respective boundary conditions and the driving data are not known in great detail, we use the simplified approach described in Konrad and others (2013), which requires glacier geometry and accumulation distribution as input, both measured at Fedchenko Glacier. The flowpath (Fig. 10) is modeled in two dimensions assuming parallel flow along the flow direction, which is a reasonable assumption for most of the profile based on the local topography. The introduction of convergence for the lower part of the flowpath (not shown here) has a negligible influence on the age distribution compared to uncertainties in accumulation rate, which was also found by Konrad and others (2013).

Fig. 10. Glacier geometry and modeled trajectories along the central flowline from Jasgulem Pass (at x = 0 m) to the central accumulation basin. The vertical bar represents a depth profile at 700 m distance from the divide (see Fig. 11).

Fig. 11. Age/depth relation at 700 m distance from the ice divide at Jasgulem pass (see Fig. 10 for the sample trajectories). The relationship is shown for the measured accumulation distribution in the basin (blue line) and different mean accumulation rates. The results are displayed as a solid line for the upper 90% of the ice column and as a dotted line down to 98% (465 m) of the ice column. The ice thickness at this location is 474 m and bedrock is indicated by the brown line at the bottom.

A major limitation exists due to the initial assumption of time-independent ice thickness and accumulation rate over the simulation period. However, earlier investigations showed that elevation change was negligible since the beginning of the new millennium and thinning was of the order of 6% of the ice thickness since 1928 (Lambrecht and others, 2014, 2018). The simulation aims at providing a reasonable estimate of a realistic age/depth relation and thus such variability might be tolerated. Calculating detailed particle trajectories (e.g. for investigating upstream catchment areas) is not a focus here, due to the considerable uncertainties in trajectory computation with the simple model approach (Bohleber, 2011).

In order to estimate the potential influence of accumulation variations (which are not known yet), we perform a simple sensitivity analysis with different mean accumulation rates of 2070, 1800 and 1300 mm w.e. a⁻¹, respectively. A mean accumulation rate of 2070 mm w.e. a⁻¹ represents spatially constant conditions found today at Jasgulem Pass, while 1800 mm w.e. a⁻¹ corresponds to the mean conditions along the flowline. Accumulation rates of 1300 mm w.e. a⁻¹ are a lower-end estimate of accumulation rates in the past. The model is not suitable to estimate the age directly at the ice divide due to differences in ice physics (no longitudinal stresses are implemented in the model). Therefore, we investigate the potential age/depth relation along-flow downstream of the ice divide and focus on a distance of 700 m (i.e. ~1.5 ice thicknesses) from the divide (Fig. 11). At this distance, the estimated age at a depth of 90% of the ice thickness (430 m) reaches ~1000 years for the modern accumulation distribution. We discuss the age variability based on the estimated age at 90% ice depth, as the uncertainties below this level make the age estimates less reliable towards the bed. For lower (but spatially constant) accumulation rates, which might be more likely for the past, the age increases to 1560 years for 1300 mm w.e. a⁻¹. An increase of the accumulation rate to 2400 mm w.e. a⁻¹ decreases the age to 860 years. At the lower end of the accumulation basin (4500 m from the divide), the expected age at 90% ice depth increases to 1200 years for modern accumulation conditions. However, this gain of 200 years in temporal coverage might be accompanied by a stronger distortion of the layer structure in the lower parts of the ice column and an increasing influence of lateral convergence. Also, the probability of a disturbance of the chronology by meltwater events increases for lower altitudes along the glacier flow. The layer thickness for our standard experiment (700 m downstream of the ice divide) decreases almost linearly from 2 to 3 m at the surface to ~12 cm at 380 m depth. Below this depth, the layer thicknesses further decrease to ~5 cm at 430 m depth.

5. Discussion

Our investigations of the snow and firn layers in the upper Fedchenko accumulation basin in combination with ice thickness measurements resulted in new insights into accumulation distribution and the potential for the retrieval of climate information from an ice core. The snow pit information enabled us to quantify the magnitude of the current accumulation rates at specific sites, while the shallow GPR surveys provided the information about their spatial variability. We were able to track several prominent layers along most of the GPR profiles, which show a stable accumulation pattern during the recent past. This pattern is characterized by a rather strong decrease in accumulation with growing distance from the western margin of the basin. Accumulation in the center of the basin is 0.57 of the amount at its western limit. The spatial distribution of accumulation during 2014–2015 (Fig. 5) shows limited variability of 1200–1400 mm w.e. a⁻¹ in the central basin compared to the 450–800 mm w.e. a⁻¹ in the sheltered regions closer to the high mountain slopes of the main range. The year-to-year changes in accumulation during 2008–2015 vary by ~800 mm, from 890 to 1700 mm. However, the standard deviation from the mean value of 1223 mm is only 260 mm in the basin and 370 mm for the mean accumulation rate of 2058 mm at the margin of the basin (S2).

The derived accumulation rates and the reanalysis data show a similar year-to-year variation. However, the spatial differences to the mean HAR value are considerable. Figure 5 represents almost the size of one HAR pixel with a mean accumulation rate of 1936 mm (extrapolated to the investigated period 2008–2015). In contrast, the measured values vary by 835 mm between S1 and S2, with the HAR precipitation being much closer to the value of 2058 mm determined at S2. This demonstrates that the modeled absolute precipitation values may differ substantially from the local mean values, even though the temporal pattern is rather consistent with the locally derived evolution. Compared to the results of Aizen and others (2009) covering the period 2002–2005, the accumulation rates seem to be rather stable for more than a decade (e.g. 2090 mm a⁻¹ compared to 2058 mm a⁻¹ for 2008–2015 at S2).

The exposed location of the Fedchenko Glacier accumulation basin at the top end of long valley systems oriented southwest results in rather high accumulation rates, compared to other locations, because the valleys enable the transport of moist air masses toward the glacier basin. Ice core investigations at Muztag Ata in the eastern Pamir revealed mean accumulation rates of 605 mm w.e. a⁻¹ at 7010 m a.s.l. between 1958 and 2002 (Duan and others, 2015). A reconstruction of mass balances for Abramov Glacier, North of Fedchenko Glacier, results in accumulation rates of ~740 mm w.e. a⁻¹ at 4700 m a.s.l. for the period 2005–2014 (Barandun and others, 2015). Further to the East in the Central Karakoram, mean accumulation rates reached ~1300 mm w.e. a⁻¹ at Hispar Dome at 5450 a.s.l. in the mid-1980s (Wake, 1989), while values at the north-eastern rim of the Karakoram were ~1060 mm w.e. a⁻¹ at a similar elevation on Urdok Glacier in 2006 (Mayer and others, 2014). While the rather low values at Muztag Ata might be due to the high altitude of the drill location, the sites in the Karakoram are comparable in height and surface topography. This comparison underlines the fact that the accumulation rates of the central basin are comparable with other sites, but the near-margin values are exceptional, even though the investigation periods are different.

Surprisingly, summer melt events reach the high accumulation basin of Fedchenko Glacier, resulting in pronounced reflection horizons dating back at least to the year 2001. But it seems that there is a rather strong gradient within the basin. While internal water accumulation is observed in the lower basin (Fig. 4d), the upper regions close to Jasgulem pass show less evidence of meltwater and clearer radar reflectors. This is consistent with our firn model experiments which indicate rather small amounts of refrozen meltwater at the upper sample locations. Therefore, preservation of accumulation layers is more likely near the upper margin of the basin, whereas the lower basin is not as suitable for ice coring due to internal water accumulation.

The flow-model based age estimate for the site 700 m downstream of the ice divide should be regarded as a minimum estimate as 2-D models have the tendency to underestimate thinning and age of lowermost layers (Konrad and others, 2013). Accordingly, we regard the model results merely as lower age constraints. A comparison with age estimates for the divide, using the model of Kaspari and others (2008) results in a similar age to depth distribution. Also, the functions in Figure 11 indicate that considerably older ice can be found in the lowermost part of the ice column, especially as ice flow is very slow and strong perturbations of the stratigraphy are not very likely even in this lower part. Together, the results suggest that a core site near the divide would contain information of more than a millennium of climate history with annual layers >5 cm in thickness, providing excellent temporal resolution and a detailed climate history. In addition, it is likely to retrieve valuable climate information for the earlier part of the Holocene, however with thinner layers.

Our investigation further illustrates the benefit of combining GPR with snow pits, shallow firn cores and simple ice flow modeling for in situ ice core reconnaissance. The main difference between this study and the study of Colle Gnifetti in the European Alps by Konrad and others (2013) is in the scale of the two glaciers and the density of available data. In this sense, our results are more closely related to earlier studies at km-scale accumulation basins in the Alps (Sold and others, 2015) and Svalbard (Pälli and others, 2002). A particular parallel exists to the exploration of ice core drilling sites in the Alaska Range, where a similar combination of GPR, glaciological and glaciochemical reconnaissance in concert with simple flow modeling was used to predict a millennial-scale ice core record at Mount Hunter (Campbell and others, 2012). Annual layer counting in subsequently drilled deep ice cores confirmed this general age span (Osterberg and others, 2017), which underlines the significance of such investigations for identifying suitable locations for the retrieval of useful climate archives.

6. Conclusion

The investigation of both spatial and temporal accumulation distribution for the central accumulation basin of Fedchenko Glacier show strong variations. The annual accumulation rates vary by a factor of almost two (1100 mm w.e. for S2) during the period 2009–2015, even though the Std dev. is ~20% (370 mm w.e.). The large maximum difference is due to the succession of a year with high precipitation in 2010 and the very dry year 2011. The general temporal pattern fits the reanalysis data, but magnitudes are considerably different. One reason is the existence of strong spatial gradients across the accumulation basin. Annual accumulation varies from 870 mm w.e. in the eastern sheltered region to 2400 mm w.e. at the western basin boundary in the balance year 2014/2015. This pattern is consistent with time and is mainly determined by the effect of higher snowfall close to the steep topography at the western margin and is not resolved by the 10 km resolution reanalysis data. The HAR data seem to show too high values, which are close to the maximum accumulation rates observed in the basin. However, it is difficult to compare both datasets, because several parts of the accumulation basin could not be sampled due to safety issues (avalanche risk) or accessibility. Still, our investigations indicate that the HAR accumulation values are considerably higher than the area-mean observations. As a consequence, long-term mass-balance calculations for the entire glacier which are based on reanalysis data require an accumulation correction based on in situ observations. Even though the year-to-year variations are high, the mean accumulation rates are rather constant since 2002, indicating stable accumulation conditions in the new millennium.

The combined study of glacier geometry and accumulation conditions suggests that this basin represents a major climate archive in Central Asia. The altitude above 5200 m is just high enough to guarantee the preservation of annual layers in the ice column, while the lowermost parts of the basin are already affected by surface melt. Snow pits and firn cores close to the ridges and passes also show summer melt layers, but no meltwater penetration to deeper levels. Based on the local conditions of accumulation and ice thickness distribution, it seems that a particularly suitable location for the retrieval of a climate record is at some distance from the divide. The maximum ice thicknesses are found in the glacier branch reaching up to Jasgulem Pass. Therefore, a promising ice core drilling location might be at a few ice thicknesses downstream of Jasgulem Pass, where the modern accumulation rate is on the order of 2000 mm w.e. a⁻¹ and the ice thickness reaches ~470 m. A considerably higher temporal resolution compared to potential drill sites in the Eastern Pamir can be expected due to the more than two times higher accumulation rates (Duan and others, 2015). At least ice from the last millennium could be recovered from the central accumulation basin of Fedchenko Glacier with annual layer thicknesses of the order of 5–10 cm. Much older ice can be expected from the lower 10% of the ice column reaching possibly back to the early Holocene. However, potential disturbances in the layer stratigraphy and the influence of basal conditions cannot be ruled out by our simple analysis. In view of the high-resolution and the long potential time span of the record, this sampling site could provide a highly valuable long-term paleoclimate record from Central Asia.