I. Introduction
It has gradually become clear that the retreat of mountain glaciers over the last 100–150 years is of a world-wide nature, e.g. Weidick (1967), Reference BjörnssonBjörnsson (1979), Reference Mayewski and JeschkeMayewski and Jesche (1979), Reference KarlénKarlén (1982), Reference MercerMercer (1982), Reference HastenrathHastenrath (1984), Reference GellatlyGellatly (1985), and Bogen and others (1988). In Figure 1, long records of glacier length are shown (Glacier d’Argentière, Rhonegletscher, and Unterer Grindelwaldgletscher in the Alps, Nigardsbreen in southern Norway, and Vatnajökull in Iceland (mean front position of five outlet glaciers of the ice cap)). The right-hand panel gives normalized records, obtained by scaling with a characteristic length L, defined as the ratio of glacier area to mean width of the tongue. Although one could think of other definitions of L, the one used here certainly corrects for differences in response due to differences in glacier geometry. Nigardsbreen is the extreme example: a large glacier with a very narrow tongue. After scaling, the glacier variations appear to be remarkably similar over the last 150 years.
It has been suggested frequently in the literature that the widespread retreat is due to the gradual increase in global temperature over the last 100 years. This warming amounts to about 0.5 K (e.g.Reference Ellsaesser, MacCracken, Walton and Grotch Elsaesser and others, 1986; Reference Jones, Raper, Bradley, Diaz, Kelly and WigleyJones and others, 1986). However, schematic statements about the relation between observed glacier retreat and slight global warming appear to have a weak foundation. To illustrate this: one could, for instance, use a statistical relation between equilibrium-line altitude (E) and temperature, based on inter-annual variations in a mass-balance series, to estimate the change in E over the last 100 or 150 years. According to an analysis by Kuhn (in press), for instance, δΤ = 0.5 K would correspond to δЕ = 40 m or so (this does not include the effect of associated increase in atmospheric counter-radiation). This is not enough to explain the retreat (see Table I). This table shows changes in equilibrium-line altitude required to explain the retreat over the last 150 years. It is based on calculations with numerical ice-flow models that fully take into account the geometry of the particular glacier. The estimate is derived from the steady-state relation between glacier length and E (which introduces an error).
With those models, several attempts have been made to simulate historic glacier variations over the last three or four centuries (Nigardsbreen: Reference OerlemansOerlemans, 1986a; Glacier d’Argentière: Huybrechts and others, in press; Rhonegletscher: Stroeven and others, in press). In these studies, the evolution of a glacier was calculated along a flow line, while forcing the equilibrium line to move up and down in proportionality with a climatic indicator like tree-ring width, (proxy) temperature and/or (proxy) precipitation. These attempts have not been successful. In particular, the substantial retreat of valley glaciers over the last 150 years does not show up in the model simulations. Of course, one can state immediately that the numerical glacier models are too simple, but it seems to be more likely that the relation between mass balance and climatic variables is the most critical factor.
Table I. Change in equilibrium-line altitude e needed to explain the difference between present-day and maximum nineteenth century glacier length
Fig. 1. Variations of some selected glaciers as measured by their length. Left: absolute scale; right: length scaled as described in the text. Data from Reference BjörnssonBjörnsson (1979); Østrem and others (1977); Reference KasserKasser (1967, 1973); Reference Kasser and HaeberliKasser and Haeberli (1979); Reference MüllerMüller (1977); Reference VivianVivian (1975); Reference HaeberliHaeberli (1985).
From a physical point of view, there are several reasons why glaciers should be very sensitive to climatic change, in particular when this change is the result of a perturbation in the radiation balance. It is important to realize that, in the ablation season, a snow or ice surface reacts in a special way. For a soil or rock surface, a change in the radiation balance will lead to higher surface temperature, and consequently to increased long-wave emission and increased turbulent-energy fluxes. When a new equilibrium is established, the net radiation balance will not be much different. In the case of a melting ice/snow surface, however, temperature cannot go up and the extra radiation is used entirely for melting! To give an order-of-magnitude estimate: an initial +1 W/m2 perturbation of the radiation budget at the surface could easily lead to an 8 W/m2 increase in emitted long-wave radiation, a 2 W/m2 increase in turbulent-heat flux, and a 9 W/m2 increase in downward long-wave radiation. (For a more detailed discussion on the climate’s response to a perturbed radiation balance and the various feed-backs involved, see e.g. Reference HansenHansen and others (1981); a further discussion is also given in the Appendix.) The associated change in surface temperature would be about 2 K. The large shift in the long-wave radiative fluxes is essentially due to the water vapour feed-back: higher air temperature leads to larger (absolute) atmospheric humidity and thus to increased emissivity. It therefore appears that for a melting glacier surface a 1 W/m2 shift of the world-wide radiation balance would lead to a 10 W/m2 change in net incoming radiation.
Another factor that makes glaciers sensitive involves the turbulent-heat flux. To first order, the turbulent-heat flux from the atmosphere to the glacier is proportional to the temperature difference between the air and ice surface. A positive perturbation of the radiation balance will lead to higher temperatures of soil or rock surrounding the glacier, and subsequently to a higher air temperature, but the glacier surface remains at the melting point. So the glacier is able to “suck” energy from its surroundings, in particular when the tongue covers a small fraction of the valley. This point has been discussed by Reference OerlemansOerlemans (1986b).
The foregoing discussion suggests that it is desirable to have at least two meteorological quantitities to act as forcing for a glacier model: air temperature and radiation budget of the glacier surface. Air temperature could be obtained from measurements (some series go back to the eighteenth century) or proxy data, but radiation is a problem. To obtain a straightforward and more universal approach, a simple climate model is used here to find perturbations of air temperature and radiation. This climate model is driven by variations in volcanic activity (reducing the amount of solar radiation impinging on the Earth’s surface) and the concentration of greenhouse gases (significant only over the last 100 years or so). In the next section, this simple climate model and the forcing functions are described. Section 3 deals with the glacier model, and in section 4 results obtained with the combined model are discussed.
2. A Simple Climate Model
The model is formulated in terms of perturbations, and has ocean mixed-layer temperature and air temperature over land (screen height, say) as dependent variables. Global average values of these quantities are considered. Compared to the ocean, heat capacity of the atmosphere and the upper soil or rock layer is very small and can be ignored when changes on a scale of decades to centuries are of interest. The simple climate model is illustrated in Figure 2. The model is driven by changes in the surface-radiation balance denoted by Q’ m (mixed layer) and Q’ l(land), giving rise to temperature perturbations Τ’m (mixed layer) and Τ’l(land). Dynamic heat fluxes between surface and atmosphere, and between air masses over land and ocean, are not explicitly modelled. Instead, the perturbation of atmospheric temperature (lower troposphere) is expressed directly in mixed layer and land temperature:
For a global study, it is natural to set p equal to the fractional coverage of the globe by ocean (i.e. about 0.7).
As illustrated in Figure 2, the atmospheric temperature plays a role in the long-wave radiation budget. The balance equation for land becomes:
The various terms have the following meaning: [1], forcing of the surface-radiation balance (short-wave or long-wave); [2], the change in long-wave radiation emitted by the surface; [3], downward long-wave radiation associated with perturbed atmospheric temperature; [4], additional downward long-wave radiation due to the water-vapour-temperature-atmospheric emissivity feed-back (here Τ a is a representative temperature for the lower troposphere, μ gives the increase of effective long-wave emissivity per degree); [5], albedo feed-back associated with increasing or decreasing ice and snow cover (Q is the annual and global mean insolation). The parameters b, γ, μ, and αl can be chosen such that results from more detailed climate models are well reproduced. This point is taken up in the Appendix.
A similar equation, but with time-dependence, can be written down for the perturbation in mixed-layer temperature:
The terms on the right-hand side are similar to those in Equation (2). The left-hand side now represents the inertia of the climate model, the associated capacity being denoted by R.
It is possible to derive explicit expressions for T’ mand T’ l. From Equations (1) and (2), we have:
where D 1 = b − α 1 Q − (1 − p) (γ + μT a).
Substitution in Equation (3) now yields:
where D m = b − α m Q − p (γ + μT a). This equation can be used for the integration in time.
The forcing functions for the glacier model now are Τ’ l (t) from Equation (4) and Z’(t), the perturbation of the radiation balance on the glacier, which is calculated from:
This expression only applies to a surface at the melting point, because the increase in long-wave emission from the surface has been omitted.
The foregoing analysis is sufficient to calculate the forcing functions for the glacier model, given proper histories of Q’. Matching of this simplified model with more sophisticated calculations on the energy balance of the climate system is best done by studying the equilibrium solutions. A discussion on this, in connection with the choice of the set of parameters used here, is presented in the Appendix.
As an illustration of the climate model consider Figure 3, in which results are shown for the case where the forcing Q’ is suddenly set at +1 W/m2 (equal over land and ocean). Part of the response in Τ’ l, is immediate — this is the direct effect of the enhanced radiation balance. The major part of the ultimate equilibrium response, however, is associated with the vapour-temperature-atmospheric emissivity feed-back and thus strongly coupled to the mixed-layer temperature through Equation (1). A similar effect is shown in the curve for Z’: a quick initial change followed by a slow increase during a few decades.
Fig. 2. Structure of the climate model. The radiative forcing is through the Qs. The climate responds by adjusting the long-wave radiative fluxes. The only prognostic variable is mixed-layer temperature, while atmospheric and land temperature are diagnoslically related to this. The example illustrates the strength of the water-vapour feed-back. A 1W/m2 forcing (square) ultimately leads to large changes in the IR fluxes (dashed squares) and a smaller change in turbulent-energy flux (thin arrow) (only shown over sea), in the case of a melting glacier surface, however, the upward IR flux does not increase.
Fig. 3. Response of the climate model to a sudden increase in the radiative forcing of +1 W/ m2.
Although not elaborated further here, the magnitude of the initial response relative to the ultimate change depends on the parameter p, of course. The initial response will be larger when p is smaller. In this way, a difference could be made between more maritime and more continental conditions. This point will be investigated in future work.
The climate model will be forced by volcanic activity. Many studies have been devoted to the volcanic activity-climate response problem, most of them being of a statistical nature (e.g. Reference Taylor, Gla–Chen and SchneiderTaylor and others, 1980; Reference Schx00F6;nwieseSchönwiese, 1984; Reference Sear, Kelly, Jones and GoodessSear and others, 1987). Reference PorterPorter (1986), in a qualitative study, has shown that in periods of increased volcanic activity the majority of glaciers tend to advance. Although the debate is still going on, most workers agree on the fact that a very large volcanic eruption causes a drop in global surface temperature of the order of 0.4 K during a year or so. Some statistical studies still show controversial results, but experiments with detailed models of radiative transfer in the atmosphere certainly support this statement (Reference Chou, Peng and ArkingChou and others, 1984). In several studies with energy-balance climate models, attempts have been made to simulate the Earth’s surface climate over the last few centuries (e.g. Robock, 1981; Reference Gilliland and SchneiderGilliland and Schneider, 1984) with reasonable success.
In these studies, the amplitude of forcing functions (volcanic activity, various solar cycles) has been optimized to give a good match with the observed temperature record. Volcanic dust is able to explain considerably more of the observed temperature variations than solar cycles.
In fact, it seems that an effect of solar variability on climate has not yet been demonstrated in a convincing way. It has now been established by the solar-maximum space program that sunspot number is related to luminosity, but that the effect is negligibly small (Reference Hoyt and McCormacHoyt, 1983). The postulated 76 year cycle associated with changes in solar diameter is highly controversial. So, in the present study, it has been decided to have only two processes contributing to the forcing function: volcanic activity, and the effect of anthropogenically produced trace gases that increase the IR emissivity of the atmosphere (carbon dioxide, methane, chlorofluorocarbons, etc.). Two series are used to estimate volcanic activity: (1) the acidity record in a Greenland ice core (Reference Clausen and DansgaardHammer and others, 1980), and (2) a series based on Lamb’s dust-veil index (DVI) for the major volcanic eruptions since 1600 A.D. (Reference LambLamb, 1970).
We follow the analysis of Reference Chou, Peng and ArkingChou and others (1984). Their calculations show that the mean radiative forcing due to the stratospheric dust load after a major volcanic eruption amounts to about −1.5 W/m2 for a year, with a tendency to somewhat larger values at high latitudes and somewhat lower values at lower latitudes (this results from the fact that aerosol reflectivity increases strongly with solar-zenith angle). In the Greenland acidity record presented by Reference Clausen and DansgaardHammer and others (1980), major volcanic eruptions are characterized by acidity levels of the order of 4 μ-equiv. H+/kg. The constant of proportionality for the radiative forcing therefore becomes 0.375 (W/m2)/ (μ-equiv. H+/kg), where the acidity is taken as a deviation from the mean value.
Lamb’s dust-veil index, although of a more subjective character than the acidity record, has the advantage that a distinction can be made between weak and strong eruptions. Weak eruptions bring little volcanic dust into the stratosphere, implying that the heating effect (through increased IR absorptivity in the troposphere) may compensate for or even exceed the cooling effect. I therefore used only eruptions with the DVI exceeding 170, and derived an annual series by using a half-life time of 1 year for the dust veil (e.g. 3 years after a major eruption the effect is reduced to 1/8). This is a commonly accepted time-scale. As a constant of proportionality for the radiative forcing, we now take 0.004 W/m2 per DVI unit. So, an eruption with DVI = 300 corresponds to a 4 μ-equiv, H+/kg acidity level.
For the greenhouse forcing, we employ the scenario (time t in year A.D.):
This includes all IR-active trace gases and broadly follows the review given by Reference WigleyWigley (1987). Here, λ. is set to 0.3 W/m2, implying the following forcing at several selected times: 1900, 0.142 W/m2 (+0.17 K); 1950, 0.233 W/m2 (+0.29 K); 2000, 0816 W/m2 (+1.00 K.). The values in brackets give the corresponding equilibrium response of the present model (see discussion in the Appendix).
The forcing functions are illustrated in Figure 4. Since the DVI is only available from 1600 A.D., the second forcing function (see Fig. 8) will be a combination of the Greenland acidity record (until 1600 A.D.) and the DVI. This requires matching, which was done by adding a constant to the DVI, in such a way that the mean radiative perturbation over the period 1600–1980 A.D. is the same for both series.
3. The Glacier Model
I use a simple continuity model applied to a stream line. It is easily derived from the continuity equation for incompressible flow, integrated over depth and width of the glacier. This type of model has now been used by various workers in slightly different forms (e.g. Budd and Jenssen, 1975; Reference Bindschadler, Harrison, Raymond and CrossonBindschadler and others, 1977; Kruss, 1984; Reference OerlemansOerlemans, 1986a).
Fig. 4. Forcing functions used in the model, expressed as perturbations of the surface-radiation balance Q’.
For a trapezoidal cross-section with varying width B and valley slope α (see Fig. 5), the prognostic equation for ice thickness H(at the centre) takes the form:
Here, Mis the surface mass balance and η =tan α.
In the basic model, sliding velocity and deformation are both directly related to the “driving stress” r, according to:
Here, h is surface elevation, g acceleration of gravity, ρ ice density, N normal load on the bed, and k1 and k 2 flow parameters.
Fig. 5. Geometry of the glacier model. The cross-section has two degrees of freedom: width and steepness of the valley walls.
Although this type of model has been shown to simulate real glacier profiles well, it is sometimes suggested that the neglect of longitudinal stress gradients would lead to large errors in the transient behaviour. To investigate this point, some tests were carried out (personal communication from W. Greuell) with a scheme in which the longitudinal stress deviator D xx was calculated along the flow line of the model glacier. The model proposed by Reference Veen van der, Veen van der and Oerlemansvan der Veen (1987) was used, which is slightly more general than the one published by Reference Alley and WhillansAlley and Whillans (1984). It should be noted that in this model the vertical variation of D xx is not explicitly considered. In the stress balance, D xx is replaced by its vertical mean value. An even simpler model for including longitudinal stress gradients has been used by Reference Budd and JenssenBudd and Jenssen (1975), but it has been demonstrated by Reference Veen van der, Veen van der and Oerlemansvan der Veen (1987) that this model is not correct.
A more detailed discussion on the effect of longitudinal stress gradients on the response of valley glaciers to climatic change will be given elsewhere (paper in preparation by W. Greuell). However, for the present discussion, the result is important: it appears that including longitudinal stress gradients has a minor (>10%) effect on the response time of a large valley glacier. In the numerical experiments described here, the simple driving-stress model described above was therefore used. In Equation (8), the overburden pressure N was set to ρgH (i.e. pressurized subglacial water is not considered). Values of k1 and k 2 used in this study are 0.7 × 10−24 Pa−3 and 10−16 s−1 m3 Pa−2, respectively. The steepness of the valley walls increases linearly with x.
In principle, the length of the model glacier is a discrete variable, because the glacier front “jumps” from point to point. A simple procedure was incorporated to obtain a smooth measure of glacier length following the method of Huybrechts and others (in press).
The formulation of the mass balance follows the equilibrium-line concept. The change in its altitude E serves as the quantity linking the mass balance of the glacier to the output of the climate model (perturbations of radiation budget and air temperature). The following formulation was employed to describe the mass balance:
In the present experiments, the following parameter values have been used: α= 0.008 year−1 and Mmax = 1.5m/year (ice depth).
Linking the equilibrium-line elevation to changes in radiation and air temperature is difficult, but it forms a crucial step in the analysis. A number of studies on the relation betweeen E and summer temperature have been carried out, but it is difficult to decide from these which part of the variations in E is actually due to air temperature and which part is due to radiation (the problem being that varying air temperature causes changes in the turbulent fluxes as well as the counter radiation). These quantities are related, of course, but long series of inter-annual variation of the radiation balance are not available to study this point. However, Kuhn’s (in press) work on long series of mass-balance measurements (of Hintereisferner) and climatic variables (of Vent) is of help here. His analysis suggests e1 = 12 m3/W and e2 = 80 m/K. (assuming that deviations in air temperature and water-vapour density are closely related). Strictly, these values apply to the central Alps only, but they certainly give an order of magnitude that is useful for the present purpose. It means that, on average, the variations in E are for about 55% due to radiative forcing (Z′) and for about 45% due to screen-height temperature (Τ′ l ]).
In the first set of experiments, the forcing will be merely through variations in E. In a second set, however, the mass-balance gradient (a) below the equilibrium line will also be changed. The reason for this is that one should expect larger balance gradients in warmer climates. The point is illustrated in Figure 6. With the assumption that melt is proportional to the integral over positive air temperatures (the degree-day method of estimating ablation is based on this, e.g. Reference Braithwaite and OlesenBraithwaite and Olesen, 1984), it appears that an increase in temperature leads to an increase in melt (black in Fig. 6) that depends on the melt in the reference case (shaded in Fig. 6). The larger the “standard amount of melting”, the more effective is an increase in temperature. The upper part in Figure 6 could represent conditions at the glacier tongue, the lower part somewhat below the equilibrium line. This implies a larger balance gradient for higher temperatures. It should also be noted that this argument applies to both the daily and the annual cycle! The formulation to deal with this in the model is:
So, below the equilibrium line, the balance gradient changes in proportion to E’ at a rate da/dE, which is set to 0.00001 (m year) −1. The steady-state geometry of the model glacier is shown in Figure 7. A 100 m grid-point spacing is used. The geometry and equilibrium-line altitude (E0 = 2700 m) were chosen to yield a typical large valley glacier for mid-latitude climatic conditions.
Fig. 6. Illustration of the effect Of a temperature increase on ablation. The increase in ablation (black) is larger when the standard ablation (shaded) is large. A temperature increase thus leads to a larger balance gradient.
Fig. 7. Equilibrium state of the standard model glacier, for E = 2700 m. The slope of the valley bed is 0.1. This glacier has a basic sensitivity of about 1400 m (length) per 100 m change in E.
4. Experiments With The Climate-Glacier Model
All integrations with the glacier model cover the period 1000–2000 A.D. Although the basic interest is in glacier fluctuations since 1600 A.D., an early start of the integration is necessary to equilibrate the model glacier. The climate model was first run with a 1 year time step to obtain annual values of the equilibrium-line altitude. Here, two cases were studied: (i) acidity + greenhouse forcing, and (ii) combined acidity/dust-veil index + greenhouse forcing, where the DVI curve has been calibrated to give the same mean forcing as the acidity series over the period 1600–1980 A.D. (see section 2).
Figure 8 summarizes the output. The upper curve shows mixed-layer temperature for case (i), the second curve the corresponding changes in equilibrium-line altitude E. Note that the mixed-layer temperature appears as a damped curve. In accordance with observations, the present simple climate model predicts a 0.5 K warming during the last 150 years. This is, of course, a direct consequence of tuning to more sophisticated models (as described in the Appendix). The lower curve shows the equilibrium-line altitude for case (ii), i.e. with Lamb’s dust-veil index. Although it is difficult to see in this figure, it will turn out that there is not too much difference between (i) and (ii).
Fig. 8. Output of the climate model. The lower two curvesgive changes in the equilibrium-line altitude. and serve asforcing to the glacier model.
Fig. 9. Model simulations. The first part of the curves (not shaded) still shows the growing phase (integration starts with no ice in 1000 A D.) Curves (a) refer to the standard case, and curves (b) to runs with changing mass-balance gradient. For one run, a close-up is given of the last two centuries.
Figure 9 then shows the response of the model glacier to these forcing functions. Curves (a) refer to the standard case, in which the equilibrium line is merely moved up and down. First of all, the simulated glacier length is in qualitative agreement with the historic record (Fig. 1): maximum extent in the early seventeenth century and half-way into the nineteenth century (with relatively little variation in between), and substantial retreat over the last 100–150 years. The amplitude, however, is too small. The difference between maximum and minimum extent is about 800 m. This converts to 0.05 non-dimensional units, and a comparison with Figure 1 reveals that this is only about half of the generally observed range.
A study of the effect of parameter variability was then carried out to see which factors would lead to a response of larger amplitude. First of all, glacier geometry appears to have some influence. The geometry chosen here (see Fig. 7) does not produce a very sensitive glacier, which is also obvious from comparing the standard sensitivity (1400 m extent per 100m change in E)to the values given in Table I. However, it was not considered useful to try many geometries and to select the most sensitive one (a better approach would be to impose the climatic signal calculated here to the glaciers mentioned in Table I. This will be done in future work).
Enhancing the sensitivity of the climate model (smaller nominal damping coefficient) enlarges the amplitude of the glacier response, but it also produces unrealistically large changes in temperature.
So, what other factors could contribute to the discrepancy between simulated and observed record? It seems that the parameterization of the mass balance forms the most critical link in the chain. As discussed earlier, there are reasons for believing that balance gradients, or at least ablation gradients, should be larger in warmer conditions. The curves (b) in Figure 9 show what happens when the mass-balance gradient is allowed to change, as described by Equation (12). Obviously, the changes in glacier response are quite substantial. While the qualitative shape of the simulated record hardly changes at all, its amplitude is almost doubled, and is therefore much closer to reality. Although this does not prove anything, it at least shows that the balance gradient may be a critical factor. The “best simulation” obtained here is also shown in a close-up in Figure 9.
It is also interesting to look more closely at the importance of the greenhouse warming. As is already obvious from Figure 4, its influence becomes comparable to the volcanic signal around 1940 or thereabouts. A run without trace-gas forcing reveals that, according to the present model, the total glacier retreat over the last 100 years is for a substantial part (about 0.5) due to the greenhouse warming. Thjs point is illustrated in Figure 10. By using only volcanic forcing (acidity in this case), the retreat since 1850 A.D. is about 600 m, but with greenhouse warming included it is about 1200 m.
A comment on time lags is also in order. The oceanic mixed layer damps and shifts the response of the climate system to any radiative forcing. Since the forcing of the glacier model is related to the direct forcing, the instantaneous response of air temperature over land, and the delayed response of the mixed-layer temperature, it is not immediately clear how the phase lag between climate forcing and glacier response will be. In Figure 11, radiative forcing, mixed-layer temperature, and simulated glacier length are plotted close together. To illustrate the behaviour of the model, two events are marked. The long period of negative radiative forcing (marked A) has a typical time-scale of 150 years or so. Here, both mixed-layer temperature and glacier length appear as smoothed replicas of the forcing. So, on this time-scale, the model is never far from equilibrium. Next, consider the positive forcing around 1860 A.D. (marked B). The maximum response of the mixed layer lags about 20 years. Although the glacier starts to retreat immediately at 1860 A.D., it does not reach a distinct minimum stand nor a maximum stand associated with the 1900 dip in mixed-layer temperature. So, on time-scales of one or a few decades, the model is grossly out of balance. These results, in fact, justify the use of the present type of model, allowing time lags between forcing, climate response, and glacier response.
Fig. 10. This illustrates the effect of the greenhouse warming. It apparently explains about half of the retreat during the last 100 years (compare vertical bars).
Fig. 11. Phase lags in the climate-glacier model. On time-scales of the order of a century (A), the model is close to equilibrium. On a scale of decades, significant shifts in the response show up (B; black dots).
Finally, some runs were carried out with steeper glaciers. The steepness affects the response time as well as the basic sensitivity of the glacier. Steeper glaciers have shorter response times, since ice velocities are larger and the height-mass-balance feed-back is weaker. In terms of glacier length, they are also less sensitive to changes in E.
First consider Figure 6 again, showing the bed profile used so far. Here, the slope of the valley floor was set at 0.1. In the additional runs, bed slopes of 0.2, 0.3, and 0.4 were used, without changing anything else. So, as a first effect, the simulated glacier will be shorter when the slope of the bed is larger. All runs were done for the “acidity + CO2 experiment”, including the variable mass-balance gradient. In Figure 12 the simulations thus obtained are put together.
As expected, the most pronounced effect is an increase of the variability on shorter time-scales. On the other hand, it is interesting to see that the largest response is found for a slope of 0.2. Apparently, the somewhat smaller basic sensitivity is more than compensated for by the shorter response time, allowing the glacier to come closer to the equilibrium state than in the case with a slope of 0.1. For slopes of 0.3 and 0.4, this effect does not show up anymore: the glaciers exhibit much smaller variations.
The fact that smaller glaciers advanced significantly at the beginning of this century, while the large ones did not, is in agreement with the results obtained here. This indicates that the model glaciers have correct response times.
Fig. 12. A comparison of simulated front positions for glaciers on beds with varying slopes. All other parameters were kept fixed.
5. Summary And Discussion
In this paper, an attempt has been made to simulate historic variations of glacier fronts by means of a simple climate-glacier model. The climate model has been driven by volcanic-activity indices, and produces perturbations of air temperature and surface-radition balance. Using these as input to a simple dynamic glacier model shows that:
the simulated record of glacier-front variations is qualitatively similar to the observed record, but the amplitude is too small;
there is not too much difference in using the Greenland acidity record or Lamb’s dust-veil index as a forcing function;
the simulation improves when mass-balance gradients are larger in a warmer climate;
low volcanic activity and increasing greenhouse warming are roughly equally responsible for glacial retreat during the last 100 years.
Maybe, the most important aspect is that the simulation presented here is better than those obtained when variations of a specific glacier are simulated with local series of temperature and precipitation (or proxy data reflecting those quantities) (Nigardsbreen: Reference OerlemansOerlemans, 1986a; Glacier d’Argentière: Huybrechts and others, in press; Rhonegletscher: Stroeven and others, in press). This indeed suggests that, on longer time-scales, the world-wide radiation balance is the decisive factor.
As mentioned earlier, on physical grounds it should be expected that the mass-balance gradient increases when the climate gets warmer. Although some authors (e.g. Reference KuhnKuhn, 1984) have considered the altitude-dependence of inter-annual variations in mass balance, a more thorough study using all available data on mass-balance gradients is needed. The same applies to the relation between meteorological parameters and the equilibrium-line altitude.
A major problem encountered here is the difference in time-scales. Many statistical relations are based on inter-annual variation within a 30 year period, say. In such a period, the geometry of a large valley glacier does not change in a significant way. However, when dealing with changes over centuries, the changing glacier geometry may play an important role. It therefore seems that a more universal theory for glacier mass balance is required, taking into account the local climatology in dependence of the degree of glacierization in a region. The development of such a theory seems to be the real challenge when it comes to relating historic (and future) glacier variations to climatic change. In addition, a natural extension of this work is to apply the simple climate model to glaciers like Hintereisferner, Nigardsbreen, Rhonegletscher, etc. These points will be taken up in future work.
Acknowledgements
I thank W. Greuell, P. Huybrechts, A. Stroeven, and R. van de Wal for the many useful discussions we had on historic glacier variations. Further helpful comments on the manuscript were given by C. Schuurmans, K. Herterich, K. Hutter, and A. Robock. I am also grateful to C Hammer, who placed the Greenland acidity data at my disposal.
Appendix
Further Discussion Of The Climate Model
To arrive at proper estimates of the parameters governing the sensitivity of the climate to “external” forcing, it is useful to carry out a further analysis of Equations (1) to (4). However, before doing this, an estimate of b (determining the long-wave emission) can immediately be obtained by linearizing the emission law of Stefan-Boltzmann around the present global mean surface temperature (about 14.5°C). This yields a value of 5.35W/(m2K) when it is assumed that the emissivity of the surface equals 1.
The next step is to calculate the equilibrium states for the simplified case where albedo feed-back and the forcing is the same over ocean and land. Multiplying Equations (2) and (3) by (1 - p) and p, respectively, yields the equilibrium equations:
Adding these equations, and using Equation (1), the following solution for Τ′ais found:
The quantity in square brackets is frequently referred to as the nominal damping coefficient. Note that for the simple stationary case considered here the perturbations of mixed-layer temperature and “land” temperature are equal to Τ′a From Equation (A-3), it can be seen immediately that climate sensitivity increases when albedo changes more with temperature (larger α) and when the water-vapour feed-back is stronger.
Values of α given in the literature vary by a factor of 2 to 3. In Oerlemans and van den Dool (1978), where a fairly detailed model for the planetary albedo was included in an energy-balance climate model, the global mean value for α turned out to be 0.006 K−1. Reference WattsWatts (1983) claimed a value of 0.0033 K−1. Results from the atmospheric GCM ice-age experiment of Reference Manabe and BroccoliManabe and Broccoli (1985) suggest α = 0.0035 K −1. Here, it was decided to use a value of 0.005 K−1 over land and of 0.003 K−1 over sea. Note that this implies an effective mean value of 0.0039 K−1.
To fix the coefficients for the water-vapour-temperature feed-back is a difficult matter, because more complicated models are normally not analysed in a way that makes possible a direct comparison with the present approach. Some informative background discussions (related to the carbon dioxide problem) have been given by Reference HansenHansen and others (1981), Reference Luther, Ellingson, MacCracken and MacCrackenLuther and ElIingson (1985), IIoffert and Flannery (1985), and also by Reference SchlesingerSchlesinger (1986). It turned out that results in broad agreement with these studies could be obtained with γ + µΤ a = 3.20 W/(m 2 К), yielding a nominal damping coefficient of 0.81W/(m 2 K) for the global equilibrium response as given by Equation (A-3). This implies that a radiative forcing of 1 W/m 2 would lead to a change in global mean surface temperature of about 1.23 K. This is within the generally accepted range.
To avoid confusion, it is useful to define the term direct radiative forcing more sharply. In connection with the carbon dioxide problem, some authors refer to this quantity as being the increase in downward long-wave radiation at tropopause level, in the case of doubling CO2 concentration say, without any feed-backs (atmospheric temperature unchanged). Others use this term for the increase in downward long-wave radiation at the Earth’s surface, which is the quantity that is needed in the present model! In the case of CO2 doubling, the former is estimated to be about 4 W/m 2 , while the latter is only about 1.5 W/m 2 (Reference HansenHansen and others (1981) gave a value of 1.1 W/m 2 , Reference SchlesingerSchlesinger (1986) a value of 2 W/m 2 ).
A final point concerns the time-scale of the climate model. It is known that the use of a one-layer model of the ocean cannot describe the penetration of a temperature perturbation with any degree of accuracy. Nevertheless, since I model the climatic signal in a very schematic way, it was not considered worthwhile to employ a multi-layer ocean model. So, the model has a single inertial component, and R is chosen such that the e-folding time-scale is 20 years (R = 6.31 × 10 8 W m −2 K −1 ). This is in line with the discussion in Reference Hoffert, Flannery, MacCracken and LutherHoffert and Flannery (1985) and Reference SchlesingerSchlesinger (1986). It should also be mentioned that the results reported in this paper are certainly not critical to the choice of this time-scale. A value of 15 or 30 years would not affect the basic conclusions.
