A full Stokes-flow thermo-mechanical model for firn and ice applied to the Gorshkov crater glacier, Kamchatka

Abstract

The Gorshkov crater glacier at Ushkovsky volcano, Kamchatka, is characterized by a large aspect ratio and special thermodynamic conditions at the bedrock caused by a locally enhanced and spatially varying geothermal heat flux. Furthermore, large parts of this glacier consist of firn rather than pure ice, which alters the rheological properties (such as viscosity and compressibility) of the glacier. We present a newly developed, thermo-mechanically coupled, three-dimensional flow model based on the finite-element (FE) modeling software Elmer, and apply it to the Gorshkov crater glacier. By assuming steady-state conditions, the present-day velocity field, temperature field, basal melting rate and age distribution are simulated. We find that flow velocities are generally small (tens of centimeters per year). Horizontal and vertical velocities are of comparable magnitude, which shows that the shallow-ice approximation is not applicable. Owing to the spatially variable volcanic heat flux, the thermal regime at the ice base is cold in the deeper parts of the glacier and temperate in the shallower parts. The measured temperature profile and age horizons at the K2 borehole are reproduced quite well, and remaining discrepancies may be attributed to transient (non-steady-state) conditions. Firn compressibility is identified as a crucial element for the modeling approach.

1. Introduction

Glaciers which develop in volcano craters are unique systems because of their particular morphologies (large aspect ratio) and thermodynamic conditions (large geothermal heat flux). In this study, we investigate the Gorshkov crater glacier at Ushkovsky volcano, Kamchatka (56804’N, 160828’E, 3903 ma.s.l.), which is the only crater glacier for which the bedrock topography has been determined so far. The glacier is situated in the summit caldera of Ushkovsky volcano and fills the concave bed of Gorshkov crater to a maximum depth of 240 m. The crater diameter is approximately 750 m. The resulting aspect ratio of ≈1/3 clearly indicates that the shallow-ice approximation (SIA) (Hutter, 1983; Morland, 1984; Greve, 1997) is not applicable. The glacier surface is gently inclined towards the northern crater rim, where the ice flows out of the crater and down the slope of the volcano. The geothermal heat flux at the crater rim is estimated to be as large as 10Wm^–2 based on direct measurements (Murav’yev and Salamatin, 1989), whereas the temperature profile of the K2 borehole close to the deepest point of the glacier indicates a far smaller heat flux of only 0.12 Wm^–2 at this position (Shiraiwa and others, 2001).

An analytical, thermo-mechanically coupled model for ice flow and heat transfer in crater glaciers was developed by Salamatin and others (2000) based on a scaling analysis which accounts for the particular geometry. The authors assumed a linear-viscous rheology for porous glacier ice and applied their model to a flowline profile of Gorshkov crater glacier from the summit of the southern crater to the outflow at the northern rim. They estimated that basal melt is up to 25% of the accumulation rate at the ice surface, and predicted a maximum age of the crater-glacier ice of more than 610 years. An attempt to simulate the dynamics of Gorshkov crater glacier numerically was made by Shiraiwa and others (2002) based on the higher-order model by Blatter (1995). In their study, the ice was treated as an incompressible, isothermal, viscous fluid, and the velocity field was computed along the flowline mentioned above. The authors found flow velocities of the order of 1–2ma^–1, generally increasing toward the northern outflow face, and demonstrated the non-applicability of the shallow-ice approximation.

Here, we refine this approach by applying a newly developed, thermo-mechanically coupled, three-dimensional flow model of the glacier. This model is based on the finite-element modeling software Elmer, developed at CSC–Scientific Computing Ltd. (Espoo, Finland), and solves diagnostically the full Stokes-flow equations for a non-linear viscous fluid. Further, since large parts of the glacier consist of firn rather than pure ice, we have implemented a compressible firn rheology and a suitable expression for the heat conductivity of firn. Our objective is to model diagnostically the three-dimensional velocity, temperature and age fields as well as the basal melting rate of the present-day glacier. Since no surface velocity measurements are available so far, comparison of computational results against data is limited to the measured temperature profile and some age horizons at the K2 borehole.

