Volume–area scaling vs flowline modelling in glacier volume projections

Abstract

Volume–area scaling provides a practical alternative to ice-flow modelling to account for glacier size changes when modelling the future evolution of glaciers; however, uncertainties remain as to the validity of this approach under non-steady conditions. We address these uncertainties by deriving scaling exponents in the volume–area relationship from one-dimensional ice-flow modelling. We generate a set of 37 synthetic steady-state glaciers of different sizes, and then model their volume evolution due to climate warming and cooling as prescribed by negative and positive mass-balance perturbations, respectively, on a century timescale. The scaling exponent derived for the steady-state glaciers (γ = 1.56) differs from the exponents derived for the glaciers in transient (non-steady) state by up to 86%. Nevertheless, volume projections employing volume–area scaling are relatively insensitive to these differences in scaling exponents. Volume–area scaling agrees well with the results from ice-flow modelling. In addition, the scaling method is able to simulate the approach of a glacier to a new steady state, if mass-balance elevation feedback is approximated by removing or adding elevation bands at the lowest part of the glacier as the glacier retreats or advances. If area changes are approximated in the mass-balance computations in this way, our results indicate that volume–area scaling is a powerful tool for glacier volume projections on multi-century timescales.

Introduction

Melting glaciers, after ocean thermal expansion, are considered to be the second major contributor to the observed sea-level rise in the 20th century (Church and others, 2001; Dyurgerov and Meier, 2005). In future climate projections, the glacier contribution is expected to accelerate due to the fast response of glaciers to global warming, and much recent and ongoing research is focused on modelling and quantifying this future contribution (e.g. Gregory and Oerlemans, 1998; Van de Wal and Wild, 2001; Raper and Braithwaite, 2006). However, modelling future glacier contributions carries a variety of uncertainties. This is due to the scarcity of glacier inventory and hypsometry data and a large spectrum of uncertainties in modelling and downscaling future climate change, in modelling mass balance and finally in assessing the glacier volume changes. The uncertainties in modelling volume changes are addressed in this paper focusing on the volume–area scaling approach proposed by Bahr and others (1997).

Numerical ice-flow models best account for the physical processes involved, but they require detailed input data on glacier surface and bed geometry and therefore can only be applied on a small number of glaciers. Hence, owing to simplicity, the volume–area scaling approach has been widely used for considering area changes in volume predictions (e.g. Church and others, 2001; Van de Wal and Wild, 2001; Radic ć and Hock, 2006) or for estimating volumes of existing glaciers (e.g. Meier and Bahr, 1996; Raper and Braithwaite, 2005). Volume and area for any glacier in a steady state are related via a power law, but under non-steady-state conditions the power-law relationship may change as the mass-balance profiles change (Bahr and others, 1997), posing a problem for the employment of volume–area scaling in modelling the response of glaciers to future climate warming. Based on experiments with a numerical ice-flow model, Van de Wal and Wild (2001) assumed such differences to affect volume predictions of a retreating glacier by not more than 20%. Pfeffer and others (1998) tested the power-law relation through the three-dimensional ice-flow modelling of synthetic glaciers in steady states, focusing on the estimation of glacier response times. Their results agreed well with the theory of Bahr and others (1997), but non-steady-state conditions were not investigated.

In this study, we apply a one-dimensional (1-D) ice-flow model to numerically generated synthetic glaciers in order to investigate the volume–area power-law relationships for both steady-state and non-steady-state conditions. The main objectives are: (1) to determine and compare the relationships for steady-state and non-steady-state conditions in order to test the validity of the power-law relationship for non-steady-state conditions, and (2) to compare volume projections derived from volume–area scaling with those derived from the ice-flow modelling.

Using synthetic glaciers has the advantage that it easily allows modelling of a large number of glaciers under defined and controlled conditions, but it must be borne in mind that conclusions on the validity of the volume-scaling approach in comparison with the ice-flow modelling are restricted to the 1-D flowline representation of glaciers as defined in our experiments. In a next step we will elaborate on this paper by considering real glaciers constrained by observational data.

