Introduction
Glaciers are potentially good climatic indicators. However, most of the source information on the dynamic response of temperate glaciers to changing climate comes from European and, to a lesser extent, North American glaciers, although European glaciers form only 3% of the world glacier area (Reference Dyurgerov and MeierDyurgerov and Meier, 1997). Hence, our knowledge of the glaciers of remote and inaccessible mountain regions, often characterized by a high degree of glacierization, is rather limited. This biased distribution of glaciological source information severely constrains present attempts to assess the impact of climate change on the global glacier volume. It also prevents a full understanding of glacier dynamics and the impact of different climatic conditions. This paper seeks to expand the current global glaciological record by presenting and analyzing data gathered on Sofiyskiy glacier, located in the Russian Altai mountains (Fig. 1), during four consecutive summer field seasons in the period 1997–2000. The climate of the Altai mountains has two important characteristics which make this study area especially interesting. First, the extreme continentality gives rise to aridity and a large diurnal and seasonal temperature range. Second, owing to the coincidence of accumulation and ablation seasons in summer, we can characterize the Altai glaciers as summer-accumulation-type glaciers. Although current knowledge of the behaviour of these glaciers is limited, they are reported to be more vulnerable to the current global warming trend than the better-known (maritime) winter-accumulation-type glaciers (Reference Fujita and AgetaFujita and Ageta, 2000). In this paper, the basic characteristics of Sofiyskiyglacier are discussed, such as surface mass balance and velocity, and the historical front variations. To improve understanding of the glacier’s dynamics, we also present a force-balance analysis and an analysis of basal conditions using a three-layer reflectivity model for electromagnetic pulses.
Fig. 1. Location map of the Altai mountains and the study area (Sofiyskiy glacier).
Study Area and Recent Glacial History
Sofiyskiy glacier (49°47’N, 87°46’E) is one of the major glaciers in the central part of the south Chuya range, a mountain chain within the Altai mountain belt, situated in the border region between Russia and Mongolia (Fig. 1). At present, the glacierized area consists of three basins (Fig. 2), i.e. the main basin (A) and two smaller basins (B and C). The latter are at present no longer connected to the main ice flow. The southwest–northeast-oriented Sofiyskiy glacier is approximately 7 km long and has an average width of 1.7 km. The head of the glacier is located near Brat mountain (3876ma.s.l.), one of the highest peaks in the range. The cumulation area has a mean surface slope of ~ 17% at the central flowline. A steep icefall is situated between 2800 and 3000 m a.s.l. Below this icefall, the glacier terminates in a relatively flat tongue.
Fig. 2. Contour map of Sofiyskiy glacier, based on the 1952 topographic map. At present, basins B and C are no longer connected with the main ice flow.
During the 20th century, several front-position markers were recorded. The earliest records are from field markers and aerial photographs (Reference Bulatov, Dik and RevyakinBulatov and others, 1967). During each of the 1997–2000 field seasons, the snout position was mapped with differential global positioning system (DGPS). The historical record is extended further back in time by radiocarbon dating of wood fragments collected in three different former front moraines of Sofiyskiy glacier (Table 1). Comparison with Maliy Aktru glacier, located in the north Chuya range some 30 km to the northwest (Fig.1), shows that both glaciers have retreated steadily since the beginning of the 20thcentury. Average retreat rates are 9.6 ma–1 for Maliy Aktru glacier and 18.3ma–1 for Sofiyskiy glacier, based on remeasurement of field markers and our own measurements of the front position. The most rapid retreat of both glaciers occurred between 1900 and 1940 (Fig. 3).
Fig. 3. Historic front positions of Sofiyskiy glacier (°) and Maliy Aktru glacier (*) during the 20th century.
Table 1. Historic front positions of Sofiyskiy glacier, measured as distance from present front position
Surface Mass Balance and Glacier Dynamics
A network of about 35 stakes was employed for determining local mass balance and ice velocity. Owing to the presence of many hidden crevasses and severe snow conditions in the accumulation area, this network was confined to the zone below 3100 m a.s.l. The position of the stakes was measured with a Leica SR 9500 DGPS. This system consists of a fixed reference station (“Camp” in Fig. 2) and a remote “rover” station. The system assures a centimeter accuracy in horizontalposition and an accuracy of 3–5 cm in elevation. Base lengths were always <10km and usually <5km, further limiting the error. However, taking into account errors due to the repeat survey, and positioning errors of the rover station due to surface melt, the overall accuracy should be somewhat lower.
Glacier mass-balance measurements are conventionally carried out at the end of the ablation season. Due to logistic restrictions, these mass-balance measurements had to be carried out some weeks beforehand. Accumulation was determined from a number of snow pits dug in 1999. Figure 4 shows the individual mass-balance measurements for each balance year. It should be noted that the ablation values for 1999–2000 are probably underestimated, as the ablation stakes were remeasured early in the summer season. Based on the individual measurements, a mean mass-balance reference profile for the period 1997–2000 was constructed using a linear parameterization (Fig. 4). Excluding the short 1999/2000 balance season, the equilibrium-line altitude (ELA) and the slope of this reference profile are 3179 m and 73.1 cm a–1w.e. per 100 m, respectively. Figure 4 also shows the mass-balance reference profile for Maliy Aktru glacier, based on a 38 year mass-balance record (1962–99; Reference Hoelzle and HaeberliHoelzle and Haeberli, 1999). The mean ablation at the glacier tongue is distinctly lower for Maliy Aktru than for Sofiyskiy glacier. This might be due to the fact that Maliy Aktru glacier lies in a very narrow and steep valley, and is therefore more shaded from direct radiation.
Fig. 4. Mass-balance measurements carried out on Sofiyskiy glacier for the mass-balance years 1997/98 (*), 1998/99 (°) and 1999/2000 (+). The solid line is the best fit of these data. The dashed line is the mean mass-balance profile for Maliy Aktru glacier based on a 38 year mass-balance record (Reference Hoelzle and HaeberliHoelzle and Ha eberli, 1999).
Based on these stake measurements, a profile of the average horizontal surface velocities along the central flowline was compiled (Fig. 5). This profile shows a maximum surface velocity of 110 m a–1 at the icefall; the average surface velocity is 35 m a–1 over the glacier tongue and 55 m a–1 over the entire ofile.
Fig. 5. (a) Surface and bedrock profile (gray shaded), surface velocity measurements (*) and modeled surface velocities using A = 0.8x10–16 (white circles) and using A = 1. 6x10–16Pa–3a–1 (black circles) along the centralflowline of Sofiyskiy glacier. (b) Driving stress τd (solid line), predicted basal shear stress τxz (dashed line) and predicted surface longitudinal stress deviator τ′xx calculated with the model along the central flowline.
Force-Balance Analysis
To gain a better insight into Sofiyskiy glacier dynamics we applied the higher-order numerical flowline model of Reference PattynPattyn (2002) in a diagnostic way. Besides shear stresses, this model takes into account longitudinal stress gradients in calculating the ice deformation, and uses Glen’s flow law as a constitutive equation. Based on the local ice geometry along the central flowline (surface and bedrock topography, width of the drainage basin), a constant flow parameter (experimentally determined as A = 0.8 × 10–16 Pa–3 a–1; Reference De Smedt and PattynDe Smedt and Pattyn, 2003), and considering a no-slip basal boundary condition, the full stress field was determined. Major components of the stress field are displayed in Figure 5. The icefall area is characterized by high driving stresses (τd > 300 kPa), and for most of the long profile τd > 100 kPa. The basal shear stress is a smoothed version of the driving stress, characterized not only by a lower amplitude, but also by phase shifts compared to the driving-stress profile. Such phase shifts are due to local pushes and pulls as a result of bedrock and surface perturbations (Pattyn, 2002). Longitudinal stretching is confined to the region upstream of the icefall. The downstream area (ice tongue) is dominated by compressive flow.
Modeled surface velocities are in agreement with the measured values upstream of the icefall. For the whole downstream section, measured surface velocities are significantly higher, especially in the icefall region. The discrepancy is probably due to neglect of basal sliding in the model (no-slip boundary condition). We carried out an experiment with a higher value of A = 1.6 × 10–16 Pa–3 a–1, hence making the ice softer (warmer). While softer ice accounts for a better agreement between modeled and observed surface velocities in the ice tongue, the discrepancy becomes larger upstream of the icefall (modeled velocity too high) and at the icefall itself (modeled velocity still too low). These results thus indicate that (i) deformational velocity alone cannot account for the spatial distribution of the velocity field of the glacier and that a significant amount of the total surface velocity must be due to basal sliding in the downstream area of the glacier, or (ii) the glacier ice is softer (hence warmer) in the ablation area compared to the area upstream of the icefall.
Radio-Echo Sounding
Radio-echo sounding (RES) measurements were carried out with a 5 MHz (central frequency) ice-penetrating radar (Reference Narod and ClarkeNarod and Clarke, 1994). The system consists of a monopulse transmitter generating 1600 V pulses across a resistively loaded 10 m dipole antenna, and an airwave-triggered oscilloscope receiver connected to a palmtop computer for digital data recording. Pulses are generated at 512 Hz and have a bandwidth of 1–200MHz. Similar systems have already been used successfully on several glaciers in Canada and Alaska (Reference Narod and ClarkeNarod and Clarke, 1994; Reference Welch, Pfeffer, Harper and HumphreyWelch and others, 1998; Reference Copland and SharpCopland and Sharp, 2001). Apart from a number of isolated RES measurements, four profiles were taken in crevasse-free zones, which reveal details of the bed topography and basal conditions. At ∼ 3450m (“X” in Fig. 2), a longitudinal (along the central flowline) and a transverse radar profile were recorded. The horizontal distance between successive measurements was -19 m in both profiles. In the ablation zone beneath the icefall, a second longitudinal-transverse profile set was recorded with a horizontal interval of ∼9.5 m (the central point is “Y” in Fig. 2).
On all measurements, static corrections were performed for the trigger delay caused by the 38 m separation between transmitter and receiver. Band-pass filtering was not carried out. Migration was calculated along the one-dimensional profiles, but is not considered in the analysis below due to the complex topography in some places. As outlined in Reference Copland and SharpCopland and Sharp (2001), migration would only serve to increase values of high bed reflection power (BRP), as it tends to increase ice thickness in zones of steep basal topography. All recorded traces present relatively clear basal reflections. An example of this is given in Figure 6, which shows a line-intensity plot for the longitudinal glacier profile in the ablation zone. The two-way travel time and amplitude of the bed reflection wave were determined using a semi-automatic bed-reflection picking routine. Assuming that the bed reflection occurs from a point directly below the midpoint between the transmitter and the receiver, ice thickness was inferred using a wave propagation velocity of 168 μm s–1 for glacier ice, a value confirmed by a common-midpoint experiment carried out in the ablation area. A density profile in the accumulation area was obtained from Reference FujiiFujii and others (2002), and a propagation velocity correction was determined by considering firn as a mixture of ice and air. According to Reference LooyengaLooyenga (1965), the dielectric permittivity for such a heterogeneous mixture is given by (e.g. Reference Hempel, Thyssen, Gundestrup, Clausen and MillerHempel and others, 2000)
Fig. 6. Line-intensity visualization of the longitudinal RES profile in the ablation area of Sofiyskiy glacier.
where εice = 3:17 is the dielectric permittivity of ice. Combining the local ice thickness with the DGPS registration of the surface position of the traces yields four basal profiles (Fig. 7). Owing to the large number of multiple basal reflections with in the reflected RES signals (most likely caused by a complex three-dimensional topography) and the relatively large distance between two successive measurements (-19 m), the accuracy of the bed position in the accumulation area varies between 15 and 45 m, a value estimated from the uncertainty with which the bed position could be determined due to multiple reflectors. I n areas with clear basal reflections, the uncertainty is typically of the order of one-tenth of the transmitted wavelength(Reference Bogorodsky, Bentley and GudmandsenBogorodsky and others, 1985), leading to an error of ±3.4 m. Our measurements suggest a maximum ice thickness of > 200 m in the accumulation zone (Fig. 7a). Bedrock signals are distinctly clearer in the ablation zone. Here, the 900 m longitudinal profile starts at the bottom of the icefall, where basal reflections are disturbed (Fig. 6 and 7d). This is most probably due to the steepness and roughness of the local glacier bed and internal structures within the ice. Beyond the icefall, the bed steepness suddenly reduces to 4°, and basal reflections become more defined.
Fig. 7. Bedrock and surface along the two transverse RES profiles on Sofiyskiy glacier in the accumulation area (a) and the ablation area (b), and along two longitudinal RES profiles in the accumulation area (c) and the ablation area (d). The cross-profiles run from south to north; the long profiles run from west to east. X and Y denote the position of the longitudinal profiles and are shown in Figure 2. Also given are BRPR (thin solid line) and IRPN (dotted line).
To gain insight into the basal layer, it is convenient to calculate the reflected power. We consider the reflected signal to be consistent with a monochromatic sinusoidal wave train (Reference Gades, Raymond, Conway and JacobelGades and others, 2000). As such, we can define the reflected power P within a time window t1 – t2 as
where Ai is the recorded amplitude of the wave. Using Equation (2), we can calculate the reflection power of both the bed reflection wave (BRP) and the internal reflection waves (IRP), which are located between the end of the direct airwave and the beginning of the bed reflection wave (Reference Gades, Raymond, Conway and JacobelGades and others, 2000). The BRP contains information on the structure and composition of the glacier bed. The IRP is a measure of the scatter and attenuation within the glacier ice. Contrary to Reference Gades, Raymond, Conway and JacobelGades and others (2000), a variable time window is used for the calculation of the IRP, which is more suitable for shallow glaciers (Reference Copland and SharpCopland and Sharp, 2001). This window extends from 0.6 ms after the beginning of the airwave to 0.1 ms before the bed reflection wave. The BRP starts at the latter point and ends at 0.4 ms after the bed reflection wave.
If we assume that losses due to polarization and scattering at the ice/bed interface are minimal or constant, and antenna characteristics and surface coupling constant for all measurements, the only controls on the BRP are ice thickness, the loss due to two-way propagation through the ice (attenuation) and the strength of the basal reflection (Reference Copland and SharpCopland and Sharp, 2001). To determine the actual BRP variations due to variations in bed properties, we thus need to account for the effect of varying ice thickness. This is done by establishing an empirical relation between the bed reflection and the two-way travel time, using all traces of the four profiles (Fig. 8). To reduce the scatter due to differences in reflection properties of different materials, data over a homogeneous section in both the accumulation and ablation areas were selected. Data from the accumulation area were subsequently corrected for density variations in firn. Best-fit curves for both areas are also drawn in Figure 8. The best fit on all data was found to be an exponential function with r2 = 0.59. Even though the accumulation data are corrected for density variations, the bulk of these data lie below the best-fit function, confirming a relatively low BRP for these measurements. Following Reference Gades, Raymond, Conway and JacobelGades and others (2000), we now define the depth-corrected reduced BRP (BRPR) as the measured BRP divided by the calculated BRP (according to the empirical relation displayed in Figure 8). Similarly, the normalized internal reflection power (IRPN) is defined as the measured IRP dividedby the mean IRP of all profiles. Although the BRP shows a high scatter, and the correlation coefficient of the best fit in Figure 8 is relatively low, the scaling operation has removed the major trend, i.e. the bias of high BRP associated with shallow ice thickness. The scatter will not hamper the interpretation of BRPR, as we will only consider the relative importance of BRPR in all four profiles.
Fig. 8. Empirical relation between the bed reflection and the two-way travel time, using traces of the four profiles on Sofiyskiy glacier from the accumulation (filled circles) and the ablation area (open circles). Best-fit curves for both areas are given dotted lines. The solid line displays the best-fit curve for all data. Values of the accumulation area were corrected for density variations in the firn to determine the ice thickness.
In the accumulation area, values for BRPR and IRPN are distinctly lower than in the ablation area, and one can argue whether these differences are related to differences in englacial conditions (attenuation) or to differences in basal conditions (reflection). RES measurements on John Evans Glacier, Canada, also demonstrate this marked difference in BRPR between accumulation and ablation area (Reference Copland and SharpCopland and Sharp, 2001). If englacial conditions are responsible for this difference, warm ice should correspond to low BRPR. This contradicts the force-balance calculations which indicate that the accumulation area, where the lowest BRPR is recorded, consists of stiff (cold) ice. It is therefore likely that BRPR patterns are related to changes in basal conditions.
IRPN variations seem to be primarily caused by differences in the snow and firn density, the ice temperature and the presence of occasional water inclusions and ice lenses. Since internal structures affect the emitted RES signal before it reaches the glacier bed, there is often an inverse relation between IRPN and BRPR. As such a relationship is not clear from our records, and IRPN values are generally low for all profiles, the BRPR appears little affected by internal attenuation processes. In the ablation area, two zones with a remarkably high basal reflection emerge from the RES profile. The first is located between 150 and 350 m (Fig. 7d) and coincides with the sudden decline of the bed slope and the confluence area with side basins A and C. A number of hyperbolic structures can be identified which make it difficult to accurately determine the ice thickness. This could indicate a complex basal topography in this area and might influence the BRPR as well. The second high basal reflection is located between 650 and 800 m (Fig. 7d) and might be related to the presence of several longitudinal crevasses and moulins on the glacier surface, enhancing surface water transport to the bed. The hyperbolic reflection structure near the glacier bed between 560 and 700 m (Fig. 6) may be related to this water transport as well. Since between 300 and 850 m along the longitudinal profile ice thickness, bed slope and englacial properties are relatively constant, the spatial variations in BRPR seem to be caused by changes in the nature of the ice/bed interface. A striking feature in the ablation area is a strong BRPR peak in the center of the transverse profile. It coincides with the downstream zone of high BRPR in the longitudinal profile and is situated topographically at the valley bottom.
Basal Conditions
To allow for a more thorough assessment of basal conditions, we applied a three-layer reflectivity model (Reference Born and WolfBorn and Wolf, 1986) to calculate theoretical bed reflectivity values for ice, a basal layer and the bedrock floor, respectively. All layers were considered to be non-magnetic, homogeneous and smooth. The theoretical bed reflectivityR of this three-layer case depends on the thickness h and the conductivity a of the middle layer, as well as on the dielectric permittivity e of each layer (see Appendix). The top layer has a value of ε1 = 3.17 for temperate ice. The occurrence of gneiss, green schist and sandstone in the Sofiyskiy valley suggests a permittivity of ε3 = 9.0 for bedrock. We performed model runs for a mixture of bedrock material and water with a varying porosity Φ. According to Reference LooyengaLooyenga (1965), the bulk dielectric permittivity is
where εw 88:0 is the permittivity for water. For a high porosity 
Equation (3) represents the permittivity of a water layer, while for a low porosity (Φ ≈ 0), the permittivity is that of ice frozen to the bedrock. For 0 < Φ < 1, the bulk dielectric permittivity equals that of a till layer with varying water saturation. The bulk conductivity is also taken as a function of the porosity (Reference Gades, Raymond, Conway and JacobelGades and others, 2000). For a given porosity Φ, K will reach a maximum corresponding to a certain layer thickness h (Fig. 9). Since water has a very high dielectric constant relative to ice and rock, high values of BRPR are consistent with the presence of a water layer, given the englacial conditions as discussed above. Besides a high reflective bed, high values for BRPR could be due to an overestimated rate of attenuation with depth, which is unlikely as we notice high spatial variability in BRPR for constant ice thickness, or could be due to multiple reflection, which only occurs in overdeep-ened areas. According to the model, the maximum TZmax
occurs for a water layer with thickness h = 1:5m. However, for thinner and/or thicker water layers, a lower R occurs. It is thus quite difficult to distinguish between a thin water layer and a saturated till layer. Low values of BRPR are consistent with ice frozen to the bed, defined by a low dielectric constant.
Fig. 9. Theoretical reflectivity 11 as a function of porosity and layer thickness. The dotted line shows the path of the maximum reflectivity Rmax.
Water forms an important component of the glacier system in the ablation area during summer. Not only is surface melting a prominent feature, the glacier surface is characterized by a large number of surface meltwater channels and moulins; the latter are abundant in the downstream part of the ice tongue. At the front of the glacier, several meltwater rivers emerge from underneath the glacier ice. It is therefore likely that a well-developed subglacial drainage system exists. However, the sparsity of ice-thickness measurements outside the central flowline area do not allow for a reconstruction of the subglacial water network, based on the subglacial hydraulic potential (e.g. Reference Copland and SharpCopland and Sharp, 2001). Only the distinct peak in BRPR in the cross-sectional profile of the ablation area, situated at the valley bottom (Fig. 7b), indicates that a subglacial water channel might be present. The presence of water — identified by spatially consistent high-reflectivity zones in the RES profiles — in the ablation area is in agreement with the force-balance calculations, which indicate that significant basal sliding should occur. Values of BRPR in the accumulation area are, in view of the above discussion on englacial conditions (attenuation) and the generally low IRPN, consistent with the force-balance analysis as well, which indicates that basal sliding does not occur. According to the reflectivity model, low BRPR values point to a bedrock reflection, although thin layers of saturated till or water exhibit low reflections as well. Nevertheless, high basal reflectors are not encountered in the accumulation area, and since englacial conditions seem constant for both RES profiles, the presence of subglacial water can be excluded.
The marked difference in basal conditions between the accumulation and the ablation area, confirmed by the RES measurements and the force-balance analysis, might indicate that Sofiyskiy glacier exhibits a polythermal regime. Although, for a polythermal glacier, the ice is relatively shallow in the accumulation area (<250 m), the short summer season and extreme cold winter conditions, typical of a continental climate, might be held responsible for these conditions.
Conclusions
A 4 year field program carried out on Sofiyskiyglacier led to a proper determination of the physical characteristics of this so-called summer-accumulation-type glacier. Like many glaciers elsewhere on Earth, glaciers in this area have retreated significantly over the last century. We calculated a mean retreat rate of 18.3m a–1 for Sofiyskiy glacier over a 100 year period. Sofiyskiy is at present characterized by a relatively low accumulation, the ELA being situated at approximately 3200ma.s.l. Its surface velocity reaches a maximum near the icefall, about 110 ma–1. A force-balance analysis indicates that significant basal sliding should occur in the ablation area (downstream of the icefall), while upstream of the icefall deformational velocities are consistent with stake measurements.
RES along several profiles on the glacier not only revealed the subglacial topography, but allowed for an analysis of basal characteristics through the BRP measured from the radar sounding data. The spatial patterns of high BRP likely indicate the presence of subglacial water in the ablation area. Low basal reflections dominate the accumulation area, which reflects a lack of subglacial water. Sofiyskiy glacier may therefore exhibit a polythermal regime.
Acknowledgements
This paper forms a contribution to the research project “Glacier Fluctuations and Climatic Change in South Siberia” (Federal Office for Scientific, Technical and Cultural Affairs), contract IN/RU/006(1998–2000), and the research project “Climatic Change in the Altai Mountains” (Scientific Research Fund, Belgium (FWO)), contract G.0120.97 (1998– 2000). The authors are indebted to I. Novikov, M. Mokhtar and M. Monte for assistance in the field, K. Matsuoka for critical remarks on a first draft of this paper, and R. C. A. Hindmarsh, A. Hubbard and F. Navarro for helpful comments during reviewing.
Appendix
Three-Layer Reflectivity Model
Consider a plane-parallel absorbing film situated between two dielectric media. The total reflectivity is then given by (Reference Born and WolfBorn and Wolf, 1986)
where η =
2πh/λ
0
, h is the layer thickness and A0 = 60 m is the wavelength of the radar signal in free space. 
are the ratios of the electric vectors, which are a function of the angle of incidence and the refractive index. Indices 1–3 refer to the layer numbers. For an electric vector perpendicular to the plane of incidence, the ratio of the electric vectors is written as
where 
and 9 is the angle of incidence. n is the refractive index, which for a non-conducting medium is defined as 
where μ is the magnetic permeability and ε the dielectric permittivity. For a conducting material, the refractive index n̂ is a function of the conductivity as well, through the attenuation index
K
(Reference Born and WolfBorn and Wolf, 1986)
where a is the dielectric conductivity and ν the frequency. In Equation (4), 0 is the phase change defined by
