Introduction

Flowline models are commonly used to simulate the evolution of valley glaciers. They have been applied to many glaciers in the world to study historical changes in glacier lengths (e.g. Reference OerlemansOerlemans, 1986; Reference Huybrechts, de Nooze, Decleir and OerlemansHuybrechts and others, 1989; Reference GreuellGreuell, 1992), the response of glaciers to climate forcing (e.g. Reference KrussKruss, 1984; Reference OerlemansOerlemans, 1997) and the behaviour of a surge-type glacier (Reference BindschadlerBindschadler, 1982). Because flowline models are essentially one-dimensional, they can be applied to glaciers for which relatively little is known. Information about bedrock and surface elevations along the central flowline, in addition to the ice-flow speed and mass balance, is sufficient to model the evolution of a glacier.

Because of their simplicity, however, flowline models tend to neglect important details of glacier dynamics and mass balance. With respect to the dynamics, the ice flux is determined only from the center-line ice thickness and surface slope along the flowline. Therefore, the effect of valley shape on the variations in flow speed across the glacier and the influence of converging or diverging flow on the evolution of the ice thickness are not considered appropriately. To incorporate such effects into the computation, flowband models have been developed and used for valley glaciers (e.g. Reference BindschadlerBindschadler, 1982; Reference Wallinga and van de WalWallinga and Van de Wal, 1998; Reference Fountain, Neumann, Glenn and ChinnFountain and others, 2004). In these models, the cross-sectional geometry is parameterized as a simple shape (e.g. trapezoid, semi-ellipse or parabola), and the ice flux through the cross-section is evaluated using an appropriate shape factor (Reference NyeNye, 1965). Glacier cross-sections are often too complex to be represented as a simple shape, and the shape factor is very sensitive to the geometrical factor used for the parameterization when the glacier is deep and narrow, as is the case with many valley glaciers. Therefore, the accurate evaluation of the ice flux in flowband models has proved to be difficult, and it is one of the key issues to consider if reliable results are to be obtained. A possible solution is to use a three-dimensional (3-D) time-dependent ice-flow model (Reference Albrecht, Jansson and BlatterAlbrecht and others, 2000; Reference Schneeberger, Albrecht, Blatter, Wild and HockSchneeberger and others, 2001). Because such 3-D modelling requires a large amount of field data and computation time, it is also important to keep the model as simple as possible for it to be a useful tool.

This paper proposes a new flowline model, which better includes valley shape effects in the computation of the ice flux. This study was motivated by new datasets obtained for Rhonegletscher, Switzerland. The bedrock geometry was measured by an ice radar, digital elevation models (DEMs) were constructed from maps from 1878 and old aerial photographs (Reference Bauder, Funk and HussBauder and others, 2007), the equilibrium-line altitude (ELA) was reconstructed from temperature and precipitation records in the region (Reference ZahnoZahno, 2004), and the measurement of ice-flow speed by Reference MercantonMercanton (1916) was reanalyzed. Although flowline models have been applied to this glacier to simulate the historical front variations and to study the sensitivity to climate change (Reference Stroeven, van de Wal, Oerlemans and OerlemansStroeven and others, 1989; Reference Wallinga and van de WalWallinga and Van de Wal, 1998), the bedrock geometry used in these studies differs significantly from that obtained by the ice-radar measurement. This work therefore uses the newly obtained basal topography. Because our first attempt at flowband modelling with a trapezoidal parameterization failed to reproduce the observed surface elevation changes (Reference ZahnoZahno, 2004), the treatment of ice dynamics was improved by using a two-dimensional (2-D) ice-flow model tuned with the measured flow rate. Instead of using a shape factor, the ice flux was calculated by integrating the flow-speed field computed by a 2-D ice-flow model applied to the glacier cross-section. The cross-section at a gridpoint was determined by the bedrock elevation data. The reliable computation of the ice flux eliminated unknown parameters in the ice dynamics and enabled us to tune the mass-balance gradients to reproduce the glacier changes of the past. In this paper, we describe an improved flowline model and show the results of its application to Rhonegletscher.

Model

Flowline models solve the one-dimensional continuity equation

where h is the ice thickness, ḃ the mass-balance rate, F the vertically integrated flux per unit width and x the coordinate along the flowline (Fig. 1). To include the 2-D effects in the model, Equation (1) is integrated over the glacier cross-section:

Fig. 1.

