268

Annals of Glaciology 46 2007

Evolution of Rhonegletscher, Switzerland, over the past 125 years
and in the future: application of an improved flowline model
Shin SUGIYAMA,1,2 Andreas BAUDER,2 Conradin ZAHNO,2 Martin FUNK2
1

Institute of Low Temperature Science, Hokkaido University, Sapporo 060-0819, Japan
E-mail: sugishin@lowtem.hokudai.ac.jp
2
Versuchsanstalt für Wasserbau, Hydrologie und Glaziologie (VAW), Eidgenössische Technische Hochschule,
ETH Zürich, CH-8092 Zürich, Switzerland
ABSTRACT. To study the past and future evolution of Rhonegletscher, Switzerland, a flowline model wasdeveloped to include valley shape effects more accurately than conventional flowband models. In themodel, the ice flux at a gridpoint was computed by a two-dimensional ice-flow model applied to thevalley cross-section. The results suggested the underestimation of the accumulation area, which seems to
be a general problem of flowline modelling arising from the model’s one-dimensional nature. The
corrected mass balance was coupled with the equilibrium-line altitude (ELA) change, which wasreconstructed for the period 1878–2003 from temperature and precipitation records, to run the modelfor the past 125 years. The model satisfactorily reproduced both changes in the terminus position and thetotal ice volume derived from digital elevation models of the surface obtained by analyses of old maps
and aerial photographs. This showed the model’s potential to simulate glacier evolution when an
accurate mass balance could be determined. The future evolution of Rhonegletscher was evaluated withthree mass-balance conditions: the mean for the period 1994–2003, and the most negative (2003) andpositive (1978) mass-balance values for the past 50 years. The model predicted volume changes of –18%,–58% and +38% after 50 years for the three conditions, respectively.

