Introduction
Briksdalsbreen (11.9 km2) is a western outlet glacier from Jostedalsbreen (487 km2), western Norway, located in a region with a predominantly maritime climate (Fig. 1). Systematic mass-balance measurements have not been carried out by the Norwegian Water Resources and Energy Directorate (NVE) at this glacier. Therefore, it has not received much attention in studies of the response of Norwegian glaciers to simulations of future climate changes (e.g. Reference OerlemansOerlemans, 1992; Reference SchulerSchuler and others, 2005; Reference Andreassen, Elvehøy, Jóhannesson, Oerlemans, Beldring and van den BroekeAndreassen and others, 2006).
Annual mass-balance measurements of selected glaciers in Norway show that western (maritime) glaciers in southern Norway have increased their mass since the early 1960s, while eastern (continental) glaciers show a simultaneous mass loss (Reference Andreassen, Elverøy, Kjøllmoen, Engeset and HaakensenAndreassen and others, 2005). In a future warmer climate it is estimated that many of the glaciers in Norway may disappear in the next 100–200 years (Reference JóhannessonJóhannesson and others, 2004; Reference Nesje, Bakke, Dahl, Lie and MatthewsNesje and others, 2008).
The lower steep glacier tongue of Briksdalsbreen is easily accessible and is one of the main tourist attractions in Norway (Fig. 2). The front has been measured annually since 1900/01 (Reference NesjeNesje, 2005; http://www.nve.no/bre) and retreated 859 m between 1900 and 2008 (data: NVE). As with most glaciers in Norway, the retreat during the 1930s and 1940s was significant due to a combination of high summer temperatures and low winter precipitation. In contrast, during the 1990s, the glacier experienced a significant advance caused by high winter precipitation (Reference NesjeNesje, 2005). Since 2000, the glacier has retreated significantly, mainly in response to high summer temperatures and a prolonged ablation season, and it is therefore of interest to make predictions of future frontal positions in response to potential climate scenarios.
The response of the terminus of Briksdalsbreen to mass-balance perturbations is relatively swift (Reference Nesje, Johannessen and BirksNesje and others, 1995; Reference NesjeNesje, 2005; Reference OerlemansOerlemans, 2007). Reference Laumann and NesjeLaumann and Nesje (2009) established an empirical relationship between the mass balance and the frontal response rate in order to calculate future frontal positions under different climate regimes. However, the correlation for this relationship is not very high (r 2 = 0.58) and one of their conclusions was that a dynamical glacier model coupled with a mass-balance model might give more reliable results. In this paper, an advance is made by simulating future frontal changes of Briksdalsbreen with a flowline model coupled with a surface mass-balance model.
Model Descriptions
The mass-balance model
In order to perform calculations with a glacier dynamical model coupled to climate, knowledge of both the glacier mass balance and the geometry is needed. However, only a few direct mass-balance measurements have been carried out at Briksdalsbreen (Reference PedersenPedersen, 1976). Mass-balance measurements on the adjacent Nigardsbreen, initiated in 1963, are therefore used. The distance from the head of Briksdalsbreen to the head of Nigardsbreen is <5 km. Nigardsbreen, which flows east-southeast, is 9.6 km long, covers an area of 48 km2 and ranges in elevation from 355 to 1950 m. A degree-day model for Nigardsbreen based on temperature and precipitation data from Bergen was established (Reference Laumann and NesjeLaumann and Nesje, 2009) (Figs 1 and 3). Ablation and accumulation are calculated on the basis of daily temperature and precipitation. The ablation at a certain elevation in a given time interval is calculated, following the procedure of Reference Andreassen, Elvehøy, Jóhannesson, Oerlemans, Beldring and van den BroekeAndreassen and others (2006). The temperature on the glacier is found by using a constant lapse rate. It is supposed that the precipitation will fall as snow when the temperature is below a specific value. The accumulation during the time interval is calculated by summarizing precipitation that falls as snow. The precipitation on the glacier is found by using both a vertical and horizontal precipitation gradient. The model yields the specific net-balance curve for the selected glacier. The calculation of the mass balances on Briksdalsbreen is performed with an accumulation correction factor (on average ∼20% less accumulation at Briksdalsbreen than at Nigardsbreen; for details see Reference Laumann and NesjeLaumann and Nesje, 2009, p. 224).
The ice-flow model
The ice-flow model has been applied in several previous related studies (Reference Stroeven, van de Wal, Oerlemans and OerlemansStroeven and others, 1989; Reference GreuellGreuell, 1992; Reference Oerlemans and BoutronOerlemans 1996, Reference Oerlemans1997a,Reference Oerlemansb). Therefore, only a short description is given here. It is basically a vertically integrated flowline model solved on a 100 m grid. The longitudinal profile (x = 0 at the head) is divided into a number of transverse profiles. The independent variable in each profile is the glacier thickness, H. The surface width, W, and cross-sections of a transverse profile area, S, are parameterized by supposing that the transverse section is a trapezoid, i.e.
W b is the cross-section width at the base of the glacier and λ represents the angle of the valley side with the vertical. λ is set to zero on the glacier plateau so the cross-section becomes a rectangle with width W b and depth H.
The continuity equation for this system becomes (constant density)
where t is time and b is specific mass balance. U is the mean velocity in the cross-section, i.e.
where τ is the local ‘basal driving stress’ for simple shearing flow. The first term on the right in the expression is the velocity contribution from internal deformation, and the second term is from basal sliding. The constants fd and fs vary somewhat in the literature. Here fd = 2 × 10−17 Pa−3 a−1, equivalent to A = 6.3 × 10−24 Pa−3 s−1 in Glen’s law (Reference PatersonPaterson, 1994), and fs = 1.8 × 10−12 Pa−3 m2 a−1 = 5.7 × 10−20 Pa−3 m2 s−1 (Reference Budd and JenssenBudd and Jenssen, 1975).
Data Input
The area distribution of Briksdalsbreen is documented on two official maps dated 1964 and 1994 from the Norwegian National Mapping Authority (Statens Kartverk). From these maps, we estimated the central flowline (Fig. 4) and set up a grid with cross-sections at 100 m intervals. To obtain stability in the calculations, time-steps of 0.05 years were chosen. An explicit scheme was used. The boundary condition at the glacier front is dependent on whether the flux out of the last gridpoint on the glacier will more than compensate for the melting in the next (Reference Oerlemans and BoutronOerlemans, 1996). If it does, the glacier advances to the next gridpoint. If the flux is smaller, the frontal position is calculated by supposing that the glacier surface is parallel to the glacier surface once the glacier has advanced one gridpoint, meaning that the relationship between the flux and frontal change is similar. This will prevent the somewhat stepwise image of the frontal changes that would otherwise appear owing to the discretization at 100 m intervals.
Limited information exists about the thickness of Briksdalsbreen. In the accumulation area there is a coarse ice-thickness map produced by the NVE, showing a maximum of ∼ 300 m. No ice-thickness data exist in the lower steep part of the glacier tongue. As a first estimate, the thickness along the flowline is calculated by means of
assuming a constant basal traction τ of 0.15 MPa. Here ρ, g and α are ice density, acceleration due to gravity and surface slope, respectively. For continental glaciers a value of τ = 0.1 MPa or less is commonly used, but for maritime glaciers commonly higher values are applied (Reference Oerlemans and BoutronOerlemans, 1996).
The results of the calculations are presented in Figure 5 showing, as expected, unrealistic variations on the basal profile largely due to the uncertainty in the surface slope taken from maps. Reference Echelmeyer and KambEchelmeyer and Kamb (1986) and Reference Kamb and EchelmeyerKamb and Echelmeyer (1986) showed that, due to longitudinal coupling, glacier surface slopes should be averaged over several ice thicknesses. Thus, in order to obtain more realistic values, the glacier surface was smoothed with a Gauss filter with standard deviation σ = 300 m (∼6 ice thicknesses) in the accumulation area and σ = 100 m in the steep part of the outlet glacier. As the results show, the corrections are minimal for the outlet glacier where α is large and the percentage change due to the smoothing is small. This is in contrast to the accumulation area where α is small and the changes are large. The resulting longitudinal profile is used as a first estimate in the calibration procedure below.
Calibration Procedure
In addition to the two maps from 1964 and 1994, data of the frontal position from 1963 onwards were used to calibrate the model. Using the glacier surface in 1963 and the calculated trapezoidal transverse sections as a starting point, the model was run with mass-balance data from 1963 to 2007 as climate input. The resulting surface profile at 1985 (Fig. 6) predicts a substantial glacier advance and a significant lowering of the glacier surface in the upper part of the accumulation area already by 1985. The simulation goes beyond the grid in 1985. This is not in accordance with the frontal position on a map from 1994, implying that the mass flux out of the accumulation area is too large in the simulation.
It is common to tune flow models by calibrating fd and fs (Equation (4)). Attempts to do this yielded unrealistic values for these parameters. We therefore chose to fix the values mentioned above for fd and fs and applied another approach to the calibration procedure.
The transverse sections in each gridpoint were par ameterized as trapezoids with height H in the ablation area (valley glacier) and as rectangles with height H in the accumulation area, resulting in transverse areas that are probably too large in the upper part of the glacier. This, we think, explains the excessive flow out of the accumulation area evident in Figure 6.
In order to simulate the frontal position correctly, it is important that the mass flux is correct. Studies of the measured glacier-front variations and the 1964 and 1994 maps show that, except for the lower parts of the glacier, there is little change in the glacier elevation distribution and therefore we assume that the glacier was close to steady state. This was used to find the most likely transverse areas, S, by supposing that the dynamical velocity of the glacier, U, in each transverse section is equal to the glacier balance velocity, U b, in the same transverse section. The balance velocity is
where Ai is the area on the glacier between two cross-sections and bi is the mean specific mass balance over the area. The subscripts refer to the cross-sections up-glacier from the one in question, and the denominator is the cross-sectional area (Equation (2)). If the cross-sections are parallel to each other, A 1 is the area of a trapezoid with parallel sides equal to W 1 and Wi −1 and height d 1 = 100 m. If not parallel, d 1 is the average distance between cross-section numbers i and i−1. Setting U b equal to U from Equation (4) with τ from Equation (5) yields a sixth-order equation in H, which is solved by trial and error, and the results are presented in Figure 5. As expected, the glacier thickness decreases in the accumulation area, whereas there are minor changes in the tongue, indicating a reasonable similarity to the assumed subglacial valley bottom there.
The model was then run with the new basal topography and calculated mass balance for Briksdalsbreen from 1963 to 2007 (Fig. 7). Except for the period subsequent to 2000, the model yields values that are in reasonably good agreement with the measured frontal positions. However, the model fails to simulate the significant frontal retreat after 2000. This retreat was most probably caused by excessive calving in Lake Briksdalsbrevatnet, combined with an extended ablation season (in the autumn) (Reference Laumann and NesjeLaumann and Nesje, 2009). A bathymetric map of Briksdalsbrevatnet produced before the glacier advance in the 1990s shows <10 m water depth in the western part of the lake close to the outlet and a maximum water depth of ∼20 m in the eastern part (close to the inlet) of the lake (Reference McManus and DuckMcManus and Duck, 1988). Calving is likely to contribute to a larger than normal retreat rate when the front is standing in the deepest eastern part of the lake. Based on Reference Benn, Warren and MottramBenn and others (2007) and references therein, the calving rate there can be estimated as ∼30–60 m a−1. However, our model has no calving module incorporated.
Experiments with the Model
Frontal time-lag of Briksdalsbreen
Previous studies based on simple statistical analysis of glacier-front data from Briksdalsbreen and mass-balance data from Nigardsbreen suggest that the maximum frontal response rate of Briksdalsbreen lags 3–5 years behind a mass-balance perturbation (Reference Nesje, Johannessen and BirksNesje and others, 1995; Reference NesjeNesje, 2005; Reference Laumann and NesjeLaumann and Nesje, 2009). This is a relatively short time-lag compared with other glaciers (Reference OerlemansOerlemans, 2007), so it is of interest to see whether the dynamical flowline model can reproduce this time-lag.
Initially, the model was run to steady state with an area distribution for a length of 8600 m, which was a position that the glacier oscillated around for several years between 1950 and 2000. The applied specific balance curve was found by displacing the mean specific net-balance curve (1950–2000) in such a way that the area-integrated mass balance was zero. A change was then introduced by increasing the specific net balance for 1 year (equal over the entire glacier). Figure 8 shows the frontal time-lag using two mass-balance perturbations (1 and 2 m w.e.). The front reacts immediately due to less melting, and maximum advance rate lags the perturbation by 4–5 years, which is in agreement with previous estimates (Reference NesjeNesje, 2005; Reference Laumann and NesjeLaumann and Nesje, 2009). The glacier will, however, continue to advance for about 10–12 years. Subsequently, the glacier should undergo a slow retreat to its initial frontal position. Several runs with different amplitudes (not shown here) show the same trend.
A possible explanation for this relatively short time-lag for maximum advance rate may be sought in the longitudinal profile together with Nye’s wave theory (Reference NyeNye, 1958; Reference Jóhannesson, Raymond and WaddingtonJóhannesson and others, 1989; Reference Van de Wal and OerlemansVan de Wal and Oerlemans, 1995). The theory describes how a steady-state glacier of constant width will react to a sudden mass-balance perturbation, b 1. The resulting equation is:
where C 0 is the kinematic wave velocity, D 0 is diffusion of the kinematic wave, t is time and h 1 is elevation change. Subscript 1 refers to the perturbed condition and 0 to the steady-state glacier.
The lower part of Briksdalsbreen is generally very steep, with a few stretches where the longitudinal profile is somewhat flatter (Fig. 6). This results in large variations in the velocity gradient, ∂u/∂x (or ∂C 0 /∂x). In our experiments, b 1 is equally distributed over the glacier with a duration of 1 year, so that immediately after the mass-balance perturbation ∂h 1 /∂t is determined by b 1 and the velocity gradients (because ∂h 1 /∂x = 0). Figure 9 shows how the elevation change (∂h1/∂t) develops along the glacier at different times. A mass-balance perturbation that is evenly distributed over the glacier may cause perturbations that will move down-glacier as kinematic waves due to glacier geometry. The glacier will move forward until the last of these perturbations, the one initiated highest on the glacier, has reached the glacier front.
To further illustrate the importance of the geometry along the longitudinal profile, we perturbed the steady-state glacier by two perturbations, h 1, in glacier thickness, one 1000 m from the head to illustrate what may occur in the upper part and another 5000 m from the head (covering the lower part of the glacier tongue). Both perturbations had a height of 5 m, a length of 100 m and lasted 1 year. Figure 10 shows the development of the perturbations. In the same figure, variations of ∂u/∂x for the initial stage along the profile are shown. The bump at 1000 m diffuses up- and down-glacier and disappears rapidly due to small angle and large diffusion. At 5000 m, in contrast, the perturbations may be affected by the geometry (which affects ∂u/∂x), giving transient and unstable elevation changes along the profile and in the frontal position.
Response time of Briksdalsbreen
The length response time of glaciers, τ L, is commonly defined as the time it takes to achieve of the difference between the new and old equilibrium position after a sudden and long-lasting change in mass balance, ∂b n (e.g. Reference PatersonPaterson, 1994). This is tested for Briksdalsbreen by running the model for 200 years with different mass-balance perturbations. The results are shown in Figure 11. For ∂b n = +0.3 m w.e. the response time is ∼52 years and for +0.6 mw.e. it increases to 60 years. Reference Oerlemans and BoutronOerlemans (1996) calculated the response time of Nigardsbreen (Fig. 1) as about 70 years for ∂b n = +0.4 m w.e., and that of Franz Josef Glacier, New Zealand, was 27 years for ∂b n = +0.5 m w.e. (Reference OerlemansOerlemans, 1997b). Reference BruggerBrugger (2007) calculated the response times of Rabots Glaciär and the nearby Storglaciären in northern Sweden as ∼215 and ∼125 years, respectively, considerably longer than we calculated for Briksdalsbreen. The theoretical estimate for volume time response given by Reference Jóhannesson, Raymond and WaddingtonJóhannesson and others (1989), τ v = H/−B, is 15 years for Briksdalsbreen. B is the specific net balance at the snout (B = −10 m w.e.) and H is what Reference Jóhannesson, Raymond and WaddingtonJóhannesson and others (1989) denote the ‘characteristic ice thickness’, here calculated as the mean thickness over the entire glacier (H = 150 m). This is smaller than found by our model, but can probably be explained by the fact that the volume response is commonly smaller than the length response, and that Briksdalsbreen has a height/mass-balance feedback, not taken into account by Reference Jóhannesson, Raymond and WaddingtonJóhannesson and others (1989) and Reference OerlemansOerlemans (1997a).
For negative changes in the mass balance at Briksdals breen, a comparable response time cannot be calculated. The steep glacier profile causes the front to retreat rapidly several hundred metres because the glacier is thinned and the lower part of the glacier is separated from the upper part, forming a regenerated glacier and making the response time theory inappropriate.
Climate scenarios and future glacier-front positions
The Intergovernmental Panel on Climate Change (IPCC) presented several scenarios for future greenhouse gas emissions (Reference SolomonSolomon and others, 2007 (hereafter IPCC, 2007)). These are used in general circulation models (GCMs) to simulate future climate development. The two most widely used were developed at the Max Planck Institute in Germany and at the Hadley Centre in the UK. These are denoted as ECHAM4/OPYC3 and HadAM3H, respectively.
These models are used for boundary conditions when downscaling to regional climate models (RCMs). In this study, eleven simulations using eight European regional climate models were used to downscale the IPCC (2007) A2 scenario for the 21st century to western Norway (Reference Sorteberg and AndersonSorteberg and Andersen, 2008). In addition, as in Reference Laumann and NesjeLaumann and Nesje (2009), we have applied downscaled climate scenarios from the Norwegian RegClim project (http://regclim.met.no/index_en.html). These are values for western Norway from a combination of two models (Hadley Centre and Max Planck Institute) based on the IPCC (2007) B2 scenario.
Temperature and precipitation changes are commonly provided in relation to 30 year meteorological normal periods. We have used measured values for the period 1961–90 (the present normal period) and scenario-estimated values for the period 2071–2100. The first period was given the year 1975 (mean of the 1961–90 period) and the second was given the year 2085 (mean of the 2071–2100 period). Calculations of temperature and precipitation for each year between 1975 and 2085 are carried out by supposing a linear trend from the first to the last period, and by displacing the 30 year periods 1 year at a time. The resulting temperature and precipitation values are used to calculate the annual net mass-balance values at four different elevations (Fig. 12).
There are relatively large differences between the values from the different climate scenarios. The largest differences are in the lower parts of the glacier, where the difference between the largest and smallest value is about 4 m w.e. It is also worth noting that the most negative values show that the equilibrium line will be located above the glacier from approximately AD 2060 onward. For a plateau glacier this normally means rapid down-wasting of a glacier that is responding to increased temperature by melting but has lost its accumulation area and therefore the mass balance is negative.
Of the 12 calculated net mass-balance scenarios, we chose for the dynamic model computation maximum (resulting in minimum retreat), minimum (yielding maximum retreat) and mean mass-balance values inferred from scenarios from IPCC SRES A2, in addition to the result from the IPCC (2007) B2 scenario. The dynamical model was then run from 1963 to 2085 with measured mass-balance data from 1963 to 2007 and scenario-calculated values until AD 2085. The frontal position relative to 1963 is shown in Figure 13a and b. Calving is irrelevant because the model indicates that the glacier has retreated from the lake and the model predicts increasing retreat up the steep valley from the lake inlet (Fig. 2).
The largest frontal retreat is suggested by the HIRHAM model from the Danish Meteorological Institute (DMI) with ECHAM4/OPYC3 as boundary model. The smallest retreat is predicted by a combination of the Hadley Centre and Max Planck Institute models for the B2 scenario. Close to this are the results from RCAO from the Swedish Meteorological and Hydrological Institute (SMHI) with ECHAM4/OPYC3 as boundary condition. In the middle is a model (CLM) run by GKSS Research Centre in Germany with HadAM3H (for details of the different climate models, see Reference Sorteberg and AndersonSorteberg and Andersen, 2008).
The result from the empirical formula in Reference Laumann and NesjeLaumann and Nesje (2009) is also drawn as a trend line in Figure 13a for comparison. This formula is based on a relationship inferred from area-integrated mass balance and glacier-front variations that does not take into account the glacier thickness distribution. The importance of the basal geometry is seen at two times in Figure 13a, when the glacier front has retreated ∼1000 and ∼3000 m from the 1963 position. This model predicts that the glacier should experience rapid retreat even though the climate forcing is not appreciably different from that in the preceding years. The model indicates that the glacier may retreat by a total of 2.5–5.0 km by the later part of the 21st century. A spectacular icefall may therefore disappear and the glacier will likely become a plateau glacier that gradually melts down because the equilibrium-line altitude is projected to be above the glacier surface.
Summary and Conclusions
Briksdalsbreen, a western outlet glacier from Jostedalsbreen, retreated 464 m between 1996/97 and 2008 (data: NVE). In order to simulate future frontal response to recently presented climate scenarios, a dynamical flowline model was coupled with a mass-balance model for Briksdalsbreen and Nigardsbreen forced with instrumental climate data from Bergen.
The model was used to calculate the frontal time-lag of Briksdalsbreen due to a short mass-balance perturbation. The front reacts immediately to an increase of the mass balance due to less melting, but its maximum advance rate lags the mass-balance perturbation by 4–5 years, in close agreement with previous studies. The glacier may, according to the model, continue to advance for 10–12 years. Subsequently, the glacier will experience a slow retreat rate close to its initial state.
The response time for Briksdalsbreen was calculated by running the model for 200 years with different mass-balance perturbations. For mass-balance perturbations of +0.3 and +0.6 mw.e. the response times of Briksdalsbreen are calculated to be 52 and 60 years, respectively.
To simulate future frontal response, twelve mass-balance scenarios were calculated and six selected to use in the dynamical model, choosing maximum, minimum and mean mass-balance values inferred from scenarios from IPCC (2007) SRES A2, in addition to the result from the IPCC (2007) B2 scenario. The dynamical model was then run from 1963 to 2085 with measured mass-balance data from 1963 to 2007, and scenario-calculated values to AD 2085. After the glacier retreats from the lake, the model simulation predicts continuing retreat up the steep valley from the lake inlet and a total retreat of 2.5–5.0 km by AD 2085 depending on which scenario is used. According to the simulation, a scenic icefall will disappear and the glacier will become a plateau glacier that gradually melts down.
Acknowledgements
We thank J. Ellingsen and E. Bjørseth for preparing the figures. Comments by R.LeB. Hooke and an anonymous reviewer helped to improve the clarity of the paper. This is publication no. A 223 from the Bjerknes Centre for Climate Research.