Sensitivity studies on model modifications to assess the dynamics of a temperate ice cap, such as that on King George Island, Antarctica

Abstract

Numerical model studies describing the dynamics of ice sheets are usually aimed at cold- ice regions of mainland Antarctica. Contrary to these ice bodies with temperatures well below the pressure-melting point at the ice surface, the ice caps on sub-Antarctic islands, such as King George Island, can be considered as polythermal or even temperate. This implies that they may contain nonnegligible amounts of water percolating through the ice matrix. We present three different modifications to a three-dimensional, thermomechanically coupled, higher-order model previously applied to cold-ice regions. We discuss the effect of these modifications on the ice dynamics obtained for diagnostic model runs. The modifications comprise changes to the enhancement factor in Glen’s flow law, the choice of negligible or non-negligible water content and the choice of a threshold altitude below which the ice surface is set to pressure-melting point conditions. All modifications lead to non-linear changes in the resulting horizontal flow velocities. The changes in velocity amplitudes obtained with these modifications compared to simulations without any modification range between 64% and ~400%. This implies that they should be considered in time-dependent simulations of temperate-ice dynamics.

Introduction

King George Island (KGI), the largest of the South Shetland Islands, is located north of the Antarctic Peninsula at 62° S (Fig. 1). It is largely governed by maritime climate conditions. Several in situ measurements have shown that the mean annual air temperature lies just below melting conditions (e.g. –2.4°C; Wen and others, 1994) with a summer maximum air temperature well above 0°C at sea level (0.9°C; Ferron and others, 2004). A possible warming (even only by 1 °C) would therefore have significant impacts on the dynamics and the mass balance of the ice cap which covers more than 92% of the island (~ 1250 km²; Braun and Hock, 2004). As in the Antarctic Peninsula, changes in the climatic and glacial systems are already being observed. This fact poses a challenging task for numerical models that are used to predict possible future behaviour of ice bodies in response to climate change. Though the ice cap on KGI is rather small compared to the large ice sheets on Greenland and Antarctica, it is advisable to simulate the ice dynamics using numerical flow models. For example, Knap and others (1996) presented a climate sensitivity study applying a twodimensional (2-D) vertically integrated ice-flow model. They found that the ice cap reacts more sensitively to changes in temperature than to changes in precipitation. Barboza and others (2004) recently published results of a 2-D numerical finite-difference model applied to a glacier on KGI. However, in their study, basal sliding processes are neglected.

Fig. 1. Satellite image of KGI. The tracks of the ice-thickness measurements are indicated (Pfender, 1999). In the lower right corner, the location of KGI is shown (personal communication from M. Braun, 2003).

The largest problem that is not considered in either model is that one has to distinguish cold and usually extended ice sheets with temperatures below the pressure-melting point from smaller and temperate ice caps with temperatures at the pressure-melting point. This is because in the case of cold ice only a negligible amount of water is found in the ice body. In warmer ice, water often moves through the ice matrix and affects the ice-flow behaviour in a non-negligible way. Hutter (1993) even suggests applying a Darcy-type permeability and treating a temperate glacier as a two-phase flow problem. But, since this approach is difficult to integrate into numerical flow models, it has not yet been followed. Here we attempt to find a more feasible way to treat a temperate ice cap numerically.

Regarding the numerical flow model previously applied to cold-ice regions in Antarctica (e.g. Paschke and Lange, 2003), there are several aspects to be considered if one wants to modify the model for a simulation of the ice dynamics of a temperate ice cap. First, a temperate ice column does not include a firn layer on top. Due to refreezing of percolating surface- and meltwaters, the ice column can be considered as almost entirely consolidated. Second, within the annual cycle, the surface temperature may reach melting conditions. The altitude up to which the surface temperature is set to 0°C has to be chosen. Third, the water content should be considered. In the present study, we modified the three-dimensional (3-D) numerical flow model developed by Sandhäger (2000), applying all these considerations. Furthermore, we tried to assess the dependency of the so-called enhancement factor in Glen’s flow law on various parameters.

These modifications are used for the diagnostic simulation of ice dynamics of an ice body with a geometry similar to that of the KGI ice cap. We present model results that reveal the sensitivity of the numerical flow model to modifications of the parameters just introduced. In this context, the importance of such parameter modifications for a time-dependent model run are discussed.

Methods

The 3-D numerical flow model