2. Topography of the Gorshkov Crater Glacier

2.1. Digital elevation model

A geodetic and radio-echosounding survey was conducted on the Gorshkov crater glacier (Shiraiwa and others, 2002). Since the data have only been accessible to us as surface and bedrock topography maps, we have manually digitized them in order to assemble digital elevation models of the ice surface and bedrock at a horizontal resolution of 20 m. The result is depicted in Figure 1, which also indicates the position of the drill site K2. A Cartesian coordinate system is defined, where the x, y and z axes coincide with east, north and vertically upwards. The corresponding unit vectors are e _x , e _y and e _z . The surface and bedrock elevations are denoted as z = h ^s(x, y) and z = h ^b(x, y) respectively.

Fig. 1. (a) Surface topography of the Gorshkov crater glacier (greyscale texture and contours, contour spacing 1 m); (b) Bedrock topography (contour, spacing 25 m) and basal relative density (greyscale texture).

2.2 Computational mesh

Using the commercial pre-processing tool Gambit, an unstructured mesh of the glacier domain has been created. It consists of approximately 4×10⁴ nodes, connected to form about 2×10⁵ linear tetrahedral volume elements and 2×10⁴ triangular boundary elements (Fig. 2). The feature of built-in size functions in Gambit is used to introduce mesh refinements towards the free surface to capture compaction effects of the firn, according to Equation (3) below. Additional increased resolutions of the mesh around the position of the drill site K2 as well as the bedrock are applied. Whereas the first of these refinements is dedicated to improvements of the post-processing steps (i.e. interpolation of profiles along the ice core), the latter resolves the large temperature gradients which occur in the case of a temperate ice base.

Fig. 2. Cut through the computational mesh of the Gorshkov crater. The black vertical line indicates the position of the borehole K2. The three boundary areas, namely the free surface, bedrock and outflow, are indicated.

3. Thermo-Mechanically Coupled Model Equations

3.1. Density profile

The volume fraction of ice in firn (or relative density) is defined as

(1)

where ρ is the variable density of the firn and ρ _i = 910 kgm^–3 is the density of pure ice. The contribution of the interstitial air to the mass is neglected. The volume fraction of ice and the porosity (volume fraction of air) _a total unity, so that

(2)

At the K2 drill site, close to the position of the maximum ice thickness, the measured density profile can be approximated by the exponential function

(3)

where

(4)

is the depth from the glacier surface, φ _a,s = 0.55 is the porosity at the surface and 훾 = 0.038 m^–1 is the compaction coefficient (Shiraiwa and others, 2001). Following Sorge’s Law, which states that the density field is in steady state and only a function of depth below the ice surface (Bader, 1954), we assume that this relation is valid throughout the whole glacier domain. Figure 1b shows the resulting distribution of the relative density at the glacier bed. Evidently, the basal density is close to that of pure ice everywhere except for the near-marginal areas where the ice is thin.

3.2 Flow law

Let t be the Cauchy stress tensor which acts on an arbitrary volume element in the firn. Like every tensor, it can be decomposed into an isotropic part –p1, where p = –(tr t )/3 is the pressure (tr denotes the trace operator i.e. the sum of the main diagonal elements) and 1 the unit tensor, and a traceless, deviatoric stress tensor, i.e.

(5)

Similarly, the deviatoric part D ^D of the strain-rate tensor D = 1/2[grad v + (grad v)^T] (symmetric part of the gradient of the velocity v; the superscript T denotes tensor transposition) is obtained by

(6)

Following the lines of Gagliardini and Meyssonnier (1997) with the definition of the invariant

(7)

we express the relations between the deviatoric and

isotropic parts, respectively, of the stress and strain rate as:

(8)

and

(9)

The quantity ɳ(T′, δ) introduced in Equation (8) is referred to as (shear) viscosity. The dependency of a(φ)and b(φ) upon the relative density has been obtained from field data following Gagliardini and Meyssonnier (1997),

(10)