A sketch of the model geometry. The dotted line is an actual glacier surface, and the shaded area is the modelled glacier cross-section.

The cross-sectional area, glacier width and flow velocity component perpendicular to the cross-section are denoted by a, w and u, respectively. The y and z axes are considered to be along the transverse and vertical directions in the cross-section (Fig. 1). The volume of ice flux, Q, through the cross-section, a, is computed using a 2-D ice-flow model. The flow model calculates the horizontal flow-speed field within the cross-section by solving the equations for balance of shear stresses and Glen’s flow law (Reference NyeNye, 1965):

where τ _ij and are the shear stress and strain rate, ρ the ice density, g the gravitational acceleration, S the surface elevation, A the rate factor and τ _e the effective stress. The basal flow speed u _b is introduced as a linear function of the basal shear stress τ_b :

To solve Equation (4), the stresses on the lefthand side were substituted by u using Equation (5) and the relationships and . The differential equation obtained was solved for u within the glacier cross-section by the finite-difference method. Equation (6) was used as the basal boundary condition, and the stress-free condition was imposed at the glacier surface. Starting from the solution of a linearly viscous flow, new effective stress and basal flow-speed distributions were calculated from the previously determined flow-speed field, and the calculations were repeated until the flow speed converged within 10^-5 ma^–1.

By coupling the ice-flow model with Equation (2), the evolution of the glacier surface along the flowline is determined as follows. Given a 2-D bedrock map, the shape of the glacier cross-section is known from the surface elevation at a gridpoint. The ice flux Q is obtained by integrating the flow-speed field computed for the cross-section with the 2-D flow model. The surface elevation and down-glacier slope are assumed to be constant across the glacier. With a mass balance given as a linear function of the surface elevation, Equation (2) is solved for da with a finite time-step dt and a grid length dx. The change in the cross-sectional area is then converted to a surface elevation change depending on the valley shape.

Model Application

Results and Discussion

Using the surface elevation observed in 1878 for the initial geometry, the glacier evolution was computed by solving Equation (2) with a time-step of 0.01 year and a grid size of 100 m. A time-step of 0.1 day was used only for the first 1000 steps to obtain a smooth surface. The surface elevation change over 125 years was simulated with the reconstructed ELA history, and the results were compared with the observations.

The first result obtained with the previously reported mass-balance gradients (Reference Chen and FunkChen and Funk, 1990; Reference Wallinga and van de WalWallinga and Van deWal, 1998) (experiment 1 in Table 1) showed too much thinning in the upper part of the glacier after 125 years (experiment 1 in Fig. 7), suggesting an underestimation of the snow accumulation. This was probably because the large catchment area was not incorporated in the flowline model. In the model, the glacier width was evaluated by using the surface elevation at the flowline, assuming a flat surface across the glacier (Fig. 1 and 5a). Because the surface is concave in the upper part of Rhonegletscher, marginal regions lying at higher elevations were not included and thus the accumulation was underestimated. The hypsometric curves of the modelled geometry differed significantly in the accumulation area from those obtained from the DEMs (Fig. 8). Taking the linear mass-balance gradient into account, the accumulation was underestimated by a factor of 3.0. This does not mean the surface area was underestimated by a factor of 3.0. The substantial discrepancy was due to the neglect of regions lying at very high elevations, where the accumulation is large. This is likely to be a general problem of flowline models arising from their one-dimensional nature.

Fig. 7.

The surface elevation of Rhonegletscher observed in and computed for 2000 with the sets of mass-balance gradients in Table 1.

Fig. 8.

Hypsometric curves constructed from the map for 1878 and the DEM for 2000 (solid line) and from the model (dashed line).

The correction of this problem is not simple, because it is unclear to which cross-section the neglected surface area belongs. Although the central surface elevation could be corrected by the elevation across the glacier, the redistribution of snowfall from steep side-walls to regions with a lower elevation by snowdrifts and avalanches makes the mass-balance regime complex. These processes are difficult to evaluate because mass-balance data are sparse in the upper part of Rhonegletscher. It is necessary to measure the spatial accumulation pattern in the upper part of the glacier to accurately include the effects of cross-sectional topography in the model. As a reasonable approximation, however, we increased the mass balance in the accumulation area at a uniform rate by multiplying a factor of 3.0 to the mass-balance equation (Equation (7)). This is equivalent to increasing the mass-balance gradient in the accumulation area by a factor of 3.0 (experiment 2 in Table 1). This simple procedure corrected the total accumulation by accounting for the fact that more accumulation was neglected in the higher altitude. It should be stressed that this factor was determined for Rhonegletscher based on the analysis in Figure 8 and thus it may not be readily adopted for other glaciers.