The original 3-D inland-ice/ice-shelf flow model was developed by Sandhäger (2000) and it has been applied to several cold-ice regions on mainland Antarctica (e.g. Sandhäger, 2003). For the present study, we only utilize the grounded-ice part of the numerical model. As already described in Paschke and Lange (2003), the 3-D stress and velocity fields are computed taking the role of normal deviatoric stress gradients into account. Sandhäger (2000) used the numerical integration scheme proposed by Blatter (1995) to solve the exact, unscaled equations. In so doing, the ice-flow velocities are calculated depending on the simulated temperature distribution.

One important boundary condition of a numerical model applied to temperate or polythermal glaciers is the treatment of basal velocities. We applied the same Weertman-type sliding law as Sandhäger (2000) and Paschke and Lange (2003). Detailed testing of different sliding laws is being addressed in ongoing and future studies.

Adaptation of the 3-D model to a temperate ice cap

Ice geometry

For the present sensitivity study, we investigate an ice body whose geometry (i.e. surface elevation, ice thickness, ice front position) is similar to that of the KGI ice cap. The surface elevation is based on a digital elevation model (DEM) published by Braun and others (2001; Fig. 2). They obtained elevations on a 50 m × 50 m grid by combining information from different large-scale maps. The observational density decreases from west to east such that the most eastern part of the island is covered only by one large-scale map. But our focus lies on the western part of KGI, so the information on surface elevations needed on a 250 m × 250 m grid for the numerical simulation is satisfactory.

Fig. 2. Surface elevation for the study region published by Braun and others (2001). The study region for the numerical simulation is indicated, as are the locations of two flowlines A and B. Slices through the ice body along these flowlines are presented in Figures 9–11.

Restrictions on the western part of the island are mainly due to the limited knowledge of the ice-thickness distribution. Except for an area of 200 km² in the western part of the ice cap where a field campaign took place during the austral summer 1997/98, the information is rather poor. N. Blindow and M. Pfender of the Institute for Geophysics in Münster, Germany, investigated the bedrock topography using radio-echo sounding with a high spatial resolution along the tracks described by Blindow (1994; Fig. 1). Due to the crevassed ice surface, the tracks could not be extended to the ice front. Therefore, the ice-thickness distribution had to be manually extrapolated towards the ice front, keeping in mind characteristic features observed in the course of the in situ measurements. The ice-thickness distribution obtained should be considered as a preliminary result of the field measurements. The final results of the ice-thickness distribution on western KGI will be published elsewhere by Blindow and others after additional postprocessing of the data. It is more important for a sensitivity study to apply the same ice geometry to all model runs than it is to use the ice-thickness distribution of a real ice cap. We therefore applied a smoothed version of the preliminary ice-thickness distribution. For this reason, we speak of an ice geometry closely resembling KGI’s ice cap and not of KGI’s ice geometry.

The location of the ice front was determined using satellite images of KGI. The northern and southern boundaries of the study region were chosen following flowlines.

Ice-density distribution

Due to the liquid water content and the refreezing surface-and meltwaters, the KGI ice cap has a firn layer different from that of cold-ice regions on mainland Antarctica. The vertical distribution of the ice density also follows a significantly different pattern than on mainland Antarctica. Simões and others (2004) published borehole measurements to 49.9m depth drilled at 690ma.s.l. They found that the density of pure ice is reached at 35 m depth. This means that less than the upper 10% of the ice column consists of unconsolidated ice. Wen and others (1998) observed the transition depth from firn to ice at a sampling site about 250 m a.s.l. to be 7 m below the surface. This strengthens the hypothesis that the transition depth from ice with lower density to ice with a pure-ice density decreases with lowering altitudes.

For the present study, we applied a constant value of pure- ice density for the lower part of the ice column and a linear decrease above a transition depth. The position of this transition depth within each vertical ice column varies linearly between 35 m at the top of the ice cap and 0m below 200 m a.s.l. Regarding the vertical resolution of the numerical simulation, one could even justify the application of pure-ice density for all gridpoints since only the uppermost one or two gridpoints of each ice column (15 gridpoints) will be affected by the aforementioned relation.

Temperature boundary conditions