Theory of Volume–Area Scaling

Several authors (e.g. Macheret and others, 1988; Chen and Ohmura, 1990) have suggested that the volume V of valley glaciers is proportional to the surface area A raised to a power of about 1.36. Bahr and others (1997) have shown the physical basis for the power-law relationship by applying exponential relationships between various glacier characteristics such as length x, width w, slope α, shape factor F and mass-balance profile b(x). These relationships are reasonable approximations to the behaviour of the ice flow and they include the flow-law exponent n, the width-scaling exponent q in [w] ~ [x] ^q , the side-drag scaling exponent f in [F] ~ [x] ^–f , the slope scaling exponent r in [sin α] ~ [x] ^–r and parameters determining the mass-balance profile. The brackets indicate characteristic values which can be assumed as the glacier’s mean width, maximum length, total area, etc. In geometric scaling analysis (Bahr and others, 1997), the exact choice of characteristic values is not critical because each type of characteristic values is scaled by constants of proportionality. Mass-balance profiles of valley glaciers can generally be expressed by:

(1)

where c ₁ and c ₀ are positive constants and m ≈ 2 (e.g. Dyurgerov, 1995; Bahr and others, 1997). This quadratic balance profile is the dominant term in a polynomial expansion of the actual mass balance. Dimensional analysis of glacier characteristics results in the following relation between [V] and [A]:

(2)

where for valley glaciers q = 0.6, m = 2 and f = 0 are suggested by empirical data, and r is either 0 for steep surface slopes or r = (1–m+n–nf)/2(n+1) for shallow slopes. Inserting these values into Equation (2), the exponent in volume–area relationship equals γ = 1.375 which is in close agreement with the exponent γ = 1.36 which has been empirically derived from many Eurasian and Alpine glaciers (Chen and Ohmura, 1990; Meier and Bahr, 1996).

Methods

Firstly, using a flowline model forced by different mass-balance profiles, we produce 37 synthetic steady-state glaciers ranging in size from 4 to 58 km². Modelled volumes and areas are used to determine the scaling exponent γ in the volume–area relationship from regression analysis. Secondly, non-steady-state conditions are modelled by imposing positive and negative mass-balance perturbations on a subset of these synthetic glaciers producing in total 24 volume evolutions for 100 years. For each volume evolution, we derive the scaling exponent γ based on the annual transient values of volume and area. Thirdly, we use the volume–area scaling approach to model glacier volume evolutions and compare results to those obtained by the flowline model. Finally, we apply several sensitivity experiments to evaluate the scaling approach when geometry changes are excluded or included in area-averaged net mass-balance computations and when the glacier is in non-steady-state condition prior to the mass-balance perturbations. We also investigate the sensitivity of results to the choice of the scaling exponent and the sensitivity of results in scenarios where the climate stabilizes after a period of perturbation.

The flowline model

We use the 1-D ice-flow model (central flowline along x) by Oerlemans (1997). We consider this model as a good reference for evaluating the scaling approach since the model has proved to perform well in reconstructions of real glacier fluctuations (e.g. Greuell, 1989; Oerlemans, 1997; Schlosser 1997). The model equations are generated from the vertically integrated continuity equation, assuming incompressibility of ice, and Euler’s equations combined with Glen’s flow. From these equations the prognostic equation for thickness H is derived as:

(3)

where b is mass-balance rate, w the width of the glacier, h the surface elevation and D the diffusivity which is equal to:

(4)

