Introduction

Valley glaciers can serve as an important indicator of climate change due to their rapid response (Reference Oerlemans and FortuinOerlemans and Fortuin, 1992; Reference OerlemansOerlemans, 1994). The recent increase in the rate of glacier melting is assumed to be contributing significantly to a rise in the global sea level (Reference MeierMeier, 1984; Reference Dyurgerov and MeierDyurgerov and Meier, 1997; Reference Gregory and OerlemansGregory and Oerlemans, 1998). Continuous monitoring of glacier mass balance can determine the influence of climate change upon glaciers. However, data regarding mass balances are spatially and temporally limited. The records of glacier terminal positions, including those obtained from geographical features, may solve these problems. However, glaciers tend to exhibit complex responses to climate change due to their dynamic processes, so fluctuations of the termini do not directly indicate the effect of the climate conditions. Several research groups have numerically modelled various glaciers to understand their behaviour over the past few decades and to predict the changes in glacier characteristics (e.g. Reference GreuellGreuell, 1992; Reference OerlemansOerlemans, 1997; Reference Schmeits and OerlemansSchmeits and Oerlemans, 1997; Reference Zuo and OerlemansZuo and Oerlemans, 1997; Reference De Smedt and PattynDe Smedt and Pattyn, 2003).

We describe a numerical model of Koryto glacier, Kamchatka Peninsula, Russia, based on field measurements obtained in 2000, and we examine the validity of the model by comparing with records of the glacier retreat during the last half-century. Finally, we try to reconstruct the regional climate since the Little Ice Age (LIA), paying particular attention to the trends in precipitation.

Koryto Glacier

Kamchatka Peninsula extends into the North Pacific and is the largest glaciated region in northern Asia, with 448 glaciers covering ∼905 km² (Reference Muravyev and NaruseMuravyev, 1999b). Kronotsky peninsula is located on the east coast of Kamchatka (Fig. 1) and has 32 lower-altitude glaciers, some of which extend as low as 250 m a.s.l. Koryto glacier is the third-largest glacier in this region, with an area of 7.55 km² and a length of 7 km from the highest point at 1200 m to the terminus at 300 m a.s.l. (Reference Muravyev and NaruseMuravyev, 1999a). It has a simple form, with a large accumulation area and a narrow ablation area. Glaciers in Kamchatka are usually covered with debris or volcanic detritus, but the surface of Koryto glacier is clean. The daily surface flow measurements in summer 2000 indicated that basal motion occurred throughout the glacier bed (Reference Yamaguchi, Naruse, Sugiyama, Matsumoto and MuravyevYamaguchi and others, 2003). Several terminal moraines were deposited between the LIA and 1930 (Reference SolominaSolomina, 1999; Reference Solomina and CalkinSolomina and Calkin, 2003).

Fig. 1. Map of Koryto glacier (54°50′ N, 161°44′ E) with 20 m contour lines, based on a 1960 map. The results of a global positioning system (GPS) survey in 2000 were used to revise the elevations of the glacier surface. AWS marks the location of the automatic weather station. Features on inset map of Kamchatka peninsula include A: Kronotsky peninsula; B: Stopozh meteorological station; C: Esso; D: Ushkovsky ice cap.

The amplitude of the annual mass balance of Koryto glacier, which is considered to be a factor in the glacier–climate regime (Reference MeierMeier, 1984), is the second largest among 50 glaciers in the Northern Hemisphere that have been studied from 1961 to 1990 (Reference Dyurgerov and MeierDyurgerov and Meier, 1999). Thus, Koryto glacier is regarded as one of the most extreme maritime glaciers and should be very sensitive to climate change.

Description of the Glacier Model

We use a one-dimensional model with the x axis along a flowline, based largely on the model of Reference OerlemansOerlemans (1997). The gridpoint spacing is 50 m, and the three-dimensional shape is taken into account by parameterization of the cross-sectional geometry at each gridpoint. Cross-section S is trapezoidal, with the valley width at surface W _s, the ice thickness H and the average slope of the valley walls on each side of the glacier γ as follows:

(1)

The values of W _s and γ were taken from the topographic maps made in 1960 by D.G. Tsvetkov and others (personal communication, 1982). The thickness H at each gridpoint was determined from the bedrock profile obtained from a radio-echo sounding survey in 2000 (Reference Macheret, Vasilenko, Glazovsky, Moshina, Miroshnichenko, Shiraiwa, Wakatsuchi and HaraMacheret and others, 2001).

The dynamic behaviour of the glacier is described in terms of changes in H that are calculated from the mass continuity equation. Since the glacier is assumed to be composed of ice of uniform density, the conservation equation for the ice volume is

(2)

where U is the depth-averaged ice speed at x and B is the net balance at the surface. Melting of ice at the bottom is neglected. The driving stress τ _d is defined

(3)

where h is the surface elevation. The surface ice speed U _s is calculated as the sum of the internal ice deformation U _d and the basal motion U _b, i.e.

(4)

where f ₁ and f ₂ are the flow parameters. In the model, the ice deformation speed U _d on the surface is multiplied by 0.8 to calculate the depth-averaged ice deformation speed (Reference PatersonPaterson, 1994).

The balance B is assumed to vary linearly with elevation as

(5)

where α is the mass-balance gradient, which is given by different values for the accumulation and ablation areas, and ELA is the equilibrium-line altitude. Climate change is included in the model through fluctuations in the ELA.

Determination of Each Parameter

The flow parameters depend on many conditions including the bed features, debris content and crystal structure of the basal ice. The contribution of basal water to basal motion is not explicitly accounted for in Equation (3). The overall effects of these conditions are incorporated into the model as tuning parameters (f ₁ and f ₂) that vary from glacier to glacier (Reference Gregory and OerlemansOerlemans, 1998).

At Koryto glacier, the distribution of the surface flow speeds at six sites along the flowline (Reference Yamaguchi, Naruse, Sugiyama, Matsumoto and MuravyevYamaguchi and others, 2003), the ice thickness and the surface profile (Reference Macheret, Vasilenko, Glazovsky, Moshina, Miroshnichenko, Shiraiwa, Wakatsuchi and HaraMacheret and others, 2001) were measured. From the data, we determined the values of f ₁ and f ₂ to best fit the distribution of surface speeds: namely, f ₁ = 8.3 × 10⁻²⁵ Pa⁻³ s⁻¹ and f ₂ = 2.8 × 10⁻²¹ Pa⁻³ m² s⁻¹, assumed to be constant through time and space.

The mass-balance gradient α in the accumulation area was estimated from the altitudinal distribution of the net balance B measured at snow pits and crevasse walls during 1995/96 and 1996/97 (Fig. 2; Reference Shiraiwa, Muravyev, Yamaguchi, Glazirin, Kodama and MatsumotoShiraiwa and others, 1997; Reference Muravyev and NaruseMuravyev, 1999a). Both periods provided a similar value for α, so we used the average value of 7.1 m w.e. a⁻¹ km⁻¹ for the accumulation area.

Fig. 2. Altitudinal distributions of net balance in the accumulation area. The black and white triangles indicate data from 1995/96 and 1996/97 respectively. The dotted line is a linear approximation of 1995/96 data with a gradient of 7.1, and the dashed line is for 1996/97 with a gradient of 7.2.