One of the results of the numerical flow model is the temperature distribution within the ice body. Nevertheless, the temperature at the ice surface has to be prescribed as a boundary condition. For cold-ice regions, most authors use the ice temperature 10m below the ice surface where the seasonal effects fade out (Loewe, 1970). On KGI, this so- called active layer is considered to reach a little deeper, down to 15 m (Wen and others, 1998). Wen and others observed ice cores at different elevations between 100 and 700 m a.s.l. and found ice surface temperatures between –1.9 and 0°C that did not correlate with the surface elevation of the respective sampling points. They assumed that the thickness of the snow–firn layer, as well as the percolation of surface water, influences the ice temperature, which may lead to increasing ice temperatures with altitude at the bottom of the active layer, but they do not verify this hypothesis.

For the KGI ice cap, only the aforementioned few borehole measurements of ice surface temperatures are available. Furthermore, only four of these points are located at the boundary of our study region, all showing ice surface temperatures between –0.2 and 0°C. All other measurements were carried out outside our study region and also at lower altitudes. Thus, we consider the air temperatures representing a linear lapse rate (e.g. Braun and others, 2004) more appropriate and apply a linear gradient for decreasing ice surface temperatures with increasing altitude of 0.0073°Cm–¹. In order to evaluate the dependence of ice dynamics on ice surface temperature, we vary the threshold altitude up to which ice surface temperatures are being set to 0°C from 0 m, 100 m, 200 m, etc. up to the full height of the ice cap. Above this threshold altitude, we apply a linear ice surface temperature gradient depending on the altitude as mentioned above. The case with 0°C ice surface temperatures over the entire ice cap nearly represents the ice temperatures measured by Wen and others (1998). Thus, the additional model runs serve to clarify the dependence of the ice dynamics on this boundary condition.

At the ice base, we apply a geothermal heat flux of 90mWm^–2. This value was chosen as it lies within the interval published by Pollack and others (1993) for Mesozoic to Cenozoic igneous rocks. Following Birkenmajer (1997), the study region on KGI mainly consists of Upper Cretaceous to Oligocene terrestrial, basaltic and andesitic rocks.

Modifications to Glen’s flow law

One method that is often applied by glaciologists to describe the rheological behaviour of an ice body is represented by Glen’s flow law. Using the following equation (e.g. Paterson, 1994), the relationship between the deformation rate (with i,j = 1,2,3) and the deviatoric stress tensor τij is given by

(1)

where (u ₁, u ₂, u ₃) represents the ice velocity field, τ the effective deviatoric stress, T ^* the temperature corrected for dependencies of the melting point on pressure (e.g. Huybrechts, 1992) and w the water content in per cent; the exponent in Glen’s flow law is typically chosen to be n = 3 for grounded ice bodies.

The two other parameters appearing in this relation are the Arrhenius factor A and the enhancement factor m. Both parameters have to be selected separately depending on the specific glaciological conditions in the study region.

The Arrhenius factor A

The Arrhenius factor depends on different conditions in the ice body, mainly on the temperature T*, but if the ice body reaches temperatures near the pressure-melting point, the water content w has to be taken into account. Greve and others (1998) attempt to distinguish different temperature regimes. Approaching T* = 0°C, the transition between a purely temperature-dependent Arrhenius factor and an Arrhenius factor that depends only on the water content has to be specified. However, the formulations given by Greve and others (1998) only show a continuous, smooth transition if the water content w = 0%. We propose to modify these relations by taking into account the influence of the water content w on the Arrhenius factor from –2°C with an increasing effect until the pressure-melting point is reached. Contrary to the water content w, the influence of the potential temperature T* on the Arrhenius factor decreases between –2 and 0°C. Similarly, the following relations are the same as for the purely temperature- dependent one proposed by Greve and others (1998) if the water content is 0%. And in the case of pressure-melting point conditions (i.e. T* = 0°C), our equations become the same as that proposed by Greve and others (1998), depending only on the water content. Since the ice cap on KGI never reaches temperatures below –10°C, the following relations are given for the interval [–10;0]°C ([A] = 10^–16Pa^–3a^–1):

(2)

Direct quantitative measurements of the water content w in the KGI ice cap are presently not available, to our knowledge. However, ground-penetrating radar (GPR) investigations on the KGI ice cap indicate finite amounts of water, which cannot be neglected if reasonable modelling results are desired. To address this problem, we chose to compare our work with investigations on a temperate glacier, Valle´e Blanche, French Alps, by Vallon and others (1976). They report data indicating a continuous fit between water content and normalized depth to the fourth power. Applying this relation to those parts on the KGI ice cap below 400 m a.s.l., and considering the ice to be cold above 675 m a.s.l., the following relations are obtained using linear interpolation between 400 and 675 m a.s.l. (Fig. 3):