where ρ is ice density and g is acceleration of gravity. Values for deformation parameter f _d = 1.9×10^–24 Pa^–3 s^–1 and sliding parameter f _s = 5.7×10^–20 Pa^–3m² s^–1 are taken from Budd and others (1979). This assumes that the vertical mean ice velocity is entirely determined by the local ‘driving stress’ τ and it has two components: one associated with internal deformation f ^d Hτ and the other with basal sliding f ^s τ/H. The ‘driving stress’ τ is proportional to the ice thickness H and surface slope ∂h/∂x. For further details about the model, the reader is referred to Oerlemans (1997). Equation (1) is solved on a 100m resolution along the flowline while time integration is done with a forward explicit scheme which is stable if the time-step is sufficiently small (e.g. 0.005 years).

Set of synthetic steady-state glaciers

We apply the flowline model to generate a set of synthetic glaciers defined as slabs of ice with uniform widths lying on a bed with uniform slope (tan α = 0:1). The model is run for 37 choices of mass-balance profile b(x) defined by different values of c ₁ and c ₀ (Equation (1)), in order to obtain a set of glaciers in steady states with different climate conditions and glacier sizes. We define the mass-balance profile as a function of elevation b(h) which is then transformed into a function of horizontal position b(x) by fitting a parabolic function. By doing this we estimate the value for parameter c_m with scaling exponent m = 2. An example is shown in Figure 1 where the mass-balance profile b(h) is approximated by the profile b(x). We consider the glacier in steady state if modelled glacier volume and area remain unchanged over a period of 100 years. In order to get the scaling exponent γ to agree with the theory of power-law relation for valley glaciers, we have chosen the following set of exponents: q = 0.6, f = 0, m = 2 and n = 3. In the theory the choice of f = 0 corresponds to the glaciers with little side drag, i.e. the glaciers with half-width much larger than glacier thickness (Nye, 1965). Although this may not be a good approximation for valley glaciers, the empirical data from valley glaciers showed that f is expected to be close to zero (Bahr and others, 1997). Therefore, to achieve f ≈ 0, we produced synthetic glaciers with large widths relative to their thickness. To obtain q = 0.6 we run the flowline model with a priori determined uniform width and produced a glacier in a steady state. Then the steady-state glacier’s length x is used to determine the glacier’s width w as:

Fig. 1. Mass-balance profile b(h) = –1.45×10^–6 h ² + 0.0085h –9.5, where h is elevation (m a.s.l.), and its approximation with the profile b(x) = –3.564×10^–8 x ² + 1.614, where x is length along the flowline (m). The profile generated the synthetic glacier with area A = 31.64 km² and V = 4.08 km³.

(5)

where c = 10m ^1–q is our choice for constant of proportionality. The flowline model is then rerun using the derived width to produce the synthetic steady-state glacier. We leave the value for scaling exponent r undefined a priori, as it is dependent on other scaling exponents and on the steady-state glacier geometry derived from the flowline model.

Model-derived volume–area relationships

From the 37 synthetic glaciers we obtain a set of values for V and A from which we determine the power-law relationship for steady-state conditions. To investigate volume–area scaling under non-steady-state conditions (transient conditions) we introduce a perturbation δb in the mass-balance profile on a subset of the 37 synthetic steady-state glaciers:

(6)

The magnitude of mass-balance perturbation db(t) increases with a constant rate:

(7)

where t = 1, . . .100. A period of 100 years is chosen because future climate-change studies are often focused on a century scale. We chose three different initial magnitudes for db(0), corresponding to climate cooling and warming: 0.005, 0.01 and 0.015 ma^–1. In total we create 24 volume evolutions of glaciers with different initial sizes (12 responding to climate warming and 12 responding to cooling). We determine a power-law relationship for each of these 24 transient volume and area evolutions by linearly regressing on logarithmic axes the modelled annual values of volume and area once the steady-state area has changed.

Volume projections using volume–area scaling

Finally, we investigate how well volume evolutions can be estimated from the volume–area scaling approach by comparing results to those obtained from the flowline model. For each of the 24 volume evolutions produced by the flowline model, we compute corresponding volume evolutions based on the scaling approach. While the flowline model calculates the thickness change for each time-step which determines the volume change, in the scaling approach the volume changes are represented by

(8)

