Introduction
Melting glaciers, after ocean thermal expansion, are considered to be the second major contributor to the observed sealevel rise in the 20th century (Reference Church and HoughtonChurch and others, 2001; Reference Dyurgerov and MeierDyurgerov 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. Reference Gregory and OerlemansGregory and Oerlemans, 1998; Reference Van de Wal and WildVan de Wal and Wild, 2001; Reference Raper and BraithwaiteRaper 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 Reference Bahr, Meier and PeckhamBahr and others (1997).
Numerical iceflow 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. Reference Church and HoughtonChurch and others, 2001; Reference Van de Wal and WildVan de Wal and Wild, 2001; Radic ć and Hock, 2006) or for estimating volumes of existing glaciers (e.g. Reference Meier and BahrMeier and Bahr, 1996; Reference Raper and BraithwaiteRaper and Braithwaite, 2005). Volume and area for any glacier in a steady state are related via a power law, but under nonsteadystate conditions the powerlaw relationship may change as the massbalance profiles change (Reference Bahr, Meier and PeckhamBahr 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 iceflow model, Reference Van de Wal and WildVan de Wal and Wild (2001) assumed such differences to affect volume predictions of a retreating glacier by not more than 20%. Reference Pfeffer, Sassolas, Bahr and MeierPfeffer and others (1998) tested the powerlaw relation through the threedimensional iceflow modelling of synthetic glaciers in steady states, focusing on the estimation of glacier response times. Their results agreed well with the theory of Reference Bahr, Meier and PeckhamBahr and others (1997), but nonsteadystate conditions were not investigated.
In this study, we apply a onedimensional (1D) iceflow model to numerically generated synthetic glaciers in order to investigate the volume–area powerlaw relationships for both steadystate and nonsteadystate conditions. The main objectives are: (1) to determine and compare the relationships for steadystate and nonsteadystate conditions in order to test the validity of the powerlaw relationship for nonsteadystate conditions, and (2) to compare volume projections derived from volume–area scaling with those derived from the iceflow 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 volumescaling approach in comparison with the iceflow modelling are restricted to the 1D 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. Reference Macheret, Cherkasov and BobrovaMacheret and others, 1988; Reference Chen and OhmuraChen 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. Reference Bahr, Meier and PeckhamBahr and others (1997) have shown the physical basis for the powerlaw relationship by applying exponential relationships between various glacier characteristics such as length x, width w, slope α, shape factor F and massbalance profile b(x). These relationships are reasonable approximations to the behaviour of the ice flow and they include the flowlaw exponent n, the widthscaling exponent q in [w] ~ [x] ^{q} , the sidedrag scaling exponent f in [F] ~ [x] ^{–f} , the slope scaling exponent r in [sin α] ~ [x] ^{–r} and parameters determining the massbalance 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 (Reference Bahr, Meier and PeckhamBahr and others, 1997), the exact choice of characteristic values is not critical because each type of characteristic values is scaled by constants of proportionality. Massbalance profiles of valley glaciers can generally be expressed by:
where c _{1} and c _{0} are positive constants and m ≈ 2 (e.g. Reference DyurgerovDyurgerov, 1995; Reference Bahr, Meier and PeckhamBahr 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]:
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 (Reference Chen and OhmuraChen and Ohmura, 1990; Reference Meier and BahrMeier and Bahr, 1996).
Methods
Firstly, using a flowline model forced by different massbalance profiles, we produce 37 synthetic steadystate glaciers ranging in size from 4 to 58 km^{2}. Modelled volumes and areas are used to determine the scaling exponent γ in the volume–area relationship from regression analysis. Secondly, nonsteadystate conditions are modelled by imposing positive and negative massbalance 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 areaaveraged net massbalance computations and when the glacier is in nonsteadystate condition prior to the massbalance 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 1D iceflow model (central flowline along x) by Reference OerlemansOerlemans (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. Reference GreuellGreuell, 1989; Reference OerlemansOerlemans, 1997; Reference SchlosserSchlosser 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:
where b is massbalance rate, w the width of the glacier, h the surface elevation and D the diffusivity which is equal to:
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^{–3}m^{2} s^{–1} are taken from Reference Budd, Keage and BlundyBudd 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 Reference OerlemansOerlemans (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 timestep is sufficiently small (e.g. 0.005 years).
Set of synthetic steadystate 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 massbalance profile b(x) defined by different values of c _{1} and c _{0} (Equation (1)), in order to obtain a set of glaciers in steady states with different climate conditions and glacier sizes. We define the massbalance 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 massbalance 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 powerlaw 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 halfwidth much larger than glacier thickness (Reference NyeNye, 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 (Reference Bahr, Meier and PeckhamBahr 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 steadystate glacier’s length x is used to determine the glacier’s width w as:
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 steadystate glacier. We leave the value for scaling exponent r undefined a priori, as it is dependent on other scaling exponents and on the steadystate glacier geometry derived from the flowline model.
Modelderived volume–area relationships
From the 37 synthetic glaciers we obtain a set of values for V and A from which we determine the powerlaw relationship for steadystate conditions. To investigate volume–area scaling under nonsteadystate conditions (transient conditions) we introduce a perturbation δb in the massbalance profile on a subset of the 37 synthetic steadystate glaciers:
The magnitude of massbalance perturbation db(t) increases with a constant rate:
where t = 1, . . .100. A period of 100 years is chosen because future climatechange 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 powerlaw 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 steadystate 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 timestep which determines the volume change, in the scaling approach the volume changes are represented by
where is annual areaaveraged net mass balance. After the volume change has been calculated at t = 0, the glacier’s area at the next timestep 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 areaaveraged net mass balance is calculated from the massbalance profile as:
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 areaaveraged net mass balance, following Reference Elsberg, Harrison, Echelmeyer and KrimmelElsberg and others (2001) and Reference Harrison, Elsberg, Cox and MarchHarrison 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. 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. The area A is assumed constant in Equation (9) but variable in Equation (8), as done, for example, by Reference Radić and HockRadić 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. 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. areaaveraged mass balance of a valley glacier becomes less negative as the glacier, in its approach to a new steady state, retreats from lowlying highablation 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 steadystate glaciers, the values obtained from the volume evolutions of the nonsteadystate glaciers, and γ = 1.375 according to the theory of Reference Bahr, Meier and PeckhamBahr 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^{2} and [0.17, 10.29]km^{3}, respectively. The strong correlation shows that the flowline model produces glacier volumes and areas that follow a powerlaw relationship. The slope of the regression line corresponds to scaling exponent γ and it equals 1.56 with r ^{2} of 0.999. Hence, it differs from γ = 1.375 derived by Reference Bahr, Meier and PeckhamBahr and others (1997) by 14%. However, although theoretically derived, the value by Reference Bahr, Meier and PeckhamBahr 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.
Volume–area relationship in nonsteady state
The iceflow model produced volumes and areas as a discrete set of values with a timestep 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 massbalance perturbation of b(0) = –0.015ma^{–1}. We decided to treat the first 50 year period as a ‘discretization spinup’ period and consider the set of V and A in the remaining period as a representative set to derive scaling exponents under nonsteadystate 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 massbalance 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.
Volume evolutions: sensitivity experiments
Figure 5 illustrates the results for the total volume change over 100 years projected by the iceflow 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^{2}, V(0) = 5.77 km^{3}) 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 steadystate 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).
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 massbalance 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 massbalance 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^{2}) 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 massbalance 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’ massbalance scaling approach, method (c), but with three different scaling exponents as derived from our 37 steadystate synthetic glaciers, the transient evolutions of the synthetic glaciers and the theoretically derived value by Reference Bahr, Meier and PeckhamBahr 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 massbalance 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 nonsteady state, i.e. their mass balance has been negative or positive for several decades prior to the massbalance 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 massbalance perturbation compared to the results for glaciers initially in steady state (Fig. 5a). This is due to the larger magnitude of the initial massbalance 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 iceflow 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 nonsteady 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 massbalance perturbation according to Equation (7) is applied beyond the period of 100 years. How much they diverge depends on magnitude of massbalance perturbation, initial size of the synthetic glacier and the method for calculating areaaveraged 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 massbalance perturbation, as employed in our previous experiments (Equation (7)), we continued the evolution for an additional 200 years, keeping the massbalance 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 areaaveraged mass balance, thus excluding the feedback between the mass balance and glacier geometry changes. Only the method with ‘conventional’ massbalance 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 areaaveraged mass balance to geometry/elevation changes as simulated by the flowline model on a multicentury 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).
Conclusions
Scaling exponent γ = 1.56 in the volume–area relationship obtained from 37 synthetic steadystate glaciers of different sizes differed from γ = 1.375 derived theoretically by Reference Bahr, Meier and PeckhamBahr and others (1997) and from the exponents (γ = [1.80, 2.90]) derived for each of 24 investigated glaciers under nonsteadystate conditions, i.e. responding to hypothetical massbalance perturbations. Exponents γ were generally larger for negative massbalance 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 areaaveraged 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 massbalance computations. The most sophisticated method accounting for area changes and considering these in the massbalance 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 iceflow model when the glaciers are initially in nonsteady 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 areaaveraged massbalance 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 massbalance 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 massbalance–elevation feedback is captured by considering area changes in the massbalance 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 iceflow model. In a next step we will test the approach on real glaciers with observational records.
Acknowledgements
Financial support for this work is provided by the Swedish Research Council for Environment, Agricultural Sciences and Spatial Planning (project No. 21.4/20030387). R. Hock is Royal Swedish Academy of Science Research Fellow supported by a grant from the Knut Wallenberg Foundation. We are very grateful to S. Raper, R. van de Wal, M. Dyurgerov and two anonymous reviewers for useful comments and suggestions.