(3)

where h denotes the surface elevation and z _n the normalized ice thickness (ranging from 0 at the bottom of each ice column to 1 at the ice surface). The values for w thus obtained also correlate with the mean water content of 0.8 ± 0.26% found for the cold–temperate transition surface of the polythermal Storglaciären, Sweden, using relative backscatter strength of GPR signals (Pettersson and others, 2004).

Fig. 3. Water content w within a vertical column depending on the surface elevation h.

The enhancement factor m

The enhancement factor m provides a means to consider characteristics such as the anisotropy of the ice. Based on laboratory experiments, Wang and Warner (1999) proposed that the enhancement factor should be chosen depending on the actual stress field. Thus, a transition from a purely compression-dominated flow regime to a purely shear- dominated one should be considered. Applying this idea to an Antarctic flowline, Wang and Warner (1999) obtained good agreement between their 2-D flowline model results and observations. Since our numerical flow model considers three dimensions in space, we extended their relation as follows:

(4)

where ms represents the enhancement factor for pure shear and m _c that for pure compression. Additionally

ranges from 1 to 0 as the stress field varies from compression (S = 0) to shear (C = 0). In three dimensions, compression C and shear S can be written as follows:

This approach represents a somewhat crude method of including a dependence of the enhancement factor on the orientation distribution function (ODF) of the c axes. This is based on the simplified assumption that the ODF is mainly determined by the current stress field. Thus, Equation (4) describes implicitly a dependence on the ODF. At least for steady-state conditions, for which the stress field does not change over time, this approach has the potential to lead to reasonable results. However, the fact that the resulting flow law does not agree with the principle of material objectivity (frame indifference) due to its dependence on individual components of the stress tensor illustrates the limitations of this method. More appropriate ways of modelling the effects of induced anisotropy have been proposed theoretically (e.g. Svendsen and Hutter, 1996). Nevertheless, we applied the aforementioned approach to investigate the general dependence of ice dynamics on a stress-dependent enhancement factor. Another reasonable choice for the enhancement factor could be either a constant value (m = const.) or one that depends on the actual temperature field. For example, Sandhäger (2003) applied the exponential relation m(T*) =0.9 exp [0.17(T* + 10)] to the Larsen ice shelves. Since one of our goals is to investigate the effect of a temperature- dependent enhancement factor, we chose a Gaussian- shaped distribution with a stronger gradient for the temperatures just below 0°C, i.e.

(5)

(5)

The values within this relation were chosen such that the resulting values for m lie in the same interval as m(τ_ij) and m = 2. This is important if a comparison of the effect of the respective enhancement factors on the resulting velocity fields is undertaken. Usually the applied enhancement factors for grounded ice bodies range between 1 and 5 (e.g. Huybrechts, 1992; Paterson, 1994; Greve and others, 1998) and they are often used as a tuning parameter to fit existing in situ measurements of ice velocities. We chose a value of m = 2 since this best represents the few in situ measurements obtained during the 1997/98 field campaign (Braun, 2001).

Description of different model runs

To determine whether the aforementioned modifications lead to significant differences in the outcome of time- dependent model runs, we present the results of different diagnostic simulations with and without each modification.

First, we varied the threshold altitude below which the surface temperature is set to melting conditions and above which a linear temperature gradient is assumed. Model runs were undertaken for threshold altitudes at sea level (i.e. a cold-ice surface), at 25, 50, 75, 100, 125, 150, 175, 200, 300, 400, 500, 600 m a.s.l. and at the top of the ice cap (i.e. a temperate ice surface). For the discussion, we will mainly focus on the results of the model runs with the threshold altitude at sea level (named T000) and at 400 m a.s.l. (named T400) since the main differences can be discussed with the aid of these two threshold altitudes.

Second, to investigate the importance of water within the ice body, we simulated the ice dynamics, including and neglecting the water content. In the former case, we applied Equation (3). The names of these model runs end with numbers 2, 3 and 4. For runs ending with numbers 6, 7 and 8 we set w = 0% for all gridpoints.

Third, three different choices for the enhancement factor m were used. Results obtained using the most often applied constant enhancement factor m = 2 end with either 2 or 6. A stress-dependent enhancement factor as described by Equation (4) was taken for model runs ending with 3 or 7, and the endings 4 and 8 denote simulations with a temperature-dependent enhancement factor (Equation (5)).