where is annual area-averaged net mass balance. After the volume change has been calculated at t = 0, the glacier’s area at the next time-step t = 1 is calculated by inverting Equation (2). The new area at t = 1 is used to calculate the mass balance (Equation (9)) and the volume change at t = 1 by again using Equation (8), and the calculation is repeated until t = 100. Annual area-averaged net mass balance is calculated from the mass-balance profile as:

(9)

where bi and a_i are discrete values of mass balance b(x,t) and surface area a(x,t ) for each spatial band (i = 1. . .n) over the glacier length, while A is total surface area. We use two definitions for annual area-averaged net mass balance, following Elsberg and others (2001) and Harrison and others (2005): If is calculated keeping surface area constant in time (equal to initial area A(t = 0)) the result is a ‘reference surface’ mass balance. If area in Equation (9) is updated for each year by volume–area scaling, we derive ‘conventional’ mass balance. Here we assume that change in total area occurs on the tongue of the glacier, so the lowest area bands are excluded if total area shrinks, or new area bands are included if total area grows. Area bands are 100 m long to be equal in size to the grid spacing in the flowline model.

For comparison, we apply three different methods to calculate volume evolution from Equations (8) and (9) differing from each other solely in whether or not area changes are included in Equations (8) and (9):

1. 1. The glacier area A is assumed constant in both Equations (8) and (9). Hence volume–area scaling is not applied. is calculated as a ‘reference surface’ mass balance using constant area A and a constant number of spatial bands in Equation (9).

2. 2. The area A is assumed constant in Equation (9) but variable in Equation (8), as done, for example, by Radić and Hock (2006). The glacier area is adjusted according to volume–area scaling, i.e. a new area is computed using Equation (2) from the volume change obtained by Equation (8), but is computed as a ‘reference surface’ mass balance using constant area (Equation (9)).

3. 3. Area changes are considered in both Equations (8) and (9). is calculated as a ‘conventional’ mass balance, meaning that A(t) and the number of spatial bands changes in Equation (9). Volume–area scaling is applied. This method partially accounts for the feedback between geometry changes and mass balance (e.g. area-averaged mass balance of a valley glacier becomes less negative as the glacier, in its approach to a new steady state, retreats from low-lying high-ablation altitudes to higher and colder altitudes).

We compare the results of each method with the volume evolution derived from the flowline model. To test how sensitive the volume projections are to the choice of scaling exponent in the scaling methods, we use the exponent γ in volume–area relation derived from the 37 steady-state glaciers, the values obtained from the volume evolutions of the non-steady-state glaciers, and γ = 1.375 according to the theory of Bahr and others (1997).

Results and Discussion

Volume–area relationship in steady state

Figure 2 illustrates the relationship between volume and area plotted on logarithmic axes for all 37 synthetic glaciers in steady states. Glacier areas and volumes span in the range [4.37, 57.88]km² and [0.17, 10.29]km³, respectively. The strong correlation shows that the flowline model produces glacier volumes and areas that follow a power-law relationship. The slope of the regression line corresponds to scaling exponent γ and it equals 1.56 with r ² of 0.999. Hence, it differs from γ = 1.375 derived by Bahr and others (1997) by 14%. However, although theoretically derived, the value by Bahr and others (1997) is largely dependent on empirical data to which their results were adjusted. Since we analyze synthetic glaciers the deviation from empirical results was expected because our synthetic glaciers have largely simplified geometry (e.g. uniform widths) and are created with the flowline model which presents a simplified approximation for glacier dynamics. Below, we aim to answer how significant this deviation is when deriving volume evolutions.

Fig. 2. Log–log plot of volume V vs area A for 37 synthetic glaciers in steady state with a regression log V = 1.56 logA – 2.11 (r ² = 0.999).

Volume–area relationship in non-steady state

