1. Introduction
Their potential for dramatic retreats has made calving glaciers an important study object in glaciology. Iceberg calving is an effective mechanism of ice loss, allowing a glacier to lose a much larger amount of ice than would occur through melting in a short period of time (Reference Van der VeenVan der Veen, 1996). Several studies have been carried out to obtain a better understanding of the calving process and to determine a ‘universal calving law’ that can explain the non-linear behaviour of any grounded calving glacier (Reference Brown, Meier and PostBrown and others, 1982);nevertheless, the fundamental physical processes that govern the calving rate are still poorly understood. Some numerical models have been formulated for individual glaciers. The most common is the water-depth model proposed by Reference Brown, Meier and PostBrown and others (1982), based on observations from 12 tidewater glaciers in Alaska, USA, and by Reference Pelto and WarrenPelto and Warren (1991) with additional data from Greenland and Svalbard. These authors found an empirical linear relationship between the annual width- averaged calving rate and the water depth at the glacier terminus. Reference Meier and PostMeier and Post (1987) pointed out that the water- depth model is not valid for the fast retreat of Columbia Glacier, Alaska, because the approach to flotation may become the controlling factor. Van der Reference Van der VeenVeen (1996) argued that the water-depth model is relevant only for glaciers that are almost in steady state. Based on the observation of the rapid retreat of Columbia Glacier and several other grounded glaciers, he suggested that the position of the calving front is controlled by water depth and ice thickness at the glacier terminus. Van der Reference Van der VeenVeen (1996) presented the height-above-buoyancy model, in which if the frontal thickness becomes less than a critical thickness H c, the glacier terminus breaks off, which leads to the frontal thickness exceeding the flotation by an amount H0. The critical thickness is given by
where ρ w is the water density, ρ i is the glacier ice density and d is the water depth at the terminus. H0 represents the minimum thickness above the flotation thickness, which is about 50 m on Columbia Glacier. The minimum thickness may be lower for smaller glaciers (Reference Van der VeenVan der Veen, 2002). Reference Vieli, Funk and BlatterVieli and others (2001) proposed a modified flotation criterion that defines the minimum thickness H0 as a small fraction q of the flotation thickness. The critical thickness is:
The factor q is smaller for glaciers that are thinner or have fewer crevasses. Reference Vieli, Funk and BlatterVieli and others (2001) applied this modified flotation criterion in a time-dependent twodimensional (vertical plane) numerical flow model which includes a water-pressure dependent sliding law. The model is capable of simulating a cycle of slow advance and rapid retreat across a basal depression. However, the modelling effort by Reference Vieli, Funk and BlatterVieli and others (2001) was concentrated on only two transient scenarios of advance and retreat over a basal depression. Therefore, a number of questions remain. First, is their model adequate for a full cycle of glacier length variations? Secondly, is it necessary to use a detailed treatment of ice flow to produce the characteristic nonlinear behaviour of calving glaciers or can any model with a calving rate related to water depth produce the same behaviour? Lastly, since their study does not consider the equilibrium states of their model, can any terminus position on a bed profile with a basal depression represent a steady state or not?
Neither the water-depth model nor the flotation model is supported convincingly by theory and so, as yet, interpretation of the data from different glaciers has not proved either model. The current study was carried out to investigate the role of the different calving parameterizations in the model predictions and to address the questions mentioned above. For this purpose, we compared three models for a full cycle of growth and decay. The cycle starts with an ice-free condition, where the glacier terminates on land, extends with further calving into water and then retreats. A minimal model based on mass conservation for the entire glacier was compared with two numerical one-dimensional ice-flow models for calving glaciers. The two ice-flow models predict the ice-thickness distribution and the ice surface velocity along the flowline in one-dimensional space and time. The first is the water-depth model, where calving is linearly related to water depth at the terminus, and the glacier terminus position changes in response to the balance between the ice velocity and the calving rate (Reference Meier and ReehMeier, 1994, 1997). The second ice-flow model is the flotation model, which applies the modified flotation criterion, Equation (2). In contrast to the water-depth model, the calving rate in the flotation model is an output quantity of the model (Reference Van der VeenVan der Veen, 1996; Reference Vieli, Funk and BlatterVieli and others, 2001). In the minimal model, the mass budget (surface mass balance minus calving flux) and the water depth at the glacier terminus determine the variation of glacier length, irrespective of the ice surface profile and the velocity field. In this model, calving rate is defined as a linear function of the water depth (Reference Brown, Meier and PostBrown and others, 1982; Reference Funk and RöthlisbergerFunk and Rothlisberger, 1989; Reference Pelto and WarrenPelto and Warren, 1991; Reference Björnsson, Pálsson and GudmundssonBjörnsson and others, 2000).
This paper first describes the three models and their calving schemes and then compares them. We have performed a systematic study of the equilibrium states of each model glacier for two different formulations of the mass balance. The first is a uniform accumulation rate over the entire glacier. This simple formulation allows the glacier to reflect the role of the bed topography without any surface mass-balance feedback. The second formulation involves the application of a common surface mass-balance function for temperate glaciers, which is a linear function of altitude. Two idealized bed geometries are used in the model experiments. One is a linear downward-sloping profile, which allows us to observe the difference between the models for the simplest possible geometry. The other bed has a profile with a Gaussian-shaped bump superimposed on a linear downward slope; this resembles a typical bed geometry for tidewater glaciers. This paper further presents an analysis of the response of each model glacier to a sudden climate change.
2. Description of Models
Three numerical models were constructed to examine the non-linear behaviour of tidewater glaciers.
2.1. Minimal model
In this model, the evolution of the glacier is obtained from the conservation of mass for the entire glacier. The glacier may gain mass as a result of precipitation in the accumulation area or lose mass due to melting and frontal calving. If we assume a constant ice density, mass conservation is equivalent to conservation of volume:
in which V is the ice volume, B tot is the surface mass balance and C is the calving flux. This model and its properties are described in more detail by Oerlemans and Nick (2006). The model defines the mean altitude of the glacier surface, h m, instead of specifying any surface profile:
in which b 0 and b f are the bed elevations at the top and terminus of the glacier respectively. H m and Hf are the mean glacier thickness and the frontal thickness. A simple assumption was adopted in this model, namely that the mean ice thickness H m and the frontal thickness H f are proportional to the square root of the glacier length when the glacier terminates on land or in shallow water. The dependence of the ice thickness on the square root of the length L is based on Weertman’s theory for perfectly plastic ice sheets (Reference WeertmanWeertman, 1961) and extensive numerical glacier model calculations (Reference OerlemansOerlemans, 2001). The frontal thickness of the glacier when it terminates in deep water is defined by using the flotation criterion. The frontal thickness exceeds the flotation thickness by a fraction q = 0.15 of the flotation thickness (Reference Vieli, Funk and BlatterVieli and others, 2001). Hm and Hf are expressed as:
where αm and αf are constants evaluated from simulations with a numerical glacier model including deformation and sliding (Reference OerlemansOerlemans, 2001).
The calving rate is assumed to be a linear function of the water depth at the glacier terminus, with a constant of proportionality c = 3.5 (Reference Brown, Meier and PostBrown and others, 1982; Reference Funk and RöthlisbergerFunk and Rothlisberger, 1989; Reference Pelto and WarrenPelto and Warren, 1991; Reference Björnsson, Pálsson and GudmundssonBjörnsson and others, 2000). Although such a relationship may not be applicable for rapid changes (Reference Van der VeenVan der Veen, 1996), it can provide a good first-order estimation of the calving mass loss. Therefore, we described the calving flux as:
The ice volume V is expressed in the model as the mean ice thickness H m times the glacier length L:
From Equations (3), (5) and (8), it follows that:
which can be easily numerically integrated.
2.2. Flotation model
A one-dimensional numerical ice-flow model was used to calculate the surface evolution and ice surface velocity along the flowline (Reference OerlemansOerlemans, 2001). For calculating the ice velocity, the shear stress is related to the strain rate according to Glen’s flow law for plane shear (Reference GlenGlen, 1955). The ice velocity is expressed as the vertical mean velocity, which is governed by both basal sliding and ice deformation. A ‘Weertman-type’ sliding law, supported by the experimental work of Reference Budd, Keage and BlundyBudd and others (1979), is applied in the model:
where H is the ice thickness, S b is the basal stress, P is the basal water pressure, g is the gravitational acceleration, and is an empirical parameter. By assuming that P is proportional to the ice thickness and that the basal stress is equal to the driving stress for plane shear, the vertical mean ice velocity U can be expressed as:
in which f d and f S are the flow parameters (Reference OerlemansOerlemans, 2001) and S d is the driving stress, which is proportional to the ice thickness and the slope of the ice surface ∂h/∂x:
The evolution of the ice surface can be determined by the vertically integrated continuity equation for incompressible ice (Reference Van der VeenVan der Veen, 1999)
where t is time, x is distance along the flowline and B is the surface mass balance. Substituting the ice velocity U into the continuity equation gives:
where b is the basal elevation and D the diffusivity:
In this model, the continuity equation (14) is solved on a numerical grid along the flowline. First, using central differencing, the diffusivity is calculated at each gridpoint. Then the ice flux HU is computed at intermediate gridpoints by interpolating values for the diffusitivity, and finally the new ice thickness is calculated using a forward scheme for the time derivative in the continuity equation (Reference Van der VeenVan der Veen, 1999). There are 250 gridpoints along the flowline, initially with a uniform distance of ∆x = 200 m apart. This distance changes with every time-step as a new grid is defined to fit the new glacier length (120 < ∆x < 240). The time-step is 0.002 year, which is small enough to satisfy the numerical stability criterion (Reference SmithSmith, 1978) for different values of ∆x.
The model domain contains boundaries at the upstream and downstream ends of the glacier. There is no ice flux into the model from the upstream boundary; therefore, the ice velocity at the first gridpoint is zero and the ice thickness at this gridpoint is extrapolated from the neighbouring points. Two boundary conditions at the downstream end are applied. First, the model assumes that the ice velocity increases linearly at the last three gridpoints; hence, the ice velocity at the terminus is extrapolated from the adjacent calculated values. Second, the terminal thickness cannot be lower than a given limit H c (Fig. 1). At each time-step, the position of the terminus is moved to the point where the ice thickness is equal to H c. The terminus may move backward (retreating (Fig. 1b)) or forward (advancing (Fig. 1c)) or stand still (equilibrium state). The actual position of the terminus is computed by interpolating between values of two neighbouring gridpoints with ice thicknesses larger and smaller than H c. Thereafter, new gridpoints are defined to fit the updated glacier length, and the ice thickness is calculated for the new gridpoints by using linear interpolation. The modified flotation criterion (Reference Vieli, Funk and BlatterVieli and others, 2001) is applied in the model if the glacier terminates in water. The given thickness H c is the same as the frontal thickness H f that is prescribed in the minimal model.
In the flotation model, the glacier length variation is mainly a consequence of the glacier thinning or thickening due to a surface mass-balance change. Calving rate, U c, is defined as the volume of ice that breaks off per unit time per unit vertical area from the glacier terminus (Reference PatersonPaterson, 1994). Assuming that the calving rate includes all the ice-loss processes at the terminus, it is expressed as:
where U f is the ice velocity at the terminus. Thus, the calving rate is a result of the glacier dynamics and the basal geometry (Reference Van der VeenVan der Veen, 1996).
2.3. Water-depth model
The third model is the same numerical ice-flow model as described in section 2.2, but with different boundary conditions at the downstream end. The glacier terminus advances (retreats) whenever the calving rate is less (greater) than the ice velocity at the terminus. For a glacier that terminates in water, the calving rate is a linear function of water depth; the calving rate is zero when the glacier terminates on land:
Another boundary condition is imposed to avoid a frontal thickness less than the flotation thickness. At each time-step, after determining the glacier length, the frontal thickness is set to a prescribed value H f (Equation (6)).
In this model, the calving rate is the controlling process for the glacier length variation. The position of the terminus is obtained from the ice velocity at the terminus, U f, minus the calving rate:
This equation describes the response of the glacier flow to increasing calving rates as the glacier terminus advances into deeper water (Reference Meier and ReehMeier, 1994).
3. Model Input
The model calculations were done for two idealized surface mass-balance formulations and two simplified basal topographies. For simplicity, we assumed a constant width for the entire glacier.
3.1. Surface mass balance
We simulated equilibrium states of a glacier for two different surface mass-balance formulations. In the first case, we assumed that the accumulation rate a is uniform over the entire glacier. Therefore, the annual mass gain at each point is the same, and equal to the length of the gridpoint ∆x times the accumulation rate:
For the minimal model, in which no surface profile is specified, we defined an annual surface mass gain for the entire glacier:
In the second case, the mass balance depends linearly on height. The height dependence of the mass balance was described as:
where ß is a constant balance gradient, h(x) is the surface height and E is the equilibrium-line altitude (ELA). For very high locations on the glacier surface, a maximum value for B was assumed. In the minimal model, the mean altitude of the glacier surface was used to compute the total surface mass balance:
Different mass-balance forcing was imposed in the models by varying a and E as a function of time.
3.2. Basal topography
The models were run for two different bed topographies. To investigate the properties of the model in the absence of any effect of a non-linear basal geometry, we ran the models first for a linear slope,
where b(x) is the bed elevation and x is the distance along the flowline. The bed elevation at the upstream end was set at b o = 220 m and the negative bed slope s = –0.015 (solid line in Fig. 2).
To represent a basal topography for a tidewater glacier with a submarine undulation at its terminus, a Gaussian- shaped profile superimposed on the linear slope profile was used in the models:
in which γ λ 340 m, x s = 40 km and σ = 10 km are the amplitude, location and the width of the bump respectively. The geometrical set-up is shown in Figure 2.
4. Results
4.1. Steady states
This section presents the model results for two types of bed topography. With each bed profile, the models were run for the two surface mass-balance formulations as described in section 3.1.
4.1.1. Linear bed
Figure 3 shows the equilibrium glacier length vs the accumulation rate for the three models on the linear bed. The equilibrium glacier length is obtained by increasing and decreasing the accumulation rate a very slowly (compared to the glacier response time). Each model run started from the ice-free condition (a = 0). A small increase in a forces the glacier to grow and reach a steady state. Since the surface mass balance is always positive, the glacier can only reach a steady state in water where the additional surface mass is compensated by the calving flux at the glacier front. The glacier tends to advance and obtains a larger equilibrium length as the accumulation rate rises. The minimal model (solid line) and the water-depth model (dotted line) behave alike because the formulation for the calving rate is the same in the two models. In contrast, the flotation model (dashed line) shows a larger equilibrium length for a < 0.8ma–1. However, for larger values of a, the sensitivity of the equilibrium length becomes much smaller than in the other two models. In the flotation model, if the glacier front is in shallow water, the frontal thickness is enough to satisfy the flotation criterion; therefore, calving hardly occurs and the glacier can grow further to reach a greater equilibrium length. In the minimal and the water-depth models, the calving rate is a linear function of the water depth and produces a calving flux, even in shallow water, causing an equilibrium state at a shorter length. On the linear bed, to become longer, a glacier must terminate in deeper water (Fig. 2). In the flotation model, when the glacier enters deeper water, a larger frontal thickness is required to fulfil the flotation criterion. In other words, if the frontal thickness is not sufficiently large to satisfy the flotation criterion, the terminus position jumps back to the location where the thickness is large enough. In this way, a glacier may lose a huge amount of ice and become unable to advance even for high accumulation rates (right end of the dashed line in Fig. 3). Decreasing a causes the glacier in all three models to retreat and reach the same equilibrium length as in the case of advance.
In the second experiment, the height-dependent surface mass balance (Equation (21) or (22)) was applied in the models. By gradually varying the ELA, E, the models simulate a full cycle of glacier initiation, advancement into deep water and retreat. Figure 4 depicts the equilibrium glacier length vs E for each model. The glacier starts to grow and reaches steady state as soon as E becomes lower than the bed elevation at the upstream end (E = 220 m). The equilibrium terminus position depends on the basal slope and the frontal boundary condition (e.g., on a steeper bed slope, a glacier may reach steady state while terminating on land). Lowering E causes the glacier to achieve a larger equilibrium length. In the flotation model, the sensitivity of the glacier to climate change is too small. Therefore, the flotation model may fail to simulate fast retreat of glaciers on a smooth basal topography (Svalbard). Again, the glacier is not capable of growing into deeper water, as a result of the model’s flotation criterion (right end of dashed line in Fig. 4). Increasing E forces the glacier to retreat and reach a smaller equilibrium length. Above a critical value E, slightly different for different models, the glacier reaches the lower branch of the steady states and vanishes entirely. The hysteresis seen in each diagram is the well-known nonlinearity due to the height–mass-balance feedback.
4.1.2. Non-linear bed
Figure 5a shows the equilibrium glacier length vs the accumulation rate on the non-linear bed profile (shown in Fig. 2). This experiment reveals the response of a glacier to the existence of a basal depression excluding the non-linear height–mass-balance feedback. For a range of values of a (~ 0 < a < 1.5) a branch of equilibrium states is produced by each model. This branch starts where the glacier terminus enters water, and ends just before the glacier terminus reaches the deepest point of the basal depression (Fig. 5b). For larger values of a, the glacier passes the deepest point and moves onto the upward bed slope. Here, the calving flux decreases significantly as the frontal water depth becomes smaller. The increasing mass input and decreasing calving flux allow the glacier to grow and pass the bump; consequently, the glacier does not reach steady state while the bed slopes upward. After passing the bump, the water depth against the calving face increases and results in a larger calving flux. This inhibits further growth of the glacier, and a new balance is established. Figure 5b presents the surface profiles of the flotation model (dashed line) and the water-depth model (dotted line) for a given value of a in two equilibrium branches. The horizontal dashed and dotted lines under the surface profile illustrate the possible equilibrium terminus positions for the flotation and water-depth models respectively. Once again, a larger value of a provides a longer equilibrium length in the minimal and water-depth models (solid and dotted lines in Fig. 5a), whereas in the flotation model the equilibrium length remains constant due to the model’s calving criterion (dashed line). Figure 5a demonstrates that for a glacier that has advanced into deep water by passing over a basal depression and has reached a new branch of steady state, a considerable lowering in the accumulation rate is required to move the glacier back to the lower branch of steady states. This hysteresis seen for each model is the direct result of the existence of a basal undulation and the dependence of the calving flux on the water depth.
In the next experiment, the height-dependent mass- balance formulation was applied to the non-linear bed. Figure 6a represents the equilibrium glacier length vs E. The solution for each model contains two regions of hysteresis (Reference Oerlemans and NickOerlemans and Nick, 2005). For the non-linear bed profile, the dependence of calving rate on water depth produces one hysteresis which is embedded in another hysteresis caused by the height–mass-balance feedback. There are two branches of equilibrium solutions, in which the terminus of the glacier is always situated on the downward slopes of the bed profile. The lower branch provides the terminus positions on the downward slope just before the deepest part of the bed depression, and the higher branch covers the downward slope after the bump into deeper water (Fig. 6b). Once more, the flotation model fails to provide a glacier advancing into deeper water.
4.2. Response to sudden climate change
We studied the effect of a sudden climate change on the modelled glaciers by shifting E after a steady state had been reached. This experiment was carried out only with the height-dependent surface mass-balance formulation and the non-linear bed. A decrease or increase in E causes an advance or retreat respectively (Fig. 7b).
4.2.1. Advancing scenario
In all three models, the glacier approaches a steady state for E = 220 m in the lower branch of equilibrium solutions. A large increase in the mass input (E = 170 m) forces the glacier to advance in deep water, to pass the basal depression and to reach a new steady state (Fig. 7a). In the first phase of advance, before the deepest part of the basal depression (point A in Fig. 2) is reached, the glacier advances slowly since the calving rate increases when the terminus enters deeper water. In this phase, the glacier develops more slowly in the flotation model (dashed line) than in the water-depth model (dotted line), as the flotation criterion requires a large frontal thickness when the glacier terminates in deeper water.
After passing the deepest part, the glacier moves into shallower water and the calving rate decreases. This allows the glacier to advance rapidly and re-establish a balance on the other side of the bump. In this second phase, the glacier in the flotation model becomes faster than the glacier in the water-depth model because the frontal thickness is already thick enough to satisfy the flotation criterion for shallower water. After passing the highest part of the basal bump (point B), the glaciers in the minimal model and the water-depth model obtain approximately the same equilibrium length, whereas in the flotation model the glacier cannot advance into very deep water and reaches a shorter equilibrium length.
The minimal model (solid line) produces the slowest advance because the length variations in the model are calculated from the mean height of the glacier surface, i.e. the length variation depends on the entire surface variation. In this model, the glacier must obtain a sufficiently large mean thickness to be able to advance. However, in the two ice-flow models, the glacier tongue can advance before the entire glacier surface is influenced by the climate change.
4.2.2. Retreating scenario
Shifting E to a higher value of 440 m after 12000 years causes the glacier to retreat and disappear entirely. First, the glacier retreats slowly up the basal bump in the shallow water. After it passes point B and enters deep water, the calving rate increases, resulting in a rapid retreat. In the first phase, while retreating up the bump, the glacier in the flotation model retreats more slowly than the glacier in the other two models because the frontal thickness is still large enough to satisfy the flotation criterion (Fig. 7a). The minimal model provides the fastest retreat in this phase, because in the flowline models it takes time before the perturbed mass balance higher up the glacier affects the glacier front. In the minimal model, this occurs very quickly because a decrease in volume, no matter where this is generated, implies that the glacier has to shorten immediately. After the glacier front passes the highest point and retreats into deeper water, the calving rate increases and forces the glacier to retreat rapidly and disappear in all three models.
5. Conclusion and Discussion
We have compared three calving models with simplified geometries and made the following observations. First, the results suggest that any model in which the loss of ice at the glacier front increases with water depth shows qualitatively the same behaviour when a submarine undulation is present in the basal topography. Secondly, in all three models, the branching of the steady states depends significantly on the bed geometry, i.e. every basal depression can lead to a new hysteresis. The equilibrium branches include only the terminus positions on the downward slopes of the bed; no terminus position on the upward slope can represent a steady state. Lastly, although the flotation model is capable of simulating retreat and advance of the glacier across the bed depression, it does not allow development of a large glacier terminating in very deep water, i.e. this model is not adequate for a full cycle of glacier length variations. One possible explanation is that the calving criterion of the flotation model causes the glacier to lose mass in an undesirable way when it advances into deeper water.
We presented the minimal model, which is effective and computationally simple, as a learning tool. This model confirms the concept of rapid retreat and slow advance of a grounded calving glacier across a basal depression which was shown in the modelling study by Reference Vieli, Funk and BlatterVieli and others (2001). In order to obtain better predictions of tidewater glacier stability with the minimal model, we recommend taking into account glacier width as well as possible changes in the basal topography as the glacier terminus retreats. The minimal model reacts too quickly to the perturbed mass balance higher up the glacier (Fig. 7a) because in this model a decrease in volume, no matter where this is generated, implies that the glacier has to shorten immediately. In view of this, the minimal model cannot be realistic for very rapid changes in the forcing.
Sediment yield is an important control on terminus fluctuation of tidewater glaciers (Reference BoultonBoulton, 1970; Reference AlleyAlley, 1991; Reference Powell, Anderson and AshleyPowell, 1991). Push-moraine banks at the glacier terminus produce restraining forces, which may affect the ice flow (Reference Van der VeenVan der Veen, 1997; Reference Fischer and PowellFischer and Powell, 1998). During the glacier advance, accumulated sediment at the glacier terminus may decrease the relative water depth (Reference PowellPowell, 1981), reducing the calving flux and leading the glacier to advance into deeper water. For simplicity, however, we have not considered the formation of submarine moraine shoal in our models.
It should also be noted that we did not take up longitudinal stress gradients in our models. We assumed that the flow is controlled by a balance between driving stress and basal drag. Van der Reference Van der Veen and WhillansVeen and Whillans (1993) showed that for Columbia Glacier, a tidewater glacier, 80% of the flow resistance is due to basal drag (the rest is due to lateral drag and, to a lesser degree, gradients in longitudinal stress). Moreover, there is no evidence that gradients in longitudinal stress are important for the large-scale flow of a glacier (Reference Van der VeenVan der Veen, 1997).
The next step is to test the flowline models for real glaciers. Columbia Glacier is an obvious candidate because a unique dataset is available for this glacier (Reference FountainFountain, 1983; Reference KrimmelKrimmel, 1987), but the models should also be tested with other tidewater glaciers such as Hansbreen, Svalbard, and Breidamerkurjökull, Iceland.