An overview of the simulations presented here is given in Table 1.

Table 1. Overview of the different model runs presented in this study. First, the threshold surface elevation was varied. Below this altitude, the ice surface temperature T_s was set to pressure-melting point conditions T_s = 0°C, and above this limit we applied a linear lapse rate. Furthermore, the water content w and the enhancement factor m were adapted in different ways

Results and Discussion

The main purpose of the sensitivity study presented here lies in elucidating the effect of various modifications to model formulations on the resulting simulated ice flow of a temperate ice cap, similar to the one on KGI.

It is well known that the warmer the ice the faster it flows. But due to the non-linear character of ice described in Glen’s flow law, the interaction between temperature and flow velocity cannot be easily assessed. The 3-D flow model utilized in our study is thermomechanically coupled, enabling the calculation of internal ice temperatures as well as ice velocities and the internal stress field. Only the temperature at the ice surface needs to be given. As described above, we varied the threshold altitude for pressure-melting point conditions. The effects of the modifications are shown in Figure 4. It is obvious that the averaged horizontal velocity increases for warmer ice surfaces for all the presented parameter combinations. But this increase is not linear; only small differences occur in the mean value between simulations with threshold altitudes at 400 and 700m a.s.l.

Fig. 4. Averaged horizontal velocity for each model run. The x axis shows the threshold altitude above sea level below which the ice surface temperature was set to the pressure-melting point and above which a linear decrease in ice surface temperature was applied. The model runs with a temperature-dependent enhancement factor react most strongly to these different transition altitudes.

Additionally, it can be seen that increases in horizontal velocities for simulations accounting or not accounting for the internal water content have a similar shape. Nevertheless, the difference between the resulting mean horizontal velocities increases with increasing surface temperature although the water content does not depend directly on the temperature but only on the ice geometry (Equation (3)). This finding will be discussed in more detail while explaining Figure 9.

But the most obvious result shown in Figure 4 is the fact that the mean horizontal velocities of the model runs with constant enhancement factor increase similarly to those with a stress-dependent enhancement factor. In contrast, the simulations using a temperature-dependent enhancement factor react most strongly to an increasing surface temperature.

Considering only one characteristic value per model run is reasonable and necessary for a first impression of the trends. But due to the aforementioned non-linear character of the ice flow, vertical and horizontal distributions should also be considered.

Figure 5 shows the horizontal distribution of the product m(τ _ij , T*, …)A(T*, w, …) at the ice base. This product affects mainly the ice dynamics as can be seen from Equation (1). Both, the enhancement factor m and the Arrhenius factor A (through the water content) are affected by our modifications. Thus, when discussing the product m(τ _ij , T*, …)A(T*, w, …), the effect of the combination of two modifications can be examined.

Fig. 5. Horizontal distribution of the product m(τ _ij , T*, …)A(T*, w, …) at the ice base for the parameter combinations given in Table 1. In the first two rows, simulations with a limit for pressure-melting point conditions at the ice surface already at sea level (T000) are shown, whereas in the lower two rows, this threshold altitude was set at 400ma.s.l. (T400).

It can be seen that there is a difference not only in the absolute value but also in the pattern of the horizontal distribution of the product m(τjj, T*, …)A(T*, w, …) depending on whether the water content is considered or neglected. Comparing the first with the second, and the third with the fourth, row of the figure indicates this difference does not depend on the surface altitude alone. Additionally, the chosen value for the enhancement factor m is important. Furthermore, the threshold altitude for the surface temperature below melting conditions has smaller effects on the pattern of the resulting fields at the ice base.

This latter result is different at the ice surface (Fig. 6). Here, the absolute value of the product m(τ_ij , T*, …)A(T*, w, …) is much smaller than at the ice base. This can be explained by considering Equation (1): The product m(τ _ij , T*, …)A(T*, w, …) contributes to the deformation rate ɛ_ij which is largest at the ice base and smallest at the ice surface. Nevertheless, the distribution pattern of this product changes, as well as the absolute value for different threshold altitudes. But for the colder-ice case, the water content leads to smaller changes than in the warmer-ice case.

Fig. 6. Horizontal distribution of the product m(τ_ij , T*,.. ,)A(T*, w, …) at the ice surface for the same parameter combinations as shown in Figure 5.