(see also Fig. 3). As the stress exponent n = 3 is applied, in the limit of pure ice φ →1 so we have a→1, b→0 and div v→0, and thus the widely used Glen’s flow law (Paterson 1994) is restored. Note also that Equation (9) relates the total pressure to the divergence of the velocity field (rate of volume change), contrary to the classical description of a compressible viscous fluid, where a thermodynamic (velocity-independent) and a viscous (velocity-dependent) contribution to the pressure can be identified. This implies that, as long as φ< 1, a hydrostatic equilibrium in the presence of gravity does not exist for firn.

Fig. 3. Dependency of the dimensionless functions a and b and density ρ upon the relative density φ

The flow rate factor A(T′) is derived from the Arrhenius law,

(11)

where temperature relative to the pressure melting point T′= T – T _m + T ₀, where A ₀ is the pre-exponential factor, Q the activation energy, R the universal gas constant (8.314 J mol^–1 K^–1), T the absolute temperature and T ₀ = 273.16 K the melting point of ice at low pressure. The pressure dependence of the pressure melting point T _m can be expressed by the linear relation T _m(p) = T ₀ – βp where the Clausius-Clapeyron constant β = 9.7456× 10^–8KPa^–1. Following Paterson (1994), we express the

pre-exponential factor as

(12)

and the activation energy as

(13)

which yields a continuous connection of the two regimes with A = 4.9×10^–25 s^–1 Pa^–3 at T′ = –10° C. The flow enhancement factor E is set by default to E = 1/3. The deviation from unity has been found to be necessary, as otherwise the upper, low-density firn region tends to yield with respect to compression. This results in too large a downward convection and consequently too low temperatures and young ice in the lower regions of the glacier.

3.3 Thermodynamic material equations

In pure ice, the specific heat (J kg^–1 K^–1) can be expressed as a function of temperature via

(14)

where c ₀ = 2127:5 J kg^–1 K^–1 and cΔ = 7:253 J kg^–1 K^–1. A representation of the thermal conductivity ⁱ (Wm^–1 K^–1) is

(15)

where ₀ = 9:828Wm^–1 K^–1 and γT= 0:0057 K^–1 according to Ritz (1987). For firn, the specific heat of interstitial air is negligible, so that

(16)

A relation linking the heat conductivity of firn κ_f(Wm^–1 K^–1) to the relative density has been provided by Sturm and others (1997):

(17)

where f₀ = 0:138Wm^–1 K^–1, _f1 = 1:010 × 10^-3Wm³ kg^–1 K^–1 and _f2 = 3:233 × 10^-6Wm⁵ kg^–2 K^–1 with the density ρ(’) evaluated using Equation (3). In order to include both the temperature and the density influence, we combine (15) and (17) and obtain

(18)

The heat flux then follows from Fourier’s law of heat conduction

(19)

and the specific internal energy is given by the caloric equation of state

(20)

In the differential form of Equation (20), the dot denotes the material derivative with respect to time t following the motion of the firn particles.

3.4 Field equations

Neglecting the acceleration terms in the momentum balance yields the Stokes equation

(21)

where g is the acceleration due to gravity (9.81ms^–2 in the negative z-direction), and the pressure p and viscosity ɳ are as given in Equations (8) and (9). In addition, by combining Equations (8) and (9), we obtain the volume balance

or

(22)

where к _cp is the volume-pressure coupling, dV an incremental, material volume, and the dot denotes the material time derivative as in Equation (20). The volume-pressure coupling vanishes in the limit ’!1, restoring incompressibility for pure ice.

From the energy balance, Fourier’s law of heat conduction (Equation (19)) and the caloric equation of state (Equation (20)), the heat-transfer equation can be derived:

(23)

The production term D: t is the strain heating. Volumetric heat sources due to radiation etc. are negligible. Since we only consider steady-state conditions in this paper, the advection (left-hand side) is balanced by conduction and strain heating (right-hand side). Note that the ice temperature is limited by the pressure melting point T _m, so that the solution of the heat-transfer equation is subject to the constraint

(24)

The age of the ice A is convected with the flow. For a given velocity field v the dating equation is of the form