With this correction, the model results showed reasonable agreement with observations in the accumulation area, but failed to reproduce the observed terminus retreat (experiment 2 in Fig. 7). Finally, we adjusted the mass-balance gradient in the ablation area so that the modelled surface evolution (Fig. 9a and b) became similar to observations (Fig. 3a and b). The computed surface-elevation changes matched well with the observations in 1980 and 2000 except near the glacier top (0 < x < 2 km) where the ice thickness was influenced by the mass-balance correction. The observed length changes were well reproduced by the modelling results (Fig. 9c). The modelled ice volume was about 10% smaller than the observation at the beginning of the simulation (Fig. 9d) because the model neglected the marginal regions in the accumulation area. This discrepancy became smaller after 1929 as the result of the thickening near the top, and slight overestimation of thickness in other parts of the glacier. Comparison of the observed and modelled geometries in 1929, 1980 and 2000 showed that the volume was overestimated by 3.1–7.5% in the region x > 2 km. The thickening at 6 < x < 9:5 km in 1929 (Fig. 9b) might be related to the steep bed slope below the point x = 9:5 km. The longitudinal stress coupling, which was not incorporated in the model, may play a significant role in this region. It is also possible that the mass-balance regime used for this study was too simple. Two different mass-balance gradients in the ablation area and their temporal changes, as suggested by Reference Wallinga and van de WalWallinga and Van de Wal (1998), might improve the match to observations. However, it is difficult to find an observational basis for the further adjustment because mass-balance data are available only for limited periods and the data do not cover the entire glacier.

Fig. 9.

(a) The surface elevation observed in 1878, and that computed for 1929, 1980 and 2000 with the tuned mass-balance gradient set (experiment 3 in Table 1). (b) The computed surface elevation changes since 1878. (c, d) The observed and computed changes in the glacier length (c) and ice volume (d).

The future evolution of Rhonegletscher is studied by assuming three constant mass-balance conditions: the mean for the period 1994–2003 and the most negative (2003) and most positive (1978) mass-balance values in the past 50 years (Fig. 10a). Although the latter two scenarios are unlikely to happen, they are used to assess the extreme limits of the glacier shape after 50 years. They are also useful to test the stability of the model under extreme conditions. The computed changes in length after 50 years are –600 and –2600m for the 1994–2003 mean and 2003 mass balance, while a 1500m increase is expected for the 1978 mass balance (Fig. 10b). The model predicts –18%, –58% and +38% volume changes for the three mass-balance conditions, respectively (Fig. 10c). With the modest scenario, the glacier does not reach a steady state after 50 years, suggesting a longer response time to current climate conditions. An interesting consequence of further glacier retreat is the formation of a proglacial lake due to the exposure of the overdeepening near the current terminus position.

Fig. 10.

(a) ELA reconstructed for the past and assumed for the future. (b) The surface elevation and (c) ice volume computed for the future.

Conclusion

We proposed and tested a new flowline model that utilizes the cross-sectional geometry in the computation of the ice flux through the cross-section at each gridpoint. The aim of the new model was to take into account the influence of valley shape on ice dynamics accurately. We applied the model to Rhonegletscher, for which a bedrock map and several DEMs of the surface have been reconstructed recently and data on the surface flow speed are available. The first result suggested that the snow accumulation was underestimated because the flowline model neglected the marginal regions of the glacier in the accumulation area. This problem is likely to be a general weakness of one-dimensional glacier models. After the correction of the mass balance, the model reproduced the observed changes in the terminus position and the total ice volume within the error range of 15%. These results suggest that the proposed model has the potential to simulate glacier evolution accurately when adequate mass-balance data are provided. It should also be noted that the reliable computation of the ice flux reduces the number of unknown parameters, which have been used as tuning parameters in earlier flowline models.

A possible future improvement of the model is more accurate treatment of the neglected accumulation in the marginal part of the glacier. Taking into account an uneven surface elevation across the glacier may improve the performance, although it spoils the model’s simplicity and requires a detailed mass-balance model. For the reliable modelling of Rhonegletscher, it is important to know the actual mass-balance and ice-flow conditions. New field data are required to obtain accurate mass-balance gradients, the coefficients for the ELA in Equation (8) and the basal flow coefficient c.