In order to evaluate the differences in the ice dynamics of a vertical ice column, we compare horizontal ice velocities at the top and bottom layers of this column and vertically averaged horizontal velocities for different combinations of the enhancement factor and the water content. Figure 7 presents the results for models with a threshold altitude at sea level, and Figure 8 for models with pressure-melting point conditions at the ice surface up to 400ma.s.l.

Fig. 7. Correlation between horizontal velocities simulated with different parameter combinations for the enhancement factor m and the water content w. All results are obtained with the threshold surface elevation at sea level (i.e. T000). Left column: basal velocities; middle column: surface velocities; right column: vertically averaged horizontal velocities. Top row: constant m vs m(τ_ij ), both considering the water content; middle row: constant m vs m(T*), both considering the water content w; bottom row: considered vs neglected water content w, both with a temperature-dependent enhancement factor m(T*). The slope of the line is indicated as B on each subplot.

Fig. 8. Correlation between horizontal velocities simulated with the same parameter combinations as in Figure 7. Here the threshold surface elevation was 400ma.s.l. (i.e. T400). The slope of the line is indicated as B on each subplot.

First, considering the results for the colder-ice case (i.e. Fig. 7), it can be seen that the data points can be fitted with a linear trend but that many points differ from this fit. This means that the introduction of a stress-dependent or temperature-dependent enhancement factor leads to nonlinear changes in the pattern of the horizontal velocity field. This non-linearity is stronger for the temperature-dependent enhancement factor (middle row, correlation coefficient of R² = 0.977 for the surface velocities and R ² = 0.983 for basal velocities) than for the stress-dependent enhancement factor (top row, R ² = 0.996 and R ² = 0.994 at the ice surface and the ice base, respectively). In the former case, the resulting velocities are in general lower than those obtained with a constant enhancement factor (i.e. the slope B is smaller than 1; the exact value is given on each subplot) whereas in the latter case the ice flows faster (i.e. the slope B is larger than 1). The non-linear changes agree well with the findings obtained by Wang and Warner (1998) with a 2-D-flowline model comparing isotropic and anisotropic rheologies.

Regarding the bottom row of Figure 7, it can be seen that the water content generally increases the ice velocity (i.e. the slope B is larger than 1).

In all three rows, it is obvious that the velocities at the ice surface are more affected by changes than the velocities at the ice base. This implies that the modifications result in non-linear effects on the velocity field in the horizontal as well as in the vertical direction. This result was also shown by Wang and Warner (1998) for their 2-D-flowline model, though they neglected basal sliding in order to discuss only flow by internal deformation. Considering the last column of Figure 7, it appears that most deviations from the linear trend are reduced when calculating the vertical mean value of the horizontal flow velocity (e.g. the correlation coefficient for the upper-right subplot reaches R ² = 0.999).

Second, considering the warmer-ice case (i.e. Fig. 8), the deviations from the linear trend are much smaller than those seen for the colder-ice case presented above, i.e. the correlation coefficient for all subplots is larger than R² = 0.997. Especially in the middle row (m(T*) vs constant m), a nearly perfect fit is obtained (R² = 0.999) because, in both cases, the temperature within the ice body is almost everywhere identical to the pressure-melting point. This implies that the dependence on temperature is no longer different from a constant enhancement factor, except for a linear factor. As already mentioned above in the context of Figure 4, the simulations with a temperature-dependent enhancement factor react most strongly to an increasing surface temperature. This can be seen by comparing the middle rows of Figures 8 and 7. In the warmer-ice case (Fig. 8), the velocities obtained with a temperature-dependent enhancement factor are smaller (only around 94%) than those with a constant m, whereas in the colder-ice case (Fig. 7) it is vice versa (around 114%).

Up to now, we have discussed the changes obtained by the applied modifications regarding one characteristic value for the entire model domain (Fig. 4), the horizontal distribution of the product m(τ_ij , T*, .. .)A(T*, w, …) (Figs 5 and 6) as well as the correlation between different parameter combinations (Figs 7 and 8). Another interesting question is how the effects are distributed within one vertical ice column. We present this in the following using vertical slices through the ice body along the flowlines A and B indicated in Figure 2.

Effect of water content