The ice-flow model produced volumes and areas as a discrete set of values with a time-step of 1 year. While volume changes occur almost immediately after introducing the perturbation δb, due to discretization, modelled length and surface area remain constant for an initial period of ~30–50 years, depending on the magnitude of perturbation. As an example, Figure 3 shows the area evolution in response to the mass-balance perturbation of b(0) = –0.015ma^–1. We decided to treat the first 50 year period as a ‘discretization spin-up’ period and consider the set of V and A in the remaining period as a representative set to derive scaling exponents under non-steady-state conditions. Figure 4 illustrates all 24 sets of V and A on logarithmic axes. Scaling exponents are derived for each of the 24 volume evolutions and they span in the range [1.80, 2.90] for γ, with a corresponding range of [–3.88, –12.01] for k in log V = γ logA + k. Larger values for γ tended to occur for negative mass-balance perturbations (warming scenario) compared to positive perturbations (cooling scenario), and γ tended to decrease with increasing initial glacier size. Scaling exponents for our set of test glaciers differ by 21% (γ = 1.80) to 86% (γ = 2.90) from the scaling exponent derived for the synthetic glaciers in steady states (γ = 1.56). One possible reason for this difference is that the glacier’s width in the transient state is not scaled with the glacier’s length according to Equation (5), meaning that the scaling parameter q in the width–length relationship may change in time since the glacier’s length changes while the width is kept constant. Thus, changes in any of the exponential relationships between glacier characteristics, such as between width and length, modify the scaling exponent in the volume–area relationship. The significance of this difference in scaling exponent γ is analyzed in the next section.

Fig. 3. Surface area evolution derived from the flowline model as a response to the perturbation of db(0) = –0.015ma^–1.

Fig. 4. Log–log plot of volume vs area for 24 volume (area) evolutions as a response to different mass-balance perturbations. Each evolution contains 50 values of V and A. For each evolution we derived a regression log V = γ log A+ k in order to determine the scaling exponent γ.

Volume evolutions: sensitivity experiments

Figure 5 illustrates the results for the total volume change over 100 years projected by the ice-flow model compared to those projected by the scaling approach. Here, we illustrate only 2 of the 24 evolutions since results in terms of sensitivity to the choice of method are similar for all evolutions. We choose the largest glacier in the set (A(0) = 38.92 km², V(0) = 5.77 km³) responding to the largest perturbation of db(0) = 0.015ma^–1. The evolutions are normalized to the initial volumes. In Figure 5a we compare three different variants of the scaling approach as described above ((a) the ‘reference surface’ mass balance with no volume–area scaling, (b) the ‘reference surface’ mass balance with scaling and (c) the ‘conventional’ mass balance with scaling). The scaling exponent γ = 1.56 as previously derived from the 37 synthetic steady-state glaciers is applied. To optimize initial conditions for each volume evolution the constant of proportionality in the volume–area scaling relationship is derived from the glaciers’ initial volume and area instead of using the constants k derived from regression analysis of each evolution (Fig. 4).

Fig. 5. Normalized volume evolutions of the largest test glacier responding to the mass-balance perturbation of db(0) = + 0.015 ma^–1 (‘cooling scenario’) and db(0) = –0.015ma^–1 (‘warming scenario’). In (a) and (c) the three methods correspond to three different ways of calculating area-averaged net mass-balance and volume changes: the ‘reference surface’ mass balance without volume–area scaling (method a), the ‘reference surface’ mass balance with scaling (method b) and the ‘conventional’ mass balance with scaling (method c). Scaling exponent γ = 1.56 is used in the volume–area relationship. In (b) volume evolutions are derived from scaling method (c) using three different scaling exponents: γ = 1.56 derived from our 37 synthetic steady-state glaciers, γ = 2.27 (1.89) derived from the transient response to warming (cooling) of this test glacier, and γ = 1.375 derived theoretically by Bahr and others (1997). In (c) the test glacier is in non-steady state prior to the mass-balance perturbation.