No direct measurements of α were made in the ablation area. We therefore estimated the value of α in the ablation area as follows. First, we calculated the altitudinal distribution of ablation B _s in the accumulation area using the measured air temperatures T _a at an automatic weather station (Fig. 1; Reference MatsumotoMatsumoto and others, 1997) with a lapse rate θ of 6.5°C⁻¹ km⁻¹. The degree-day factor D varies with the surface conditions from 4.7 mm w.e. °C⁻¹ d⁻¹ for snow to 7.0 mm w.e. °C⁻¹ d⁻¹ for ice (Reference Matsumoto, Naruse, Konya, Yamaguchi, Yamada and MuravyevMatsumoto and others, 2004). We then reconstructed accumulation B _w in the accumulation area using B _w = B + B _s. Next, assuming that the altitudinal gradient of accumulation η was constant from the accumulation area to the ablation area, the altitudinal distribution of B in the ablation area can be estimated from η and B _s calculated using T _a, θ and D. Thus, we determined α in the ablation area as 11.0 m w.e. a⁻¹ km⁻¹.

The model requires fluctuations in the ELA to introduce the influence of climate change, but ELAs on Koryto glacier are only known for the years 1960, 1971, 1982, 1984, 1996 and 1997 (Reference Haeberli and MüllerHaeberli and Müller, 1988; Reference Muravyev and NaruseMuravyev, 1999a). A method of estimating the glacier ELA is therefore required.

In general, ELAs depend on summer air temperature and winter precipitation, since the ELA is determined by the balance of summer ablation and winter accumulation. The nearest meteorological station to Koryto glacier is Stopozh meteorological station, about 50 km southwest of the glacier (Fig. 1 inset, point B). The mean summer air temperature T averaged over June–August and the winter precipitation P over October–May were measured from 1941 to 1995. Here, we assumed that the glacier ELA is a function of T and P:

(6)

where P′ indicates the ratio of winter precipitation to the average winter precipitation (1941–90) and ΔT represents the deviation of the summer temperature from the average mean summer temperature (1941–90).

Approximate values of the three constant factors (k₁, k₂, k₃) were obtained empirically using four pairs of meteorological data acquired from Stopozh meteorological station and the ELA (1960, 1971, 1982 and 1984). The values obtained for k₁, k₂ and k₃ were −746, 101 and 1491, respectively. The validity of Equation (6) using these parameters is examined below.

Testing the Model

Since each parameter in the model and the factors to estimate the ELA were derived independently of the fluctuation to the glacier terminus, we test their validity by comparing the calculated glacier profile and length with its measurements. The terminus positions of Koryto glacier were measured in 1960, 1971, 1982, 1997 and 2000, and the glacier surface profiles were surveyed along the flowline in 1960 and 2000. Using lichenometry, the latest, clearly defined terminal moraine was identified at 550 m downstream of the 2000 terminus with a formation date of 1930 (Fig. 3; Reference SolominaSolomina, 1999). Because the terminal moraine is well established, we assumed that the glacier had been stable at that position for some years before 1930. We calculated the steady-state profile of the glacier when its terminus was located at the 1930 moraine, and used this profile and terminal position as the initial conditions in the calculation.

Fig. 3. Distribution of terminal moraines at Koryto glacier, modified from Reference SolominaSolomina (1999) using the GPS data surveyed in 2000.

With this initial condition and the fluctuation of mass balance in terms of the inferred ELA, we calculated the surface profile and terminal position for each year from 1930 to 2000. Note that since the ELAs were only estimated as far back as 1941 (when meteorological measurements began at Stopozh), we simply assumed that the ELAs from 1930 to 1940 were equal to the average for 1941–50. In the same way, we compensated the lack of ELA data from 1996 to 2000 with the average value for 1990–95.

We compared the calculated glacier surface profiles with the measurements. Figure 4 represents comparisons in 1960 and 2000. Although the calculated profiles are overestimated in the upper reaches in 1960 and 2000, and near the terminus in 1960, the calculated and measured surface profiles generally showed good agreement. We next compared the calculated glacier lengths with the measurements (Fig. 5); they agreed well. Thus, these two tests indicated that the model provided a plausible simulation of Koryto glacier fluctuations.

Fig. 4. Calculated surface profiles from the model (thick curves) and profile measurements (thin curves) from (a) 1960 and (b) 2000.