In order to demonstrate the effects of internal water on the ice dynamics, we compare results for models that explicitly treat water content w in the calculations with those that neglect internal water. Figure 9 compares the ratios of horizontal ice velocities for models with and without internal water. Regions with a ratio larger than 1 are those where the resulting horizontal flow velocities for models with internal water are faster than those neglecting internal water. It can be seen that neglecting the water content leads to lower velocities in most parts of the ice body, except for small parts near the ice divide. Furthermore, our results demonstrate that the effect of water is highest at the top of the ice body and smallest at the ice base. Regarding Equation (3) which describes the distribution of the water content and Figure 3, the reverse is to be expected.

Fig. 9. The effect of the water content is shown using the ratio between two resulting horizontal velocity fields on vertical slices along flowlines A (left column) and B (right column). First row: constant enhancement factor m with threshold altitude at sea level (T000_2 vs T000_6). Second row: temperature-dependent enhancement factor m(T*) with threshold altitude at sea level (T000_4 vs T000_8). Third row: constant enhancement factor with threshold altitude at 400 m a.s.l. (T400_2 vs T400_6). Fourth row: temperature-dependent enhancement factor m(T*) with threshold altitude at 400ma.s.l. (T400_4 vs T400_8). The locations of the two flowlines are indicated in Figure 2. The grey is the underlying bedrock.

Another aspect that can be seen in the slices in Figure 9 is that the effect of the water content on horizontal ice velocity is larger in the colder-ice case (upper two rows) than in the warmer-ice case (lower two rows). Regarding the colder-ice results, there are parts in the ice body where the flow velocity with respect to the water content reduces to 83% of the velocity field obtained without water. In other parts of the same slice, the velocity obtained with internal water content reaches 111% of the no-water velocities. In the warmer-ice case, these intervals lie between 100% and 110%. These findings confirm the result seen when comparing Figures 7 and 8, that the correlation between the data points and the linear fit in the warmer-ice case is better than in the colder-ice case. The differences in resulting ice velocity within the ice body due to the water content are distributed more evenly in the warmer-ice case than in the colder-ice case.

Effect of enhancement factor

Figure 10 also shows ratios in horizontal ice velocities along flowlines A and B. In this case, we discuss the effect due to the choice of the enhancement factor. That is, the ratios in horizontal velocities from model runs with a temperature-dependent enhancement factor (first and second row) or a stress-dependent enhancement factor (third and fourth row) against those obtained with a constant enhancement factor are calculated. Values larger than 1 indicate that the constant enhancement factor leads to lower velocities than the stress- or temperature-dependent enhancement factor.

Fig. 10. The effect of the choice of enhancement factor is shown using the ratio between two resulting horizontal velocity fields on the same vertical slices along flowlines A (left column) and B (right column) as in Figure 9. First row: threshold altitude at sea level both for temperature-dependent enhancement factor m(T*) (T000_4) vs constant enhancement factor m (T000_2). Second row: the same choices for the enhancement factor for the threshold altitude at 400 m a.s.l. (T400_4 vs T400_2). Third row: threshold altitude at sea level both for stress- dependent enhancement factor m(τ _ij ) vs constant enhancement factor m (T000_3 vs T000_2). Fourth row: the same choices for the enhancement factor for the threshold altitude at 400 m a.s.l. (T400_3 vs T400_2).

Regarding the effects of applying a temperature-dependent enhancement factor (upper two rows), it is obvious that for most regions in the colder-ice case the horizontal velocities are smaller (down to 64%) than applying a constant enhancement factor. Only at the ice front, the velocity for m(T*) is up 5% for slice A and 14% for slice B larger than with a constant m. In contrast, for the warmer-ice case, all regions show higher velocities for m(T*) than for a constant m (103–115%). This change was noted for Figure 4, which shows a stronger increase of the averaged value for the horizontal velocities applying the temperature-dependent enhancement factor than for the two other possible choices for m.

Considering the stress-dependent enhancement factor vs the constant enhancement factor, all parts of the ice body show higher velocities for m(τ _ij ) than for constant m. The intervals vary between 103% and 137%.

Comparing now the two colder-ice cases, it is obvious that the temperature-dependent enhancement factor leads to larger differences regarding the absolute interval as well as the horizontal and vertical distributions of the ratios. This can be explained by the fact that Wang and Warner (1999) applied the stress-dependent enhancement factor to their model domain to describe the anisotropy of the ice body. These anisotropy effects are larger for cold ice and reduce if the ice temperature reaches pressure-melting point which is the case in our study region.