Results based on methods (a) and (b) closely follow the evolution curves from the flowline model in the first ~50 years for both the warming and the cooling scenario, while those from method (c) notably deviate somewhat earlier (Fig. 5a). However, by the end of the 100 year period the volume change obtained from method (c) is smallest. Note that the volume response is not symmetrical for positive and negative mass-balance perturbation of equal magnitude for methods (b) and (c) which include scaling. This is due to the exponential character of the volume–area relationship (Equation (2)). Results from the whole set of 24 evolutions showed that by increasing the magnitude of the mass-balance perturbation, the sensitivity to the choice of the scaling method increases, i.e. the difference between the projections derived by the flowline model and the scaling methods increases. Also, smaller glaciers in the set (A<20km²) are more sensitive to the choice of scaling method. However, these differences in total volume change over 100 years for the whole set of 24 evolutions are within the range of 12% of initial volume when method (c) is applied and 16% and 23% when methods (a) and (b) are applied, respectively. Thus the smallest differences are produced by method (c). This was the expected result since scaling method (c) is the most sophisticated of the three, taking into account area changes and considering these in the mass-balance computations.

The next sensitivity test quantifies the uncertainty in volume projections due to different values of scaling exponent γ in volume–area scaling. For this purpose we use the ‘conventional’ mass-balance scaling approach, method (c), but with three different scaling exponents as derived from our 37 steady-state synthetic glaciers, the transient evolutions of the synthetic glaciers and the theoretically derived value by Bahr and others (1997). Results for the largest glacier in the set are shown in Figure 5b. In the total set of 24 evolutions the 100 year volume changes derived by the scaling method with three different scaling exponents differ from each other less than 6%. The difference tends to decrease with decreasing mass-balance perturbations or increasing initial glacier size. These results suggest that applying scaling exponents that differ up to 86% yields differences not larger than 6% in derived volume changes on a century timescale. This difference may be considered negligible in comparison to the range of differences due to the choice of the scaling methods and the range of uncertainties due to the approximations in the flowline model and volume–area scaling approach. Additionally, applying the scaling exponent γ derived from each transient evolution of the synthetic glacier produced the volume projections that followed most closely those from the flowline model. This is to be expected since the scaling exponent is calculated directly from the relationship between transient volumes and areas produced by the flowline model.

So far we have evaluated the scaling approach for synthetic glaciers that are initially in steady states. In the next sensitivity experiment we compute the volume evolutions for glaciers that are initially in non-steady state, i.e. their mass balance has been negative or positive for several decades prior to the mass-balance perturbation. The results for the largest synthetic glacier in the set are shown in Figure 5c. Initial mass balance for the glacier with warming scenario is b(0) = –0.54ma^–1, and it is perturbed with db(0) = –0.015ma^–1, while for the cooling scenario the values are b(0) = 0.55 ma^–1 and db(0) = 0.015 ma^–1. All scaling methods show a stronger response to the mass-balance perturbation compared to the results for glaciers initially in steady state (Fig. 5a). This is due to the larger magnitude of the initial mass-balance perturbation. In addition, deviations between the different projections are much larger. For all 24 evolutions, the differences between the 100 year volume changes obtained from the ice-flow model and the volume changes from methods (a) and (b) are up to 47% and 74%, respectively, while the volume changes from method (c) differ up to 16% of initial volume. Thus, method (c) produced the best approximation of 100 year volume evolutions derived from the flowline model for the synthetic glaciers initially in non-steady state. We also assumed different scaling exponents in scaling method (c), as done above, and derived the 100 year volume changes which differed by <12%. As in the experiment above, applying the scaling exponent γ derived from each transient evolution produced the closest volume projection to that obtained from the flowline model. We expect the volume projections derived from the scaling approach to continue to diverge from those derived by the flowline model if the mass-balance perturbation according to Equation (7) is applied beyond the period of 100 years. How much they diverge depends on magnitude of mass-balance perturbation, initial size of the synthetic glacier and the method for calculating area-averaged mass balance.