Fig. 5. Modelled (solid line) and measured (triangles) glacier lengths. The steps in the solid curve are a result of the 50 m grid resolution.

Discussion

Causes of the terminal retreat acceleration since the mid-20th century

The terminal retreat accelerated after the 1970s (Fig. 5). Reference Yamaguchi, Naruse, Sugiyama, Matsumoto and MuravyevYamaguchi and others (2003) argued that a cause of the acceleration was the decrease in winter precipitation. To validate this hypothesis, we conducted two simulations under the following conditions. Scenario I involves actual fluctuations in mean summer temperatures and an assumed constant winter precipitation after 1975 with the average value before 1975 (Fig. 6a), while scenario II describes actual fluctuations in winter precipitation and an assumed constant mean summer temperature after 1975 with the average value before 1975 (Fig. 6b). In Figure 6c, the simulated and measured glacier lengths are shown.

Fig. 6. Fluctuations in glacier lengths under different scenarios. (a) Scenario I: winter precipitation after 1975 assumed equal to the average value before 1975. (b) Scenario II: mean summer temperature after 1975 assumed equal to the average value before 1975. (c) Calculated glacier lengths from the model (solid curve and dotted curve) and measurements (triangles).

The results under scenario I show that the glacier terminus remains almost constant after 1982. By contrast, the results under scenario II exhibit a similar fluctuation to the measured glacier termini, which represent a rapid retreat after 1982. These simulated results support the hypothesis that the rapid retreat of the glacier terminus since the 1970s is most likely due to the decrease in winter precipitation.

The standard deviations s of P′ and ΔT from 1941 to 1995 were calculated as 0.22 and 0.70°C, respectively. The influences of the fluctuations in both parameters on the ELA were estimated using Equation (6). The fluctuation in ELA caused by s of P′ is ±164 m, whereas that caused by s of ΔT is ±71 m, which is less than half of the former value. These results indicate that the ELA of Koryto glacier has mostly been influenced by the fluctuation in precipitation during the latter half of the 20th century, which supports the previous hypothesis.

Reconstruction of the fluctuation in ELA since the LIA

Using the model with the data on terminal moraines (Fig. 3), we now estimate the evolution of the ELA from the LIA to the present. Terminal moraines were lichenometrically dated by Reference SolominaSolomina (1999). The oldest formed in 1711, which corresponds to the LIA. Although the reliability of the lichenometric method is difficult to evaluate, the accuracy is considered to be of the order ±15–20% (Reference Solomina and CalkinSolomina and Calkin, 2003). The positions of the moraines were determined with the global positioning system (GPS) in the 2000 survey, and these positions and dates were used to estimate the glacier lengths and retreat rates. Figure 7a clearly indicates that the retreat rate has not been constant since the LIA.

Fig. 7. Input data and results of the model calculation since the LIA. (a) Terminal retreat rates (solid lines) estimated from the moraine positions (triangles). (b) Modelled ELAs (dashed lines) since 1711 with the variation in the calculated glacier lengths (solid lines) and the moraine positions (triangles). (c) Decadal means of the early-summer temperature (May–June) at Esso inferred from tree-ring analyses (Reference Gostev, Wiles, d’Arrigo, Jacoby and KhomentovskyGostev and others, 1996). These data are shown as deviations from the average values (1941–90). (d) Modelled winter precipitation rate P′. The solid lines are the reconstructed decadal means of winter precipitation from the model, normalized to the average (1941–90) values.

Using the retreat rates of the glacier, we are able to estimate the fluctuation in the ELA by the inverse method; that is, we modelled various ELAs by trial and error to fit the calculated to the measured retreat rates. First, we reconstructed the steady-state profile of Koryto glacier when its terminus was located at the oldest terminal moraine (1711). We then adopted this steady-state profile as the initial condition for this calculation. To calculate variations in the ELA, we assumed that the ELA was constant during the period from one moraine to the next. Except for the initial condition, the glacier did not reach a steady state, i.e. the glacier geometry changed with time even when the terminus was located at a moraine.