In the warmer-ice cases (second and fourth rows), the absolute values of the ratios are smaller than in the colder-ice case, but their horizontal and vertical distribution stays non-linear. As already mentioned within the discussion about the water content, it can also be seen here that the differences of the horizontal flow velocities are largest at the ice surface and smallest at the ice base.

Effect of surface temperature distribution

It has been mentioned several times that the modifications to the numerical model have larger effects in the colder-ice case than in the warmer-ice case. We now compare these two cases. Figure 11 shows the ratio of the horizontal velocities obtained for the warmer-ice case against results for the colder-ice case. First, the expected fact that warmer ice flows faster than colder ice can be observed. But this difference is much larger in the higher altitudes (up to 398%) than near the ice front. This is because the ice temperature at sea level is at, or close to, the pressure-melting point conditions in all cases. Thus the temperature increase mainly affects the upper regions. The same is true for the ice base, which is almost not touched by changing the threshold altitude. Comparing the two rows in Figure 11, it can be seen that the differences of the horizontal flow velocities are larger applying the temperature-dependent enhancement factor than using the constant enhancement factor.

Fig. 11. The effect of the threshold surface elevation for the temperature at or below pressure-melting point is shown using the ratio between two resulting horizontal velocity fields on the same vertical slices along flowlines A (left column) and B (right column) as in Figures 9 and 10, but with a different colour scale. First row: constant enhancement factor m both for a threshold altitude at 400 m a.s.l. (T400_2) vs one at sea level (T000_2). Second row: the same choices for the threshold altitudes but for a temperature-dependent enhancement factor (T400_4 vs T000_4).

Conclusions

We modified a 3-D, thermomechanically coupled ice-flow model formerly applied only to cold-ice regions on mainland Antarctica (e.g. Paschke and Lange, 2003; Sandhäger, 2003) in order to consider polythermal and temperate ice. The effects of the different modifications on the numerical model results are presented and discussed for the case of an ice cap similar to that on western KGI.

The three quantities that were modified are the enhancement factor m which is chosen either to be constant or to depend either on temperature or on stress conditions in the ice body; the choice of considering or neglecting the internal water content w of the ice cap in Glen’s flow law; and the threshold altitude up to which the ice surface temperature T _s is set to pressure-melting point conditions.

We find that:

1. 1. The inclusion of the water content leads to higher flow velocities.

2. 2. When considering the temperature dependence of the enhancement factor m(T*), smaller (higher) velocities are obtained in colder-ice (warmer-ice) cases compared to results from model runs with a constant enhancement factor. This implies that the averaged velocity of all gridpoints increases faster with the threshold altitude for the temperature-dependent enhancement factor simulations than for those with constant or stress-dependent enhancement factors.

3. 3. Considering the dependence of the enhancement factor m(τ _ij ) on the stress field leads to small increases in the horizontal flow velocity compared to results obtained with a constant enhancement factor. These differences are larger for the colder-ice cases than for the warmer-ice cases. This might be explained by reduced anisotropy for ice reaching pressure-melting point conditions.

4. 4. The choice of the threshold altitude for pressure-melting point conditions at the ice surface has larger effects in the upper part of the ice body since the region close to the ice front can be considered as temperate ice in all cases. In general, increasing threshold altitude leads to higher flow velocities.

5. 5. The differences of the flow velocities as well as of the product m(τ _ij , T*,.. .)A(T*, w, …) are not entirely linear. They vary with the horizontal as well as the vertical position within the ice body. This implies that these changes cannot be described by just using any constant factor within Glen’s flow law.

All these findings underline the necessity of modifying cold- ice models before applying them to a temperate ice body. The observed differences are larger for the colder-ice cases, i.e. an ice body that is not entirely temperate but includes sections with ice temperatures between –5 and 0°C. This is certainly the case for the KGI ice cap where so-called cold spots are expected near the ice divide, based on radar measurements. Additionally, ice caps on the Antarctic Peninsula may also be suitable to be treated with our modified model.

We showed that it is crucial to carefully choose the respective parameters before simulating the ice dynamics of temperate and polythermal ice bodies. The presented findings will be used in future time-dependent model runs when predictions about the impacts of climate change on temperate and polythermal ice caps or glaciers are required. Nevertheless, a general description of which parameter combination is the best cannot be given for all temperate or polythermal ice caps. This sensitivity study is more intended to describe the possible effects that one or other parameter combination can have on the simulated patterns of the respective velocity fields.