Our final sensitivity test evaluates the scaling methods for hypothetical scenarios where the climate stabilizes after the initial period of perturbation. To that end, we derived volume projections for a 300 year period applying a cooling or warming scenario for the first 100 years while afterwards the climate is kept stable. Thus, after an initial period of 100 years with mass-balance perturbation, as employed in our previous experiments (Equation (7)), we continued the evolution for an additional 200 years, keeping the mass-balance perturbation equal to the perturbation at t = 100. The results are shown in Figure 6 for a cooling and warming scenario applied on the largest glacier in the set. For both scenarios the volume evolutions derived from the flowline model reach new steady states. This is not the case for scaling methods (a) and (b) which keep the surface area constant in the calculations of area-averaged mass balance, thus excluding the feedback between the mass balance and glacier geometry changes. Only the method with ‘conventional’ mass-balance calculation, method (c), is able to simulate the approach of the glacier to a new steady state. Although method (c) produces 100 year volume changes which deviate up to 12% from the changes derived by the flowline model, it is the only one of the three scaling methods that is capable of simulating the response of area-averaged mass balance to geometry/elevation changes as simulated by the flowline model on a multi-century timescale. For our synthetic glacier with uniform width, this feedback is simulated by subtracting (adding) area bands on the tongue of the glacier as the glacier retreats (grows) due to warming (cooling).

Fig. 6. Evolution of area-averaged mass balance b (a, c) and normalized glacier volume V (b, d) derived from the flowline model and from the scaling methods. Initial perturbation is db(0) = 0.015 ma^–1 (cooling scenario) and db(0) = –0.015ma^–1 (warming scenario). Scaling methods (a) and (b) are based on ‘reference surface’ mass balances, and scaling method (c) is based on ‘conventional’ mass balances as also in the flowline model.

Conclusions

Scaling exponent γ = 1.56 in the volume–area relationship obtained from 37 synthetic steady-state glaciers of different sizes differed from γ = 1.375 derived theoretically by Bahr and others (1997) and from the exponents (γ = [1.80, 2.90]) derived for each of 24 investigated glaciers under non-steady-state conditions, i.e. responding to hypothetical mass-balance perturbations. Exponents γ were generally larger for negative mass-balance perturbations (warming scenarios) than for positive perturbations (cooling scenarios), and γ tended to decrease with increasing initial glacier size. However, the range of differences in scaling exponent by up to 86% is shown to make negligible differences, <6%, in 100 year volume changes derived from the scaling approach.

Volume projections on a century timescale differed within the range of 12–23% of initial volume from the flow model results, depending on the method by which the area-averaged net mass balance is calculated, i.e. whether or not volume–area scaling is applied and area changes obtained from volume–area scaling are included (‘conventional’ mass balance) or excluded (‘reference surface mass balance’) in the mass-balance computations. The most sophisticated method accounting for area changes and considering these in the mass-balance computations resulted in the smallest differences (up to 12%) in projected volume changes over 100 years. This method best agreed with the projections by the ice-flow model when the glaciers are initially in non-steady state or when the climate is assumed to stabilize after a period of perturbation. In fact, the method is capable of simulating the glacier approaching a new steady state by simulating the feedback between area-averaged mass-balance and glacier geometry/elevation changes resulting from retreat or advance of the glacier. This feedback is captured by excluding area from or adding area to the lowest part of the glacier. In contrast, neglect of volume– area scaling and neglect of area changes in the mass-balance computations fails to simulate this feedback and the approach to a new steady state.

Although based on a set of synthetic glaciers of highly simplified geometry, our results are promising for use of volume–area scaling in glacier volume projections provided that the mass-balance–elevation feedback is captured by considering area changes in the mass-balance computations. Our approach to add area to and remove area from the lowest elevation bands of the glacier seems to be able to capture these processes sufficiently well to obtain results comparable to those from the ice-flow model. In a next step we will test the approach on real glaciers with observational records.