(25)

As none of the material parameters depend on the age of the ice (which would be the case if impurity layers, such as ash from a volcanic eruption, were included in the rheological behavior) there is no feedback to either the Stokes or the heat transfer equation. Consequently, it can be run in a post-processing step. Since steady-state of the flow field is assumed, the first term on the left-hand side, әA/әt is dropped.

3.5 Boundary conditions

As depicted in Figure 2, boundary conditions for three different areas must be provided. These are the free surface, the bedrock and the outflow area at the northern margin of the glacier.

3.5.1. Surface

Since we assume steady-state conditions and only carry out diagnostic simulations with prescribed geometry, it is neither required to prescribe the surface mass balance (accumulation-ablation rate) nor a kinematic boundary condition for the evolution of the free surface. The stress-free condition

(26)

serves as a dynamic boundary condition, where n is the outer normal unit vector (pointing into the atmosphere), and the subscript s denotes values taken at the free surface. The surface temperature is kept at the constant value T∣_s = -16:6°C which represents the average temperature at a depth of 1.79 m at the drill site B1 over a period of two years measured in 1995/96 (Shiraiwa and others, 2001).

The boundary condition for the dating equation is given by A∣_s = 0, representing fresh snow falling on the free surface. As the free surface is the only inflow boundary of the volume, no other boundary condition for this variable has to be set.

3.5.2. Bedrock

The geothermal heat flux, which enters the ice body from the underlying bedrock, is assumed to be lowest for the deepest areas and to increase exponentially with a previously undetermined factor m with increasing bedrock elevation h _b,

(27)

The minimum and maximum bedrock elevations are = 3640m and = 3900 m. The lower limit of the geothermal heat flux, = 0.12Wm^–2, is chosen to reflect the measured temperature gradient at the bottom of the drilled ice core at K2 (Shiraiwa and others, 2001) whereas the highest assumed value, = 10Wm^–2 is set according to earlier estimations (Murav’yev and Salamatin, 1989) that have recently been validated by direct measurements at the rim of the crater (Murav’yev, personal communication, 2006). As the exponent m was initially unknown, runs with different values, starting from m = 1 with increments of unity up to m = 7 have been conducted. From the comparison of results with estimations of the bedrock area covered by ice below the pressure melting point (‘cold base’) (Shiraiwa and others, 2001), the most reasonable range of this parameter could be constrained to 2 ≤ m ≤ 6. In the case of a linear distribution (i.e. m = 1) the entire bedrock is at the pressure melting point (‘temperate base’), and the other extreme value (m = 7) leads to results with cold-base conditions covering too much of the bedrock up to the rim. A detailed discussion of the results obtained with parameter values m = 2, 4 ,6 follows in Section 5.

In case of a cold base the geothermal heat flux determines the normal derivative of the basal temperature. However, if the pressure melting point is reached, basal melting sets in and the basal melting rate M ^b is

(28)

where L = 3.335×10⁵ J kg^–1 is the latent heat of ice and n the outer normal unit vector (pointing into the bedrock). As no indication of meltwater storage at the bedrock from the K2 drilling site could be found, we assume complete drainage. As a consequence, the volume flux M _b /ρ equals the normal velocity v|_b . n at the bottom (note that the outer normal unit vector n points into the bed). Further, we impose a constraint of positive normal velocities (i.e. re-freezing of ice is neglected), so that

(29)

As for the tangential component of the basal velocity v|_b, we employ no-slip conditions: that is, this component is assumed to be negligible.

3.5.3. Outflow

In view of the lack of data on the mass flux and/or the velocity distribution of the outflow zone, the outflow velocity is prescribed directly,

(30)

where n e _y is the outer normal unit vector in the horizontal plane (pointing away from the glacier). We use the fixed value v _out ≈ 0.39ma^–1 which has been determined from a balance of the average accumulation mass flux at the free surface, j ^s = 570 kgm^–2 a^–1 (Shiraiwa and others, 2001), and an average of the basal melting rate (Equation (28)) obtained from an initial run (with no-flux conditions imposed), via

(31)