INTRODUCTION
Flowline models are commonly used to simulate the evolution of valley glaciers. They have been applied to manyglaciers in the world to study historical changes in glacierlengths (e.g. Oerlemans, 1986; Huybrechts and others, 1989;
Greuell, 1992), the response of glaciers to climate forcing(e.g. Kruss, 1984; Oerlemans, 1997) and the behaviour of asurge-type glacier (Bindschadler, 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 andmass balance. With respect to the dynamics, the ice flux isdetermined only from the center-line ice thickness andsurface slope along the flowline. Therefore, the effect ofvalley shape on the variations in flow speed across theglacier and the influence of converging or diverging flow on
the evolution of the ice thickness are not consideredappropriately. To incorporate such effects into the computation, flowband models have been developed and used forvalley glaciers (e.g. Bindschadler, 1982; Wallinga and Van
de Wal, 1998; Fountain 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 (Nye, 1965). Glacier cross-sections are often
too complex to be represented as a simple shape, and theshape factor is very sensitive to the geometrical factor used
for the parameterization when the glacier is deep andnarrow, 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 toconsider if reliable results are to be obtained. A possible
solution is to use a three-dimensional (3-D) time-dependentice-flow model (Albrecht and others, 2000; Schneeberger
and others, 2001). Because such 3-D modelling requires alarge amount of field data and computation time, it is also
important to keep the model as simple as possible for it to bea useful tool.
This paper proposes a new flowline model, which better
includes valley shape effects in the computation of the iceflux. This study was motivated by new datasets obtained forRhonegletscher, Switzerland. The bedrock geometry wasmeasured by an ice radar, digital elevation models (DEMs)
were constructed from maps from 1878 and old aerial
photographs (Bauder and others, 2007), the equilibrium-line
altitude (ELA) was reconstructed from temperature andprecipitation records in the region (Zahno, 2004), and themeasurement of ice-flow speed by Mercanton (1916) was
reanalyzed. Although flowline models have been applied tothis glacier to simulate the historical front variations and to
study the sensitivity to climate change (Stroeven and others,
1989; Wallinga and Van de Wal, 1998), the bedrock
geometry used in these studies differs significantly from thatobtained by the ice-radar measurement. This work therefore
uses the newly obtained basal topography. Because our firstattempt at flowband modelling with a trapezoidal parameterization failed to reproduce the observed surface elevation changes (Zahno, 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 ashape factor, the ice flux was calculated by integrating theflow-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

Sugiyama and others: Evolution of Rhonegletscher

269

Fig. 1. A sketch of the model geometry. The dotted line is anactual glacier surface, and the shaded area is the modelled glaciercross-section.

in the ice dynamics and enabled us to tune the mass-balance
gradients to reproduce the glacier changes of the past. In thispaper, we describe an improved flowline model and show
the results of its application to Rhonegletscher.

MODEL
Flowline models solve the one-dimensional continuity
equation@h@F¼ b_ 
,ð1Þ
@t@xwhere h is the ice thickness, b_ the mass-balance rate, F thevertically integrated flux per unit width and x the coordinate
along the flowline (Fig. 1). To include the 2-D effects in themodel, Equation (1) is integrated over the glacier crosssection:@a@Q¼ w b_ 
,
@t
@x
ZZ
Q¼
u ðy, zÞ dy dz:

ð2Þ
ð3Þ

a

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 thecross-section, a, is computed using a 2-D ice-flow model.
The flow model calculates the horizontal flow-speed fieldwithin the cross-section by solving the equations for balanceof shear stresses and Glen’s flow law (Nye, 1965):@yx @zx
@S
,
þ
¼ g
@x
@y
@z

ð4Þ

"_ ij ¼ Aen1 ij ,ð5Þ

where ij and "_ ij 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 ub is introduced as a linear function of thebasal shear stress b :ub ¼ cðyÞb :ð6Þ

To solve Equation (4), the stresses on the lefthand side were
substituted by u using Equation (5) and the relationships"_ yx ¼ 1=2ð@u=@yÞ and "_ zx ¼ 1=2ð@u=@zÞ. The differential
equation obtained was solved for u within the glacier crosssection by the finite-difference method. Equation (6) wasused as the basal boundary condition, and the stress-free
condition was imposed at the glacier surface. Starting fromthe solution of a linearly viscous flow, new effective stressand basal flow-speed distributions were calculated from theFig. 2. Contour maps of Rhonegletscher showing (a) the glacier bed
topography and (b) the surface topography observed in 1878. The
contour intervals are 100 m and the coordinates correspond to theofficial Swiss coordinate system. The flowline used in this study isindicated in (a) with the gridpoints located every 1 km. The
locations of the cross-profiles used for the flow model calibration
are indicated with the double circles.

previously determined flow-speed field, and the calculations
were repeated until the flow speed converged within
105 m a–1.
By coupling the ice-flow model with Equation (2), theevolution of the glacier surface along the flowline isdetermined as follows. Given a 2-D bedrock map, the shapeof the glacier cross-section is known from the surface
elevation at a gridpoint. The ice flux Q is obtained byintegrating the flow-speed field computed for the crosssection with the 2-D flow model. The surface elevation anddown-glacier slope are assumed to be constant across theglacier. With a mass balance given as a linear function of thesurface elevation, Equation (2) is solved for da with a finite
time-step dt and a grid length dx. The change in the crosssectional area is then converted to a surface elevation
change depending on the valley shape.

MODEL APPLICATION
Rhonegletscher
The proposed model was applied to Rhonegletscher, a 10 kmlong valley glacier in the Swiss Alps (Fig. 2). The historicalchanges of this glacier have been well documented since the17th century by paintings, photographs, maps and DEMs.
For this study, DEMs of past glacier geometries wereobtained on a 25 m mesh by digitizing maps constructedduring the periods 1874–82 (Mercanton, 1916) and 1928–30(Jost, 1936) and by the photogrammetric analysis of aerial
photographs taken in 1959, 1980, 1991 and 2000 (Zahno,
2004; Bauder and others, 2007). The bedrock map shown in
Figure 2a was constructed from a 2003 ice-radar survey(Zahno, 2004). We defined the flowline to be along thevalley bottom (Fig. 2a), taking the x axis from the up-glacierend of the study region. A comparison of the surface
elevations for 1878 and 2000 along the flowline (Fig. 3)

270

Sugiyama and others: Evolution of Rhonegletscher

Fig. 3. (a) The surface elevation of Rhonegletscher observed in
1878, 1929, 1980 and 2000. It decreased steadily from 1878 to
2000. (b) The surface-elevation changes since 1878.

showed that the glacier retreated by 1700 m and decreasedin thickness by >50 m at the current terminus. The surface
flow speed was measured during the period 1874–1910 bytracking stones distributed on the glacier (Mercanton, 1916).
Six flow-speed profiles across the glacier were generated
from these measurements and were used for the calibrationof the 2-D flow model. The annual mass-balance data wereavailable only for the periods 1874–75, 1885–1910 (Mercanton, 1916) and 1980–82 (Chen and Funk, 1990). The
ELAs determined from these measurements were used toreconstruct the ELA history for the region.

Ice dynamics
First, the parameters of the 2-D flow model were determined
for Rhonegletscher. In Glen’s flow law (Equation (5)), a rate
factor of A ¼ 75 MPa–3 a1 (determined for nearby Unteraargletscher (Gudmundsson, 1999)) was used, and theexponent n ¼ 3 was assumed. The linear coefficient cðyÞin the basal flow law (Equation (6)) was calibrated at six
gridpoints, x ¼ 3, 5, 8, 9, 10 and 11 km (double circles inFig. 2a), where the measured flow-speed profiles wereavailable, so that the computed surface speed profiles
reproduced the measured flow speeds (Fig. 4). The surface
elevation in 1878 was used for this calibration. The crosssection was discretized with a 20  20 mesh using thecoordinate transformation described by Blatter (1995).
Matching of the computed and measured surface flow
speeds was achieved by adjusting cðyÞ at the basal
gridpoints through trial and error. The six sets of calibrated
coefficients were interpolated along the x direction assuming that the basal flow speed was zero at the upstream
margin of the glacier and cðx, yÞ was constant below the flux
gate at x ¼ 11 km. Although cðx, yÞ was assumed to be
temporally invariant, the effect of the changing glaciergeometry on the basal flow speed was taken into account by
the stress-dependent basal flow law (Equation (6)).

Fig. 4. (a–f) The ice-flow speed measured and computed for the
cross-sections at 3, 5, 8, 9, 10 and 11 km from the glacier top. The
upper part of each panel shows the measured surface flow speeds(Mercanton, 1916) (open circles), computed surface flow-speed(thick line) and basal flow-speed (thin line) profiles. The lower partof each panel shows the computed horizontal flow-speed fields.
The contour intervals are 10 (a–d), 20 (e) and 2 m a1 (f ).

Second, the calibrated flow model was used to calculate
the volume ice flux Q at gridpoints located at intervals of
1 km as a function of the ice thickness and surface slope. At
x ¼ 6 km, for example, Qðh, @S=@xÞ was obtained as in
Figure 5b by repeating the flow model computation with thecross-sectional bedrock geometry in Figure 5a. The computed ice flux was interpolated to the intermediate gridpoints so that the flowline model could obtain Q for every
time-step using the thickness and slope at the gridpoints. The
relationship between the central ice thickness and the crosssectional area (Fig. 5c) was used to convert the areal change
to a thickness change.

Mass balance
The mass balance was provided as a linear function of thesurface elevation S:@ b_b_ ¼ ðS  EÞ:
ð7Þ@z
The ELA E was reconstructed for the period 1878–2003 fromthe temperature and precipitation records measured at
nearby meteorological stations. We used the mean valuesof two stations at Andermatt (15 km from Rhonegletscher;
1448 m a.s.l.) and Reckingen (9 km from Rhonegletscher;1332 m a.s.l.) for this purpose. The available mass-balance
observations (Mercanton, 1916; Chen and Funk, 1990)
were utilized to correlate E (m) to the summer temperatureSugiyama and others: Evolution of Rhonegletscher

271

Fig. 6. The summer temperature (June–August) (a) and annual
precipitation (b) measured at Andermatt and Reckingen; and (c) the
ELA reconstructed by Equation (8) (thin line) and observed by
Mercanton (1916) and Chen and Funk (1990) (asterisks). The boldlines are 5 year running means.

Fig. 5. (a) The cross-sectional geometry of the bedrock at 6 km fromthe glacier top (thick line). The observed (dashed line) andmodelled (solid line) surface profiles for 1878 are indicated.
(b) The ice flux through the cross-section (m3 a1 ) as a function of
the ice thickness and surface slope. (c) The cross-sectional area of
the glacier as a function of the ice thickness at the gridpoint.

(June–August) Ts ( C) and the annual precipitation P (mm a1 )by multiple regression (Fig. 6). The best-fit relationship was
given by:
E ¼ 2392:2 þ 42:348Ts  0:0936P

ðn ¼ 27; r ¼ 0:71Þ:ð8Þ

We start the modelling with the mass-balance gradients_ð@ b=@zÞreported by Chen and Funk (1990) and Wallinga
and Van de Wal (1998), and we then tune the balance
gradients to reproduce the observed glacier changes of thepast as discussed below. The sets of mass-balance gradientsused in this study are summarized in Table 1.

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 thesurface elevation at the flowline, assuming a flat surface
across the glacier (Figs 1 and 5a). Because the surface isconcave in the upper part of Rhonegletscher, marginal
regions lying at higher elevations were not included and thus
the accumulation was underestimated. The hypsometriccurves of the modelled geometry differed significantly in theaccumulation area from those obtained from the DEMs
(Fig. 8). Taking the linear mass-balance gradient into
account, the accumulation was underestimated by a factorof 3.0. This does not mean the surface area was underestimated by a factor of 3.0. The substantial discrepancy wasdue 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 onedimensional nature.
The correction of this problem is not simple, because it isunclear to which cross-section the neglected surface area
belongs. Although the central surface elevation could becorrected 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 massbalance regime complex. These processes are difficult toevaluate because mass-balance data are sparse in the upper

RESULTS AND DISCUSSION
Using the surface elevation observed in 1878 for the initialgeometry, 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 theobservations.
The first result obtained with the previously reported
mass-balance gradients (Chen and Funk, 1990; Wallinga
and Van de Wal, 1998) (experiment 1 in Table 1) showed too

Table 1. The sets of mass-balance gradients used in this study
Experiment1
2
3

Accumulation area

Ablation area
>2100 m<2100 mm w.e. m–1

m w.e. m–1

m w.e. m–1

0.00703.0  0.0070
3.0  0.0070

0.00760.00760.01500.01500.0170272

Sugiyama and others: Evolution of Rhonegletscher

Fig. 8. Hypsometric curves constructed from the map for 1878 andthe DEM for 2000 (solid line) and from the model (dashed line).
Fig. 7. The surface elevation of Rhonegletscher observed in andcomputed for 2000 with the sets of mass-balance gradients inTable 1.

part of Rhonegletscher. It is necessary to measure the spatialaccumulation pattern in the upper part of the glacier toaccurately include the effects of cross-sectional topography
in the model. As a reasonable approximation, however, weincreased the mass balance in the accumulation area at a
uniform rate by multiplying a factor of 3.0 to the massbalance equation (Equation (7)). This is equivalent toincreasing 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 thehigher altitude. It should be stressed that this factor wasdetermined for Rhonegletscher based on the analysis in
Figure 8 and thus it may not be readily adopted for otherglaciers.
With this correction, the model results showed reasonableagreement with observations in the accumulation area, butfailed 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 themodelling results (Fig. 9c). The modelled ice volume wasabout 10% smaller than the observation at the beginning of
the simulation (Fig. 9d) because the model neglected themarginal 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 otherparts of the glacier. Comparison of the observed andmodelled geometries in 1929, 1980 and 2000 showed thatthe volume was overestimated by 3.1–7.5% in the regionx > 2 km. The thickening at 6 < x < 9:5 km in 1929 (Fig. 9b)
might be related to the steep bed slope below the pointx ¼ 9:5 km. The longitudinal stress coupling, which was notincorporated in the model, may play a significant role in thisregion. 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, assuggested by Wallinga and Van de Wal (1998), might

improve the match to observations. However, it is difficult tofind 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.
The future evolution of Rhonegletscher is studied byassuming three constant mass-balance conditions: the mean
for the period 1994–2003 and the most negative (2003) andmost positive (1978) mass-balance values in the past 50 years(Fig. 10a). Although the latter two scenarios are unlikely tohappen, they are used to assess the extreme limits of theglacier shape after 50 years. They are also useful to test thestability of the model under extreme conditions. The
computed changes in length after 50 years are –600 and–2600 m for the 1994–2003 mean and 2003 mass balance,
while a 1500 m increase is expected for the 1978 massbalance (Fig. 10b). The model predicts –18%, –58% and+38% volume changes for the three mass-balance conditions, respectively (Fig. 10c). With the modest scenario, theglacier does not reach a steady state after 50 years, suggestinga longer response time to current climate conditions. Aninteresting consequence of further glacier retreat is theformation of a proglacial lake due to the exposure of theoverdeepening near the current terminus position.

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 ofvalley shape on ice dynamics accurately. We applied themodel to Rhonegletscher, for which a bedrock map andseveral 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 wasunderestimated because the flowline model neglected themarginal regions of the glacier in the accumulation area.
This problem is likely to be a general weakness of onedimensional glacier models. After the correction of the massbalance, the model reproduced the observed changes in theterminus position and the total ice volume within the error
range of 15%. These results suggest that the proposed modelhas 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 havebeen used as tuning parameters in earlier flowline models.

Sugiyama and others: Evolution of Rhonegletscher

273

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

Fig. 9. (a) The surface elevation observed in 1878, and thatcomputed 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 computedchanges in the glacier length (c) and ice volume (d).

suggestions from S. Marshall and J. Kavanaugh. We are also
grateful to M. Sharp, the Scientific Editor of this paper. This
research was partially supported by the Japanese Ministry of
Education, Science, Sports and Culture, Grant-in-Aid forYoung Scientists (Start-up), 18840002, 2006–2007.

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

REFERENCESACKNOWLEDGEMENTS
We thank the members of the 2003 ice radar campaign from
VAW, ETH Zürich, at Rhonegletscher. The manuscript was
substantially improved by insightful comments and helpful

Albrecht, O., P. Jansson and H. Blatter. 2000. Modelling glacier
response to measured mass-balance forcing. Ann. Glaciol., 31,
91–96.
Bauder, A., M. Funk and M. Huss. 2007. Ice-volume changes of
selected glaciers in the Swiss Alps since the end of the 19th
century. Ann. Glaciol., 46 (see paper in this volume).
Bindschadler, R. 1982. A numerical model of temperate glacierflow applied to the quiescent phase of a surge-type glacier.
J. Glaciol., 28(99), 239–265.
Blatter, H. 1995. Velocity and stress fields in grounded glaciers: a
simple algorithm for including deviatoric stress gradients.
J. Glaciol., 41(138), 333–344.
Chen, J. and M. Funk. 1990. Mass balance of Rhonegletscher
during 1882/83–1986/87. J. Glaciol., 36(123), 199–209.
Fountain, A.G., T.A. Neumann, P.L. Glenn and T. Chinn. 2004. Can
climate warming induce glacier advance in Taylor Valley,
Antarctica? J. Glaciol., 50(171), 556–564.

274

Greuell, W. 1992. Hintereisferner, Austria: mass-balance reconstruction and numerical modelling of the historical lengthvariations. J. Glaciol., 38(129), 233–244.
Gudmundsson, G.H. 1999. A three-dimensional numerical modelof the confluence area of Unteraargletscher, Bernese Alps,
Switzerland. J. Glaciol., 45(150), 219–230.
Huybrechts, P., P. de Nooze and H. Decleir. 1989. Numericalmodelling of Glacier d’Argentière and its historic front variations. In Oerlemans, J., ed. Glacier fluctuations and climaticchange. Dordrecht, etc., Kluwer Academic Publishers, 373–389.
Jost, W. 1936. Die seismischen Eisdickenmessungen am Rhonegletscher 1931. Denkschr. Schweiz. Naturforsch. Ges., 71,
26–42.
Kruss, P. 1984. Terminus response of Lewis Glacier, Mount Kenya,
Kenya, to sinusoidal net-balance forcing. J. Glaciol., 30(105),
212–217.
Mercanton, P.L., ed. 1916. Vermessungen am Rhonegletscher /Mensuration au glacier du Rhône: 1874–1915. Neue Denkschr.
Schweiz. Naturforsch. Ges., 52.
Nye, J.F. 1965. The flow of a glacier in a channel of rectangular,
elliptic or parabolic cross-section. J. Glaciol., 5(41), 661–690.

Sugiyama and others: Evolution of Rhonegletscher

Oerlemans, J. 1986. An attempt to simulate historic front variations of Nigardsbreen, Norway. Theor. Appl. Climatol., 37(3),126–135.
Oerlemans, J. 1997. Climate sensitivity of Franz Josef Glacier, New
Zealand, as revealed by numerical modeling. Arct. Alp. Res.,29(2), 233–239.
Schneeberger, C., O. Albrecht, H. Blatter, M. Wild and R. Hock.
2001. Modelling the response of glaciers to a doubling in
atmospheric CO2: a case study of Storglaciären. Climate Dyn.,
17(11), 825–834.
Stroeven, A., R. van de Wal and J. Oerlemans. 1989. Historic
front variations of the Rhône Glacier: simulation with an ice
flow model. In Oerlemans, J., ed. Glacier fluctuations andclimatic change. Dordrecht, etc., Kluwer Academic Publishers,
391–405.
Wallinga, J. and R.S.W. van de Wal. 1998. Sensitivity of Rhonegletscher, Switzerland, to climate change: experiments with a
one-dimensional flowline model. J. Glaciol., 44(147), 383–393.
Zahno, C. 2004. Der Rhonegletscher in Raum und Zeit: neue
geometrische und klimatische Einsichten. (Diplomarbeit, ETH
Zürich.)