Figure 7b illustrates the calculated ELAs since 1711. The changes in the calculated glacier length and the measured position of terminal moraines are also depicted. This figure indicates that the ELA was sometimes located at a lower position than in the previous period, even though the glacier had retreated since the LIA. Such a situation occurred when the retreat rate was smaller than in the previous period (e.g. 1765–1845, 1860–1917), together with the delay action due to the response time of the glacier. According to the model, the response time at Koryto glacier is approximately 120 years, which agrees with the order of the value estimated from the mean surface slope (Reference Haeberli and HoelzleHaeberli and Hoelzle, 1995).

Reconstruction of precipitation since the LIA

Equation (6) indicates that the ELA is determined by a combination of the summer temperature (ablation) and the winter precipitation (accumulation). However, in the 18th and 19th centuries, there were no temperature records near Koryto glacier. As a result, we used a proxy temperature in early summer (May–June) that was estimated from tree-ring data (Fig. 7c; Reference Gostev, Wiles, d’Arrigo, Jacoby and KhomentovskyGostev and others, 1996) at Esso in central Kamchatka (Fig. 1, point C). We then reconstructed the winter precipitation using the following procedure. Decadal averages of the proxy temperature and their deviations ΔT from the average value (1941–90) for early summer in the Esso region were calculated. The same changes were assumed for the summer temperature between Esso and Koryto glacier. Equation (6) was then used to estimate the fluctuation in P′.

Figure 7d shows the reconstructed winter precipitation, from which two tendencies were obtained after careful examination. First, the winter precipitation from the mid-19th to the mid-20th century was about 90% of the present precipitation (1941–90). Second, despite a considerable fluctuation in the precipitation, the average precipitation for the period from the 18th to the mid-19th century was almost identical to the present values.

An ice core with a depth of 212 m was retrieved from the summit of Ushkovsky ice cap (3900 m a.s.l.) in central Kamchatka (Reference ShiraiwaShiraiwa and others, 2001) (Fig. 1, point D), and analyzed for annual precipitation (Reference Shiraiwa and YamaguchiShiraiwa and Yamaguchi, 2002). We compared the modelled winter precipitation to the core data (Table 1). Each value of P′ indicates the mean precipitation over each period from one moraine to the next, normalized by the average value for the period 1941–90.

Table 1. Comparison of P′ calculated from model and by reconstructed annual precipitation from an ice core at Ushkovsky ice cap

Table 1 indicates the following tendencies. Between 1860 and 1930, P′ estimated from the model shows a precipitation decrease (0.88), and that derived from the ice core at Ushkovsky ice cap exhibits a similar slight decrease (0.93). These results suggest that a decrease in winter precipitation of about 10% since the mid-19th century is quite likely to have occurred over Kamchatka Peninsula. The values of P′ between 1765 and 1860 show a different trend: modelled P′ is almost the same as at present (0.99) while P′ at Ushkovsky ice cap indicates a slight increase (1.06). The causes of the slight difference in trends between modelled and actual data are not clear. Possibilities include the fact that the data from Ushkovsky ice cap represent annual and not just winter precipitation, but also the difference in precipitation regime between the high mountains and the low-altitude glacier.

Conclusions

A glacier model was developed using field data from Koryto glacier and tested against additional data from the glacier. The results showed that the model correctly predicted the terminus positions as well as the change in the surface profiles for the relevant period.

Employing the model in several numerical experiments, we determined that the rapid decrease in glacier length since the mid-20th century has likely been caused by the recent decrease in winter precipitation. We also concluded that the winter precipitation in the area of the glacier may have decreased by about 10% from the mid-19th to the mid-20th century. This trend is approximately consistent with the results of an ice-core analysis from Ushkovsky ice cap.

This study demonstrated that numerical modelling was effective in the reconstruction of winter precipitation variations in a valley glacier where data on glacier fluctuations and temperatures are available.