Here, the areas of the free surface, the bedrock as well as the outflow are denoted by the symbols A _s, A _b and A _out respectively. We have chosen this approach as alternative options of prescribing an external normal stress lead to exaggerated accumulation rates at the free surface. Extension of the computational domain beyond the crater into the ice-covered caldera as well as including the evolution of the density profile would reduce the uncertainty introduced by condition (31).

For the thermodynamic boundary condition, we assume the no-flux condition

(32)

4. Finite Element Model Elmer

The model equations detailed in the previous sections are solved numerically with the open-source multi-physics package Elmer (see http://www.csc.fi/elmer) which is based on the finite element method (FEM). Since the code takes care of scaling internally, the simulations are conducted using SI units irrespective of the numerical values of the several quantities.

The system is initialized by computing the flow depth d for all nodes, which is non-trivial due to the unstructured grid. This is carried out by applying a standard Galerkin method (Girault and Raviart, 1986) to the elliptic boundary-value problem

(33)

which simply reflects Equation (4) (subscripts s and b refer to the surface and the base, respectively). Thereafter, the relative density ’(d ) is set according to Equation (3).

The iteration loop consists of solving the heat-transfer Equation (23) and the Stokes Equation (21). By integrating Equation (23) with a test function, over the whole volume V and applying Fourier’s law of heat conduction (Equation (19)), we obtain the weak formulation of the problem

(34)

In this case we choose piecewise linear functions, i.e. ϕ ∊ P ₁, which are also used for discretization of the unknown temperatures, T = Σ_iϕ_i T _i leading to the standard Galerkin formulation of the problem. After discretization and linearization, the two integrals on the left-hand side contribute to the system matrix M. The surface integral on the right-hand side introduces Neumann type (i.e. external heat flux) boundary conditions into the system. Since the solutions of the flow and the temperature field are treated sequentially, the integral over the heat production can be identified as the force vector F. With the solution vector T containing all nodal values of the unknown temperatures, the linearized system takes the simple form

(35)

As the resulting diffusion-convection equation shows tendencies towards instability for large velocities, the Stabilized Finite Element Method (Franca and Frey, 1992) is used to guarantee convergence.

The system given by Equation (35) still lacks the constraint imposed by Equation (24). We would like to stress the fact that the approach of simply resetting all entries of the solution vector exceeding T _m obtained by solving the unconstrained system (35) does not lead to a correct solution for the temperature field. A consistent method of introducing the inequality (24) is described as follows.

Nodes with values from the solution vector of the (k–1)th iteration T ^(k–1) showing a value larger than or equal to the constraint, T_i ^(k–1) ≥ T _m, are marked as ‘active’ (starting from a complete set of ‘inactive’ nodes at the first iteration). Here T_i ^(k–1) is the i th component of the solution vector T ^(k–1).

An ‘active’ node is switched back to ‘inactive’ only if its residual is (in case of an upper limit being imposed) larger than zero (R _i > 0).

The system described in Equation (35) is assembled for the new iteration step k: M. T ^(k) = F. If the i th node is identified as ‘active’, the system is manipulated such that all entries of the i th row in the system matrix are set to zero, except for the diagonal entry M_ii to which the unit value is assigned. Simultaneously, the i th entry of the force vector is set to the local pressure melting point, F_i = T ^m. From a numerical point of view this is equivalent to setting a Dirichlet condition in the bulk of the domain, leading to the manipulated system

(36)

for which the solution T ^(k) is computed.

The residual vector is then obtained by using the original system matrix and force vector

(37)

If the i th node has been marked as ‘active’, its component of the residual vector R_i can be interpreted as the additional cooling needed to comply with Equation (24).

The iteration above is continued as long as convergence of the solution is achieved. Compared to the previously applied Uzawa method (Girault and Raviart, 1986), this new algorithm proves to be faster and more robust.

The discretization of the Stokes Equation (21) is equivalent to that demonstrated for the heat transfer equation. Contrary to the earlier method, stabilization is introduced by the Method of Residual Free Bubbles (Baiocchi and others, 1993).

The final step consists of solving the dating Equation (25). Since it is convection-dominated, it requires a different numerical approach. To this end, the Discontinuous Galerkin Method for convection-reaction equations (Brezzi and others, 2004) has been applied. This approach treats the system locally for each element and treats the resulting boundary terms such that stability as well as consistency is maintained.

The flowchart of a complete steady state solution is given in Figure 4.

Fig. 4. Flowchart of the numerical simulations. Symbols behind leftwards pointing arrows indicate the needed input variables, whereas rightwards pointing arrows depict solution variables for the different modules.

5. Results and Discussion

For the reference run (m4), we have applied the compressible firn model and the rheological parameters as given in Section 3, as well as a geothermal heat-flux exponent of m = 4 (see Equation (27)). In addition, in order to obtain some information about the robustness of our results, two series of sensitivity studies were completed.

In the first series of tests, we varied the values of the parameter m, heat-flux exponent, in Equation (27). In run m2, all other variables were set as for m4, but with m = 2. In run (m6), we set m = 6.

The second series of tests involved variation of the rheological behavior. For run (a), parameters were as for (m4) but b(’) 0 (i.e. ‘incompressible firn’). For run (b), parameters were as for (m4) with the exception of a(’) 1, b(φ) 0 and E = 1/(3×’) (i.e. ‘porous ice’). For run (c), the same parameters as (m4) were employed except for a(’) 1 and b(φ) 0 (i.e. ‘pure ice’).

Figure 5 shows the horizontal and vertical near-surface velocities of the reference run (m4), and Figure 6 shows the velocity field for a S-N transect through the K2 borehole. Naturally, the viscous glacier ice is drained towards the northern outflow. Horizontal velocities are largest where the glacier is thick (compare with Fig. 5), whereas the largest vertical (downward) velocities occur further upstream in the southern part of the glacier due to firn compaction and intensive basal melting. Absolute values of the velocity are of the order of 10’s of centimeters per year and therefore quite small compared to typical valley glaciers. This is partly a consequence of the gently inclined surface, and partly due to the drag of the trough-shaped bedrock, transferred into the ice body by resistive longitudinal stresses. Note also that the horizontal and vertical velocities are generally of the same order of magnitude, which clearly demonstrates the failure of the shallow-ice approximation for the given geometry.

Fig. 5. Near-surface velocity (at 10m depth) for the reference run (m4) (a) absolute value of the horizontal (x–y plane) velocity; (b) absolute value of the vertical velocity (z-direction).

Fig. 6. Velocity for the reference run (m4): vertical transect in S-N direction through the K2 borehole. Note that no vertical exaggeration has been applied.

The thermal conditions at the glacier bed are shown in Figure 7. For the reference run (m4), the area of cold basal conditions makes up approximately one-quarter of the glacier area. The transition line between cold and temperate basal ice follows roughly the 3765 m contour of the basal topography in the interior of the glacier, whereas close to the outflow face, cold ice also occurs at higher bedrock elevations. As to be expected, the extent of basal temperate ice strongly depends on the assumed value for the heat-flux exponent. For run (m2), the area of cold basal ice is reduced to a small spot situated in the very deepest part of the bed trough. By contrast, for run (m6), the area of cold basal ice extends approximately 50 m further up the bed trough compared to run (m4). Lacking the downward convection of cold temperatures caused by the compressibility of firn, the runs (a), (b) and (c) produce a temperate base for the entire glacier area.

Fig. 7. Regions of cold and temperate basal ice (inside and outside of the white contour lines, respectively) for the reference run (m4), the run (m2) and the run (m6). The bedrock elevation is underlaid as greyscale texture.

Computed profiles of the temperature and age of the ice as well as the northward and vertical velocity components at the position of the K2 borehole are shown in Figure 8. The measured temperature profile T(d) at K2 (Shiraiwa and others, 2001), against which we check the computed profiles, has been approximated by the piecewise-linear function

Fig. 8. Profiles of simulation results at the K2 drill-site: (a) temperature; (b) age profile and (c) vertical and northward velocity. Measured temperature profile and age points by Shiraiwa and others (2001).

(38)

Extrapolated to the base (d = 240 m), this yields a basal temperature of –2.768C. Also, four age horizons are available for the K2 borehole, which have been dated by assigning ash layers to known volcanic eruptions (Shiraiwa and others, 2001).

Comparing the results of runs (m2), (m4) and (m6), the temperature profile at K2 becomes warmer with decreasing heat-flux exponent m, and run (m2) even produces a temperate base at K2 (see also Fig. 7). This is not surprising, because a smaller value of m means that the geothermal heat flux increases faster with bed elevation, and therefore generally warmer conditions prevail. The shape of the three simulated profiles is slightly concave from the top down to approximately two-thirds of the depth, which is the normal behavior for steady-state conditions in the accumulation area of a glacier due to downward advection of cold firn from the surface (Paterson, 1994). However, the measured profile does not show any concavity in this part, and is therefore up to 28C warmer. This may be a hint that the real glacier is not in a steady state.

The simulated age profiles agree quite well with the measured horizons, in that the four data points are within the range of age profiles of the three runs. For the three uppermost horizons, the agreement is best for the reference run (m4), whereas the lowest horizon (year 1400) is closer to the profile of run (m2). This may also be a consequence of non-steady-state conditions. We are not able to make a prediction for the age of the basal ice, because the combination of steady state, a cold base and no-slip makes the basal age unconstrained.

The velocity profiles of the three runs are similar. The horizontal, northward velocities (towards the outflow face) show the typical shear-flow profile and reach maximum values of slightly less than 0.5ma^–1 at the surface. The vertical velocity profiles are almost linear below 20 m depth. By contrast, at shallower depths the gradient is much larger, which is an immediate consequence of the large firn compressibility for small densities. As discussed above, the magnitude of horizontal and vertical velocities is comparable.

An important lesson to learn is that accounting for the compressibility of firn is crucial. Runs (a), (b) and (c) all have the property b(φ) = 0, which expresses compressibility in the flow law (Section 3.2), and the resulting temperature, age and velocity profiles at K2 fall virtually together. It is evident that the simulated temperatures are too high, and the simulated ages are far too large; the latter resulting from strongly reduced vertical velocities as a direct consequence of the now incompressible rheology. Therefore, the results of these three runs must be discarded as unrealistic, and we conclude that a compressible firn rheology is essential for obtaining reasonable simulation results.

Conclusion

A new dynamic-thermodynamic model for glacier flow has been presented. It is based on the finite-element package Elmer (CSC–Scientific Computing Ltd., Espoo, Finland) and solves the full Stokes-flow and heat transfer equations for a non-linear viscous, compressible fluid in three dimensions. The model has been applied to the Gorshkov crater glacier, which is characterized by a large aspect ratio between depth and diameter, a relatively low surface temperature and intensive volcanic heating from below. Under the assumption of steady-state conditions, we have simulated the present-day velocity field, temperature field and age distribution. Flow velocities are generally small – of the order of tens of centimeters per year – and of comparable magnitude in the horizontal and vertical directions. As a consequence of the highly variable geothermal heat flux, the basal ice is cold in the deeper parts of the glacier and temperate in the shallower parts. The measured temperature profile and age horizons at the K2 borehole have been reproduced quite well, even though some discrepancies point out that the steady-state assumption may be simplistic. Firn compressibility was identified as a crucial element for the modeling approach.

Future work will comprise simulating the density profile instead of simply prescribing it. In order to do so, the model equations need to be complemented by a three-dimensional, compressible mass balance. Also, it is highly desirable to conduct a new field campaign on Gorshkov crater glacier and obtain measurements for the surface velocity at different locations. Comparison of our results against data would be very helpful for further constraining the simulation set-up. As a long-term objective, simulations of non-steady-state scenarios over the last 1000 years or so should be striven for. However, this introduces the need for the climatic forcing and basal conditions of the glacier, and it is largely unclear how these input fields can be provided. A possible way to tackle this problem is by applying methods of inverse modeling, provided that sufficient data of the state of the glacier are available.