Introduction
Over the past decade it has become evident that subglacial water is widespread beneath the Antarctic ice sheet (Reference SiegertSiegert, 2000; Reference Siegert, Carter, Tabacco, Popov and BlankenshipSiegert and others, 2005; Reference Smith, Fricker, Joughin and TulaczykSmith and others, 2009) and affects the ice dynamics on sub-annual timescales (Reference Fricker, Scambos, Bindschadler and PadmanFricker and others, 2007; Reference Stearns, Smith and HamiltonStearns and others, 2008; Reference Fricker and ScambosFricker and Scambos, 2009). These observations dealt with water on horizontal scales of kilometres and vertical scales of tens of centimetres to metres (i.e. lakes and hydraulic connections between lakes). Here we focus on the effect of a thin water layer on the dynamics of ice sheets.
Underneath a thick ice sheet the basal temperature locally reaches the pressure-melting point and melt rates of the order of millimetres per year may generate a subglacial water layer. This basal water lubricates the base. A thought-experiment helps us to understand the consequences. Let us begin with an ice sheet frozen to the base. Large ice deformation at the base produces heat (strain heating) when ice is frozen to the base. The heating may lead to a temperate ice layer at the base. In this temperate layer the deformation is enhanced due to the reduced viscosity caused by microscopic water inclusions; temperate ice is soft. Ice with a temperate base is lubricated by subglacial meltwater, which enhances the sliding velocity. Basal sliding additionally leads to frictional heating at the base. A further consequence of the initiation of sliding is that the reduction of shear decreases the heat source of internal strain heating and, over time, the basal layer cools. When heat transport by geothermal and frictional heat into the ice does not counterbalance this cooling, the ice again becomes frozen to the bedrock.
The effect of subglacial water is thus twofold: it affects the sliding and the thermal regime, and in both these ways affects the ice flow. Subglacial water is therefore a crucial component in the dynamic evolution of ice sheets and needs to be incorporated into ice-sheet models. Here we investigate the effect of water at the ice base on sliding, deformation and the thermal regime.
Modelling the ice flow at a horizontal scale of the order of the ice thickness, and thus precise predictions of mass balance, requires the solution of the momentum-balance equations for all stresses (full-Stokes), since high gradients in geometry and the change of basal sliding properties arise on short length scales (Reference PattynPattyn and others, 2008). Studying the effect of subglacial water on the dynamics of ice streams requires coupling the hydrological system with the sliding properties, hence the velocity at the base of the ice sheet depends on water flux or water pressure. This simplest approach to subglacial water is incorporated into the full-Stokes model TIM-FD3 (Reference KleinerKleiner, 2010).
For the application we have chosen an area in which a variety of ice-sheet elements occur: large and medium ice streams, small glaciers, a mountain range with rock outcrops and stagnant flow. A suitable area with these characteristics, which is also reasonably surveyed, is located in western Dronning Maud Land, Antarctica, where two major ice streams drain into the adjacent Riiser-Larsen and Brunt Ice Shelves. An additional reason for choosing this area is that a previous modelling study (Reference Humbert, Kleiner, Mohrholz, Oelke, Greve and LangeHumbert and others, 2009) optimized the settings for the ice shelves, assuring us that we have a good represention of the buttressing effect of the ice shelves.
Model Description
The new three-dimensional (3-D) numerical ice-flow model, TIM-FD3, estimates the ice thickness and the 3-D flow regime (velocity, strain rate and stress fields), as well as the temperature distribution in response to environmental/boundary conditions. The model is implemented as a ‘multi-layer model’ (cf. Reference HindmarshHindmarsh, 2004), but with all stresses included in the momentum balance. Due to the reduced number of variables that need to be solved by linear sparse matrix solvers, the model is relatively fast and allows long integration times.
In the following discussions, x and y are the horizontal dimensions, z is the vertical dimension (positive upwards) and t is the time. The elevations of the upper and lower ice surface, zs = h (x, y, t) and zb = h (x, y, t),respectively, vary spatially as well as temporally, where H is the ice thickness. For this reason we use terrain-following coordinates, (ξ, ƞ, ζ, Ʈ) in which the transformation is implemented by
The new vertical coordinate, ζ, maps the local ice thickness onto unity such that ζ = 1 at the surface and ζ = 0 at the base. Although we change the coordinates, we still use vectors and tensors in terms of their components relative to the Cartesian base vectors. Partial derivatives according to the Cartesian coordinates must therefore be transformed to the new coordinates by applying the chain-rule
and similarly for ∂/∂y, ∂/∂z and ∂/∂t (cf. Reference Greve and BlatterGreve and Blatter, 2009, p.91–92). Many of the coordinate derivatives, ∂ξi,∂xj, vanish and
where all derivatives are affected by the vertical scaling. Note that, although ξ = x, ƞ = y and Ʈ = t, the entire spatio-temporal coordinate transformation needs to be considered.
Stokes flow
The flow of ice is modelled using the full-Stokes set of equations containing the conservation of mass and momentum. For an incompressible fluid they read in coordinate-free form as
and
where div u is the divergence of the velocity vector, u , with components, ui ∈ {u, v, w}, S is the Cauchy stress tensor, with components σij, i,j ∈ {x,y,z},and ρg is the volume force due to the density, ρ, and acceleration due to gravity, g . It is assumed that (1) body forces other than those arising from gravity can be neglected, due to the very slow flow and high viscosity, (2) density variations can be neglected, except in the gravitational term, ρ g , and (3) inertial (convective) terms can be ignored compared with gravity. The Cauchy stress tensor is given by
where pI and T are the isotropic and deviatoric parts of S , respectively, I is the identity tensor and tr (.) is the trace operator. The isotropic pressure, p, is in general not equal to the hydrostatic pressure as used in the shallow-ice approximation (SIA). Following Reference WhillansWhillans (1987) and Van der Veen and Whillans (1989a), the full-stress components from Eqn (9) can alternatively be partitioned into a resistive, R , and a cryostatic, L, component, so that
Where R is a tensor with components Rij, i,j ϵ {x, y, z} and
is the mean density of the overlying ice mass in the selected ice column. Both representations (Eqns (9) and (10)) of the Cauchy stress tensor are used in the text below. The deformation of polycrystalline ice is modelled using the commonly used Nye generalization (Reference NyeNye, 1957) of the Glen–Steinemann power-law rheology (Reference GlenGlen, 1955; Reference SteinemannSteinemann, 1954, 1958) given by
where
D
is the strain-rate tensor with components 
, i,j ϵ {x, y, z}, μ is the effective viscosity, n = 3 is the power-law exponent, A is the flow-rate factor, 
is the effective strain rate (second invariant of
D
) and E is the flowenhancement factor that parameterizes other physical contributions (e.g. impurities or damage). To avoid an infinite viscosity for vanishing effective strain rates we add a small 
to in Eqn (12). The flow-rate factor is parameterized as an Arrhenius relationship (Reference Paterson and BuddPaterson and Budd, 1982; Reference PatersonPaterson, 1994):
where A0 is a temperature-independent material constant, Q is the activation energy for creep, R is the universal gas constant, T is the temperature, T* is the homologous temperature (absolute temperature corrected for the dependence of the melting point on pressure) and γ is the Clausius– Clapeyron constant. The values for A 0, Q and R are taken from Reference PatersonPaterson (1994). If the ice reaches the pressure-melting point (T pmp), the rate factor in Eqn (13) is replaced by a rate factor depending also on the microscopic water content, W, following Reference DuvalDuval (1977) and Reference Lliboutry and DuvalLliboutry and Duval (1985):
as commonly used in numerical models (e.g. Reference Aschwanden, Bueler, Khroulev and BlatterAschwanden and others, 2012; Reference Sato and GreveSato and Greve, 2012). We do not solve for the water content explicitly, but assume a water content of W = 0: 01 (1%) for each gridpoint that is at the pressure-melting point. The chosen value is the upper limit of the validity range of Eqn (14) (Reference Lliboutry and DuvalLliboutry and Duval, 1985). In our simulations we use W ¼ 0 and W ¼ 0: 01, taken as the upper and lower limits of the water content. We further use the terms cold and temperate ice rheology in the following text and refer to Eqns (13) and (14) respectively.
By integrating the vertical component of Eqn (8) from the bottom to the ice upper surface and assuming a traction-free ice surface, the pressure at a certain position inside the ice can be found as (Reference Van der Veen and WhillansVan der Veen and Whillans, 1989b)
with vertical resistive stress,
accounting for bridging effects. The total isometric pressure, p, therefore consists of a dynamic part and a cryostatic part (last term in Eqn (15)). Based on the definition of the strainrate tensor,
each component of the tensor can be derived using general tensor analysis methods (e.g. Reference ArisAris, 1989) as
with ξi ϵ {ξ, ƞ, ζ,) xi ϵ (x, y, z} and ui ϵ (u, v, w} Together with the divergence of the stress tensor
where ei ϵ {ex, ey, ez } is the Cartesian basis, the horizontal components of the momentum balance (Eqn (8)) can be written in terms of the velocity field as
And similarly for the velocity component, v:
An expression for the vertical velocity, w, is obtained through the incompressibility condition (Eqn (7)), as
We also reformulate Eqn (16) as
Where the deviatoric stress components can be found from Eqns (12) and (18).
Heat transport
The general heat transfer by advection and diffusion with internal heat sources can be written as
where cp (T) and k (T) are the specific heat capacity and thermal conductivity of ice. Strain heating,
due to the deformation of ice, is the only heat source considered in the interior of the ice.
We further follow the analysis of Reference PatersonPaterson (1994) and neglect horizontal diffusion of heat, which is small compared with horizontal advection, and consider spatial variations in the thermal conductivity only vertically. Using the transformation rules for the partial derivatives (Eqns (3–6)), Eqn (24) can be written in coordinate-dependent form as
With
and
where w * is the effective vertical velocity, K is the vertically scaled diffusivity and S is the heat source term. The implicit dependency of k on temperature (Table 1) makes Eqn (26) nonlinear and reinforces downward advection due to the additional velocity term, 1/(ρcpH2 )∂k /∂ζ The remaining part,(ucx + vcy + ct ) /=H, accounts for the spatial–temporal coordinate transformation as a vertical grid velocity. The ice temperature is always kept below or equal to the pressure-melting point. Thus all heat that exceeds this limit is used for melting ice.
Used constants and model parameters
If the base is at the pressure-melting point and a basal layer of temperate ice exists, then the temperature gradient in the layer prevents the inflow of frictional and geothermal heat, so all that heat melts ice at the base as
where L is the specific latent heat of fusion, qgeo is the geothermal heat flux into the ice and q frict = u b τb is the frictional heat.
It is assumed that meltwater generated within a temperate basal ice layer drains away to the base immediately and contributes to the total basal mass balance according to
where ζ CTS is the cold transition surface (the interface between temperate and cold ice). Thus, a b = –(M b + M i) is the total basal mass balance (m ice eq. a–1, positive for freezing) at a given grounded ice position.
Evolution of the ice thickness
The evolution of the ice thickness is a direct consequence of the mass conservation expressed by the incompressibility, Eqn (7). The local horizontal ice flux that determines the evolution of the local ice thickness is obtained by integrating Eqn (7) from the bottom, z b, to the ice surface, z s,
where as = as (x, y, t) and ab = ab (x, y, t) are the surface and basal mass balances, respectively. They are given by the kinematic boundary conditions at the ice surface and base,
and
prescribing the motion of these free surfaces relative to the motion of ice, where u, v and w are taken at the surface and at the base, respectively. The signs of a s and a b are chosen to be negative for mass loss due to, for example, surface ablation or basal melting. At all lateral boundaries the ice thickness is prescribed. The calving front position is assumed to be constant in time, so no kinematic boundary condition needs to be given at the calving front.
Boundary conditions
Basal interface
At the ice base, z b = h(x, y, t) – H(x, y,t) or ζb = 0, shear stresses act in the tangential plane, and the normal stress must be compensated from the bedrock or water at the base
and
where 
is the mean density of the ice column
t
and
b
are the tangent and bi-normal unit vectors at the base. Although t and b could be arbitrary orthogonal vectors in the tangential plane at the lower surface, they are chosen to be aligned in the x-z plane and y-z plane, pointing into positive x- and y-directions, respectively. With this alignment the unit normal vector at the base is n = b × t and points outside the ice. The basal drag parameter, β2
(x, y, t) e.g. Reference MacAyeal, Bindschadler and ScambosMacAyeal and others, 1995), is an always-positive scalar field that in general varies in space and time and needs to be prescribed or parameterized by a sliding relation (see Eqn (56), below). The signs in Eqns (35) and (36) account for the orientation of
t
and
b
, so basal stresses oppose ice motion.
To obtain boundary conditions related to the horizontal velocity field we first use Eqn (37) to compute the pressure inside the ice,
and then eliminate the pressure term in Eqns (35) and (36), thus:
For the floating part of the base β2 = 0 (no basal traction). The basal vertical velocity component can be derived from the kinematic boundary condition (Eqn (34)) at the lower surface:
In this study we do not evolve the ice geometry, so ∂zb/∂t = 0 for all time-steps.
The temperature boundary condition at the grounded ice base depends on whether the ice is at the pressure-melting point or not. If the grounded ice base is below the pressure-melting point the following Neumann boundary condition applies:
The frictional heating, q frict, needs to be considered in this case only if basal sliding below the pressure-melting point is allowed. A second condition applies to basal ice at the pressure-melting point when the geothermal heat flux plus the heat flux generated by basal friction heating (sliding) is larger than the heat flux from the bed into the ice (Reference Payne and DongelmansPayne and Dongelmans, 1997),
and reads as the Dirichlet boundary condition. If the condition in Eqn (43) cannot be fulfilled, again Eqn (42) applies as the Neumann condition for the temperate ice base. The ‘if’ condition in Eqn (43) allows the temperature at the base to be less than the pressure-melting point during thermal evolution, even if the pressure-melting point was reached once. As a common approximation in ice-sheet and glacier models (Reference Payne and DongelmansPayne and Dongelmans, 1997; Reference PattynPattyn, 2003, 2010), the geothermal heat flux is assumed to be temporally constant and independent of heat transport in the bedrock, so changes in heat storage in the underlying bedrock cannot affect the basal heat budget in the ice model. We have further assumed that 
At the ice/ocean interface the ocean surface temperatures, T|b
= Tocn
(x, y, t) are prescribed as the Dirichlet boundary condition.
Upper interface
At the upper surface, zs = h(x, y, t) or ζs = 1 the tangential (wind) and normal (pressure) atmospheric stresses are assumed to be small compared with typical stresses in the ice and can therefore be neglected:
and
The tangential and bi-normal unit vectors, t and b , at the surface have a similar orientation to those at the base. In this alignment the surface normal unit vector is n = t × b pointing into the atmosphere. The boundary conditions for the horizontal velocity components can be obtained again by replacing the pressure (taken from Eqn (46)) in Eqns (44) and (45). They read as in Eqns (44) and (45), with β2 = 0 and cx, cy as well as u evaluated at the surface.
At the ice surface the mean annual air temperature, T |s = T atm (x, y t) is prescribed as the Dirichlet boundary condition.
Calving front
As for the upper and lower surface of the ice, a stress boundary condition is needed to solve the momentum equation (Eqn (8)). Since the convective transport of momentum across the calving front is very small, the traction is continuous at this boundary. Neglecting the atmospheric stresses, and assuming that dynamic stresses as well as viscous drag can be ignored, the pressure distribution in the ocean is hydrostatic. The dynamic boundary condition therefore reads as (Reference Weis, Greve and HutterWeis and others, 1999)
where p ocn is the hydrostatic pressure in the ocean and ρsw is the density of sea water. We assume that the vertical shear stresses, σ xz , σ yz , in the ice shelf far away from the grounding line can be neglected, as well as bridging effects (pressure is cryostatic), and approximate the velocities at the calving front using the vertically integrated shallow-shelf approximation (e.g. Reference MacAyeal, Rommelaere, Huybrechts, Hulbe, Determann and RitzMacAyeal and others, 1996).
Inflow
Along the lateral boundaries (lb) of the computational domain, Dirichlet boundary conditions for velocities are specified. We chose SIA velocities as (Reference Payne and DongelmansPayne and Dongel-mans, 1997)
where sliding at the base is neglected ( u b = 0). The horizontal heat transport is neglected (Reference PaynePayne, 1995; Reference Payne and DongelmansPayne and Dongelmans, 1997; Reference PattynPattyn, 2003) and the temperatures are prescribed.
Subglacial Hydrology Model
The temporal evolution of the water layer thickness, H w, is given by the continuity equation
where 
is the depth-averaged water velocity vector (two-dimensional). The subglacial water moves in the direction of decreasing hydraulic potential, Φ, a function of elevation potential and water pressure. Assuming that basal water pressure is equal to the overlying ice pressure, the hydraulic potential can be written as (Reference ShreveShreve, 1972)
where ρ
w is the density of water. The subglacial water is assumed to flow in a thin film of water (a few millimetres thick). In this thin-film model, water flow is assumed to be laminar; thus becomes (Reference WeertmanWeertman, 1966, 1972)
where μ w is the viscosity of water.
It is further assumed that the timescales for the water flow are much smaller than the timescales for the ice flow, so the water system reaches steady state (∂H
w
=∂t = 0) in every time-step of the ice-flow model. At steady state, the basal melting rate must balance the water flux divergence and the steady water flux, 
can be obtained by integrating the basal melt rate over the whole drainage area, starting at the hydraulic head in the direction of the hydraulic potential gradient. We use different grid-based flux-routing schemes that are widely used in the literature (e.g. Reference Le Brocq, Payne, Siegert and AlleyLe Brocq and others, 2009; Reference PattynPattyn, 2010; Reference Carter and FrickerCarter and Fricker, 2012) to perform this integration. The flux routing is based on the assumption that the outgoing water flux from an area is equal to the incoming flux plus the melt rate in the area. In contrast to the work of Reference Le Brocq, Payne, Siegert and AlleyLe Brocq and others (2009) and Reference Carter and FrickerCarter and Fricker (2012), where they used a prescribed basal melting rate, we compute the basal melt rate at every time-step. Hence, changes in the basal thermal regime, caused by changes in basal sliding conditions, feed back into the amount of meltwater available at the base.
We have implemented three different flux-routing schemes (Reference Quinn, Beven, Chevallier and PlanchonQuinn and others, 1991; Reference Budd and WarnerBudd and Warner, 1996; Reference TarbotonTarboton, 1997) in our flow model that differ in the number of neighbours and in the direction of flow. An overview of the main characteristics of the three schemes is given by Reference Le Brocq, Payne and SiegertLe Brocq and others (2006).
Basal Sliding
Basal sliding is described by a Weertman-type sliding law (Reference WeertmanWeertman, 1957, 1964):
Where
are the velocity and basal shear stress vectors in the basal plane and 
is the basal normal pressure, taken as the ice overburden pressure. The stress and pressure exponents are chosen as p = 3 and q = 2 (Reference GreveGreve, 2005). The basal sliding parameter, C
b, depends on, for example, bed material and roughness, and is modified to allow sliding below the local pressure-melting point, T
pmp (Reference Budd, Jenssen, Veen and OerlemansBudd and Jenssen, 1987; Reference Hindmarsh and Le MeurHindmarsh and Le Meur, 2001; Reference GreveGreve, 2005; Reference Pattyn, de Brabander and HuyghePattyn and others, 2005):
where the exponential decay, v (sub-melt sliding parameter), measures the temperature range in which a significant amount of sliding is allowed below the pressure-melting point. According to Weertman’s theory of sliding, the sliding velocity depends on the thickness of the basal water layer (Reference WeertmanWeertman, 1969; Reference Weertman and BirchfieldWeertman and Birchfield, 1982). We therefore extend the sliding law to account for the subglacial water layer thickness as
where 
is the typical scale of the layer thickness. The approach is similar to that of Reference Johnson and FastookJohnson and Fastook (2002), but here the increase of sliding velocity due to basal water is limited to a factor of 10, due to the exponential saturation function. We expect that if a sufficient water layer thickness is reached, more meltwater would not lead to additional sliding. In this study we take 
as a constant (Table 1), assuming the same subglacial material and roughness throughout the whole model domain.
In the limit of the SIA, Eqn (52) would read as Dirichlet boundary conditions for the horizontal velocity components, since the bedrock is assumed to be horizontal for each position at the base, thus 
In general the slope of the bedrock cannot be neglected. We therefore rewrite Eqn (52) as
that can be used together with Eqn (53) to replace the righthand sides in Eqns (35) and (36). Due to the nonlinear coupling between the velocity and the shear stress at the base in Eqn (56), the velocity field needs to be computed in an iterative way. The basal sliding parameter, β 2, varies over time, depending on the local subglacial water layer thickness and ice temperature, as well as basal shear and normal stresses.
Numerical Solution Scheme
Although we only consider a fixed geometry in this paper, we present the full solution scheme of the model, for the sake of completeness. The equations are discretized on a co-located grid (thus all variables are located at the same grid position) that is, in general, non-equidistant in the vertical and horizontal directions using a second-order finite-difference approximation (Reference Payne and DongelmansPayne and Dongelmans, 1997; Reference PattynPattyn, 2002). In the terrain-following coordinate system, the grid is given by N = Nξ Nƞ Nζ gridpoints connected through orthogonal lines.
Due to the nonlinear coupling of the horizontal velocity fields, a successive substitution of 
and 
is used. The finite-difference forms of the transformed differential equations (Eqns (20) and (21) together with their boundary conditions, Eqns (39) and (40) evaluated at the upper and lower surfaces) can be written in matrix-vector notation as
Where l is the iteration index, and U* = (U* 1,, U* N)T, V* = (V* 1,, V* N)T are the preliminary solutions based on the previous iteration. The new velocity vector is obtained by an under-relaxation scheme,
and analogously for V l + 1, where α = 0: 8 (Reference Pimentel, Flowers and SchoofPimentel and others, 2010) is the relaxation parameter. During each iteration, l, of the nonlinear iteration scheme we solve Eqns (57) and (58) using an iterative linear solver for large sparse matrices (ILU0-GMRES; Reference SaadSaad, 2003). Based on the new estimates for the horizontal velocity field we then compute the vertical velocity, w, in each ice column using Eqns (22) and (41) in their discretized forms and a direct band diagonal solver (Reference Press, Teukolsky, Vetterling and FlanneryPress and others, 1992). A very similar approach is used for Rzz (Eqn (23)). Convergence of the linear solver is assumed if
and analogously for V*. We stop the nonlinear iteration process if
or
At this point the Stokes equations are solved for a given temperature and geometry.
The advective terms, u∂T /∂ξ, v∂T /∂ƞ, and w* ∂T /∂ζ, in Eqn (26) are discretized using the hybrid difference scheme (HDS; Reference SpaldingSpalding, 1972) which switches the discretization between the second-order central difference scheme (CDS) and the first-order upwind difference scheme (UDS) according to the local cell Peclet number. Thus 
where ℒHDS is the one-dimensional spatial discretization operator of the advective terms and s is a weighting factor depending on the Peclet number, Pe, as
The HDS is applied to all spatial dimensions and avoids numerical instabilities arising from the CDS for high velocities (Pe > 2); it also limits the amount of artificial diffusion due to UDS for low velocities (small Pe). The conductive term, K∂2T /∂ζ2, in Eqn (26) is discretized using a second-order central difference scheme for the second derivative, denoted here as CDS2. Expressions for the discretization schemes (CDS, UDS and CDS2) on non-equidistant grids can be found in the literature (e.g. Ferziger and Perić, 2002). The time-stepping is performed using the semi-implicit Crank–Nicolson scheme, where the temperature at the next time-step depends on the temperature, the source term (Eqn (29)) and the material properties, cp and k, at the actual time-step. We use a constant time-step of 1 year and an asynchronous coupling, where the temperature module is run repeatedly for every time-step coupled to the ice-dynamic module (Stokes) every fifth step (1: 5 coupling). The Crank–Nicolson scheme allows much larger time-steps, but we found an effective time-step of 1 year to be a good compromise between temporal resolution and computational time. The smaller time-steps in the temperature module can further be seen as a sub-iteration scheme to resolve the temperature dependency of the thermal diffusivity and heat capacity.
The hydrological component of the numerical model is coupled to the temperature at the end of each time-step. Basal ice temperatures, the resulting basal melt rates and the ice geometry from the flow model are used as input for the hydrology module, which computes the steady-state water-layer thickness as an input quantity for the ice-dynamic module.
Although ice-thickness evolution is already implemented in TIM-FD3, and has been tested in the Ice Sheet Model Intercomparison Project for Higher-Order Models (ISMIPHOM; Reference PattynPattyn and others, 2008) and also in the Marine Ice-Sheet Model Intercomparison Project (MISMIP3d; Reference PattynPattyn and others, 2013), all simulations performed in this study are based on a fixed geometry. Several studies (e.g. Reference Durand, Gagliardini, Zwinger, Le Meur and HindmarshDurand and others, 2009; Reference Gladstone, Lee, Vieli and PayneGladstone and others, 2010; Reference PattynPattyn and others, 2012), based on idealized geometries, have shown the importance of a very fine horizontal resolution of the grounding zone (transition between grounded and floating ice) to capture the ice evolution well in these areas. Running the model in the geometry evolution mode would therefore require very high-resolution data of the bedrock topography under the ice and, even more difficult to observe, the bedrock topography under the floating ice shelves.
Study Area
The study area is located at the margin of western Dronning Maud Land and the eastern part of Oates Land, Antarctica, at 71–78.58° S, 2–30° W (Fig. 1). The region consists of the Brunt Ice Shelf (BIS) in the southwest, the Riiser-Larsen Ice Shelf (RLIS) in the northeast and their drainage areas. The BIS and the RLIS are separated by the Stancomb-Wills Ice Tongue (SWIT), the floating extension of the Stancomb-Wills Ice Stream (SWIS) that crosses the Caird Coast (CC).
Map of the Brunt Ice Shelf, the Riiser-Larsen Ice Shelf and the adjacent grounded ice parts of western Dronning Maud Land, including feature names and the outline of the model domain (dashed line). Underlain is a Mosaic of Antarctica (MOA) image (Reference Haran, Bohlander, Scambos and FahnestockHaran and others, 2005) acquired between late 2003 and early 2004. The red arrows in the lower part of the figure indicate the origin and orientation of the local coordinate system used in this study. Different glaciological units are marked with coloured lines. The inset shows the location of the area at the eastern margin of the Weddell Sea, Antarctica.
The BIS and SWIT areas are highly heterogeneous: between the relatively thick and fast SWIT and the considerably thinner and slower BIS, an area exists that is made up of a discontinuous mass of ice blocks that have calved at the grounding line and are frozen together by a melange of sea ice, small icebergs and accumulated snow (Reference Hulbe, Johnston, Joughin and ScambosHulbe and others, 2005; Reference Humbert, Kleiner, Mohrholz, Oelke, Greve and LangeHumbert and others, 2009). Two major rift systems evolve under the stress regimes between the SWIT and the RLIS (rift 1) and between the SWIT and the BIS (rift 2). In the rift areas the ice is thinner and formed of a mechanically weaker mixture of meteoric and mainly marine ice. The margin of the BIS is pinned by the relatively small McDonald Ice Rise (Reference Humbert, Kleiner, Mohrholz, Oelke, Greve and LangeHumbert and others, 2009).
The RLIS area is relatively homogeneous compared with the BIS/SWIT area and contains three major ice rises near the calving front (Lyddan Island, Kvitkuven and Skjoldet). The ice shelf in this area is further stabilized by several small pinning points inside the ice shelf, as well as near the calving front. Further important features are the ice streams Plogbreen and Veststraumen, that are responsible for most of the mass discharge over the grounding line in this region. A grounded area (ice rise/ice rumple) is located in Vest-straumen at about 74° S, 17.5° W and diverts the flow as can be seen on satellite images (Fig. 1). The two ice streams are separated by the dome-shaped Högisen.
Model set-up
The study area is relatively well covered by satellite altimetry measurements of surface elevation from the Ice, Cloud and land Elevation Satellite (ICESat)/Geoscience Laser Altimeter System (GLAS) (Reference ZwallyZwally and others, 2005). Data from 2003 were used, which cover approximately the same time period as the Mosaic of Antarctica (MOA) image (November 2003 and February 2004; Reference Haran, Bohlander, Scambos and FahnestockHaran and others, 2005). Ellipsoid heights were converted to freeboard heights using the World Geodetic System 1984 (WGS84) version of the Earth Gravitational Model EGM2008 (Reference Pavlis, Holmes, Kenyon and FactorPavlis and others, 2008). Ice-thickness datasets from expeditions, which are collected in the BEDMAP database (Reference Lythe and VaughanLythe and others, 2001), have to be assembled. We included additional airborne electromagnetic radar (EMR) measurements of the ice thickness, taken by the Alfred Wegener Institute between 1994 and 2008 (Reference Steinhage, Nixdorf, Meyer and MillerSteinhage and others (2001) give details of the method) as well as electromagnetic radiation (EMR) measurements provided by the British Antarctic Survey (personal communication from D.G. Vaughan, 2009). Older Russian EMR measurements taken in 1988/89 are used only for the grounded part of the ice and only in areas where no other data are available.
Line segments of the calving front and grounding-line position were manually digitized by visual inspection of the MOA image and existing datasets (Antarctic Digital Database (ADD), MOA) as well as surface elevations and TerraSAR-X images. Small pinning points with an area less than the grid size were manually increased to be represented by at least one gridpoint. In the freely floating part (5 km away from the grounding line) surface elevations from ICESat/GLAS, together with the floating conditions, are used to compute the ice thickness. This was necessary because only a few ice-thickness profiles exist in the floating areas. The required distribution of (x, y) was computed from locations where surface and ice-thickness measurements were available.
Previous model studies of the BIS–SWIT system have shown that the effect of increased damage, and therefore reduced stiffness in the rift areas, needs to be incorporated in simulations to better approximate the observed flow field (Reference Hulbe, Johnston, Joughin and ScambosHulbe and others, 2005; Reference Humbert, Kleiner, Mohrholz, Oelke, Greve and LangeHumbert and others, 2009; Reference Khazendar, Rignot and LarourKhazendar and others, 2009). We therefore reduce the associated stress-enhancement factors, E s = E-1/n , to E s = 0: 3 at rift 1 and E s = 0: 1 at rift 2. The positions of the two major rifts are the same as given by Reference Humbert, Kleiner, Mohrholz, Oelke, Greve and LangeHumbert and others (2009). Areas of ice melange located between the SWIT and the RLIS and between the SWIT and the BIS are assumed to be less stiff (E s = 0: 5) than the meteoric ice areas. For all other areas (grounded and floating) we use a constant enhancement factor of E s = 1: 0.
Annual mean temperatures described by Comiso (2000a), were used to provide the surface temperature field together with the surface mass balance (Reference Van de Berg, Van den Broeke, Reijmer and Van MeijgaardVan de Berg and others, 2006) and are applied at the ice upper surface. At the base the basal geothermal heat flux of Reference Shapiro and RitzwollerShapiro and Ritzwoller (2004) is used, unless stated otherwise. These boundary conditions were taken from the ALBMAPv1 dataset (Reference Le Brocq, Payne and VieliLe Brocq and others, 2010) and interpolated to the used grid positions.
Under the assumption of a constant geometry (∂H/∂t = 0), the basal mass balance, a b, under the ice shelf is derived from the ice-thickness evolution (Eqn (32)) and accounts for the divergence of the horizontal ice flux, and the surface mass balance, a s. For the free-floating ice the vertical mean velocity equals the surface velocity, since no vertical shearing occurs, so the basal mass balance can be estimated from observed velocities (Reference HumbertHumbert, 2010; Reference Neckel, Drews, Rack and SteinhageNeckel and others, 2012). Surface velocities derived from interferometric synthetic aperture radar (InSAR) (Reference Rignot, Mouginot and ScheuchlRignot and others, 2011) are used for this purpose (Fig. 2). Short distance variations arising from strong ice-thickness gradients (e.g. ice melange area) or velocity gradients between fast- and slow-flowing units (shear margins) were spatially averaged, in order to obtain a smooth distribution of a b.
Surface velocities derived from satellite radar interferometry, mainly from Advanced Land Observing Satellite (ALOS) Phased Array-type L-band synthetic aperture radar (PALSAR) measured in 2007, 2008 and 2009 (Reference Rignot, Mouginot and ScheuchlRignot and others, 2011).
In general, the thermodynamic conditions at the ice/ocean interface are very complex and depend on, for example, currents and vertical mixing in the ocean. Boundary conditions at the ice-shelf base should therefore be derived from an ocean model or observations. An adequate method to quantify surface fluxes and subsequent ocean melting rates at the ice/ocean interface (Reference Holland and JenkinsHolland and Jenkins, 1999) would require a realistic simulation of the sub-ice-shelf ocean circulation, based on a precisely known geometry and temporal evolution of the ice-shelf cavity, to accurately capture the basal melt rate under the ice shelf (Reference Schodlok, Menemenlis, Rignot and StudingerSchodlok and others, 2012). Our steady-state estimate of a b shows highest melt rates close to the grounding line at the SWIT (11 m w.e. a–1) and Veststraumen (5 m w.e. a–1) and close to the calving front of the SWIT (5 m w.e. a–1). Refreezing rates up to 4 m w.e. a–1 occur mainly in the small areas between individual icebergs in the melange area between the BIS and the SWIT. The estimated total average is a b ¼ 0: 19 ± 1: 1 m w.e. a–1 (net melting) and is significantly lower than the simulated melt rates of Thoma and others (2006; 0.88 m a–1) or Timmermann and others (2012; 0.94 m a–1).
The final model domain consists of Nξ = 281 by Nƞ = 213 horizontal nodes forming a 2.5 km × 2.5 km grid, and Nζ = 21 unequally spaced vertical layers. Model parameters and the used constants are summarized in Table 1.
Results
Spin-up and cyclic behavior
The model was run for 200 ka without basal sliding, under the present-day climate, to achieve a steady-state temperature and velocity distribution. Initial temperatures for the spin-up were chosen to lead to a cold ice base in the whole domain. To illustrate the thermal evolution during the model spin-up (SUt), the mean ice temperature and the temperate ice area fraction (TIAF) are shown in Figure 3 (red and black curves). Approximately 28% of the grounded area is at the pressure-melting point, and a mean ice temperature of –19°C is reached at the end of the spin-up. After ~100 ka, variations in temperature and TIAF become small, but never vanish. The observed period is ~2 ka, with varying amplitudes in a small band of ± 0.1°C and ±0: 5% TIAF. If we compare the thermal evolution of this spin-up with a similar simulation (SUc), where we used the cold ice rheology everywhere, the oscillations of the ice temperature and also of the TIAF vanish after ~100 ka (pale red and grey curves in Fig. 3). The oscillations in SUt are therefore related to feedbacks between the thermal and the dynamical subsystems and need to be considered in the later analysis. The mean temperature and TIAF at the end of SUc are slightly higher than in SUt.
Temporal evolution of the mean pressure-corrected ice temperature (red) and the temperate ice area fraction at the base (TIAF, black) during the 200 ka spin-up using 1% water content where the ice is at pressure-melting point (SUt) and the corresponding quantities in pale red and grey for the cold ice rheology applied everywhere (SUc).
The distributions of the basal homologous temperature and the resulting basal melt rates at the end of the spin-up, SUt, are shown in Figure 4 (upper row).
The homologous temperature at the base (upper left) and basal melt rate (upper right) at the end of the 200 ka spin-up simulation, SUt, are used the compute the basal water flux (lower left) and the basal water layer thickness (lower right). Here we have used the Reference Budd and WarnerBudd and Warner (1996) flux-routing scheme, BWt, as an example. Contours of the hydropotential are overlaid as thick black curves every 5 MPa. In the upper left panel the –0.1°C contour is shown in light grey.
Temperate areas are mainly located close to the grounding line and in the ice-stream areas. These regions are subject to basal melt rates of a few centimetres per year, with largest values (>10 cm a–1) close to the grounding line of the SWIS. At the outer margins of the temperate ice areas some grid nodes also show very low freeze-on rates. Temperate ice areas with significant melt generally coincide with areas of ice-stream flow. The simulated cold area in the SWIS at x ≈ 270 km, y ≈ 200 km is located where observed surface velocities (Fig. 2) are very low. In our simulations the location of this area, hereinafter referred to as a ‘sticky spot’, varies over time and is, together with an area near the grounding line of CC, predominantly responsible for the fluctuations shown in Figure 3. To eliminate the effect of small temporal variations and to investigate the mean behaviour over time, we compare model results as the temporal mean of the last 10 ka, unless stated otherwise. We further vary the basal geothermal heat flux, the stress-enhancement factor and the inflow velocity at the domain lateral boundary independently, to illustrate the model sensitivity to the standard parameters used in this study. The results are summarized in Table 2.
Comparison of the temperate ice area fraction (TIAF) at the base, the area mean basal homologous temperature, T * b, the volume mean velocity magnitude, |u |, and the mean surface velocity difference, Δu s = u sim – u obs, if uobs ≥ 5 m a–1, and the correlation, r, between the observed, u obs, and simulated, u sim, surface velocities for different spin-up simulations without basal sliding. (All quantities except the mean velocity are computed for the grounded ice part only)
We find that the model results are not sensitive to the inflow velocities in the range of zero velocity up to doubled SIA velocities. The lateral boundaries mostly affect gridpoints close to the boundary, since the velocities in that area (far away from the grounding line) are generally low. Since the horizontal advective heat transport for low velocities is not very efficient, we assume a low sensitivity to the prescribed temperatures at the margin. The area of temperate ice at the base is very sensitive to the applied basal heat flux and ranges between 15% and 56% for the 40 and 80 mW m-2 simulation, respectively. The TIAF of the standard spin-up is below the value given by the 60 mW m-2 simulation. This value corresponds well with the mean geothermal heat flux of 56 ± 4 mW m-2 in the model domain for the Reference Shapiro and RitzwollerShapiro and Ritzwoller (2004) dataset that we apply in this study. Mostly the mountainous areas do not reach the pressure-melting point in the 80 mW m-2 simulation because of the very low surface temperatures and the relatively thin ice cover.
Differences in the mean temperatures (not shown here) between all spin-up simulations are small. They vary between –18.9°C (SUc) and –19.6°C (80 mW m-2) and are therefore only useful to characterize the transient behaviour during the spin-up, but not to compare simulation results. This is also found for the mean velocity that varies between 35.9 m a–1 (SUc) and 46.5 m a-1 (80 mW m-2) with a value of 41.4 m a-1 for the standard spin-up.
The simulation performed with a constant enhancement factor, E = 1, for the whole modelling domain shows only a minor reduction in the mean ice velocity. This is due to a higher buttressing effect of the SWIT, since no shearing or decoupling at the two major rifts is considered in this model run. The mean surface velocity differences, Δu s = u sim – u obs, where u obs are the velocities given by Reference Rignot, Mouginot and ScheuchlRignot and others (2011) (Fig. 2), are between –5 and +8 m a–1 and have standard deviation (stddev) of ~20 m a–1. The correlation between observed, u obs, and simulated, u sim, surface velocities in the grounded ice area is ~0.8 for all simulations.
Reference simulation
Starting from the standard spin-up simulation, the reference simulation (REFt) is integrated for another 30 ka, where we allow for basal sliding computed with the Weertman-type sliding law without basal water (f w = 0 in Eqn (55)). In addition to REFt we continue the spin-up without changes in the set-up in the control simulation (CONt) to distinguish between sliding- and deformation-related processes.
A direct comparison of the simulated surface velocities, u sim, and observed velocities, u obs, is shown in Figure 5 for all grounded and floating gridpoints. The simulated velocities are slightly higher than the observed velocities. The mean difference is 19.4 m a–1 (±78.7 m a–1 stddev) and the overall correlation is very high (r = 0: 98). We compare the velocities only where the observed velocities are ≥5 m a–1, since the very low observed velocities show a large spread in magnitude and direction (cf. Fig. 2). If we compare only the grounded points, the number of observations is reduced from 34 366 to 22 103 gridpoints and the mean difference reduces to 2.35m a–1 (±22.2 m a–1 stddev). The correlation between observed and simulated velocities reduces to r = 0: 83 but is still high (Table 3).
Scatter plot of the simulated, u sim, and observed, u obs, surface velocities separated for the floating (grey) and the grounded (black) gridpoints. The inset shows the histogram of the differences between observed and simulated velocities using bins of 25 . The statistical quantities quoted in the figure apply to the grounded and floating part of the ice. Simulated velocities are averaged over the last 10 ka of model integration.
The mean surface velocity difference (Δu s = u sim – u obs, if u obs ≥ 5 m a–1) (m a–1), the standard deviation, stddev (m a–1), and the correlation, r, between the observed, u obs, and simulated, u sim, surface velocities are given for all sliding simulations. Quantities are calculated for the grounded ice only (gr.) and for the entire domain (all) averaged over the last 10 ka model integration
The simulated basal homologous temperatures, the basal sliding parameter, β
2, the basal sliding velocities and the surface velocities are shown in the left column of Figure 6, from top to bottom. As the sticky spot moves along the same path again and again it is not visible as such in the 10 ka averages. The modelled sliding velocities are relatively low because of the low sliding parameter, . Highest sliding velocities (~25 m a–1) occur close to the grounding line at CC, and in the southwestern part of Veststraumen, where geometry gradients and thus basal shear stresses reach their maximum. In general, sliding in the model occurs where the known ice streams are located. Only a few areas close to the grounding line are not sliding due to the low basal temperature (e.g. at Högisen, Mannefallknausane, and northeast of Plogbreen). High values for the basal sliding parameter, β
2 (dark blue in Fig. 6), represent areas of high basal resistance, controlled by the basal temperatures, that do not lead to significant sliding velocities even in areas of sub-melt sliding (mostly in light blue). Along the grounding line of CC, where surface gradients and therefore basal shear stresses are very large, β
2 is low. Although the same bed material and roughness, , as well as the same sliding mechanism (p, q) is assumed for the entire model domain, the sliding parameter, β
2, varies spatially over several magnitudes.
Comparison of the reference, REFt, simulation (left column) with sliding including Reference Budd and WarnerBudd and Warner (1996), BWt, flux routing (right column) as temporal mean of the last 10 ka at the end of the sliding experiment. From top to bottom: basal homologous temperature; sliding parameter, β 2; basal velocity (if > 10-3 m a-1) and surface velocity. The velocity is represented as magnitude (colour-coded) and flow direction (arrows). The –0.1°C temperature contour is shown in light grey in the uppermost panels. In the lowermost panels the 200 m a–1 contour is drawn as a thick black curve on top of the other contours.
The highest modelled surface velocities of the grounded ice can be found close to the grounding lines of the SWIS (~500 m a–1) and both branches of Veststraumen (~250 m a–1). Compared with the high jump in velocity at the grounding line of the SWIS, the transition between the grounded and floating parts of Veststraumen is very smooth. Surface velocities at the calving front of the SWIT are slightly less than 1500 m a–1 and much higher than the surface velocities at the calving front of Veststraumen (~600 m a–1). The simulated velocities in the RLIS are higher than the observed velocities, especially close to the calving front. The simulated velocities in the BIS area are lower than the observed velocities. The differences between simulated and observed velocities that arise in the scatter plot (Fig. 5) for u obs ≤ 300 m a–1 and 300 m a-1 ≤ u sim ≤ 750 m a–1 are located mainly in the melange area close to rift 1 between the SWITand the RLIS and near the calving front of the RLIS. With respect to the statistical quantities and the overall flow pattern, the simulated and observed velocities agree fairly well.
Reference simulation REFc, where we apply the cold ice rheology everywhere, leads to slightly different results. The mean difference from observed surface velocities decreases to 5.4 m a–1 (± 71.2 m a–1 stddev) for the entire domain, while the absolute mean difference for the grounded part increases. In REFc the simulated surface velocity is still below the observed surface velocity and has a mean difference of –4.4 m a–1 (±19.7 m a–1 stddev).
Basal sliding including subglacial water
In this subsection, a set of sensitivity studies investigating the influence of different basal water routing schemes on the modelled ice flow is described. The set-ups are comparable with the reference simulation, but with Budd and Warner (BWt), Quinn (QUt) and Tarboton (TAt) flux routing incorporated in the sliding law (f w ≥ 0 in Eqn (55)).
Starting from the subglacial water flux and the derived water layer thickness at the end of the spin-up (lower row in Fig. 4), both quantities are recomputed every time-step, as well as the sliding parameter,β 2.
The basal homologous temperatures, basal sliding parameters, β 2, basal sliding velocities and surface velocities for an experiment with BWt routing are shown from top to bottom in the right column of Figure 6. Compared with the reference simulation (REFt), shown in Figure 6 (left column), the BWt set-up results in lower temperatures in the area of the sticky spot and close to the grounding line at CC. In most areas β 2 is reduced, due to lubrication by subglacial water, but along CC the reduced temperatures lead to larger areas of increased2. Especially at CC the areas of enhanced basal motion coincide well with flow features visible on the satellite background image and show individual outlet glaciers in more detail (basal velocities in Fig. 6). Although the TIAF for simulations including subglacial water decreases compared with the reference simulation, the mean velocity increases for the cold and temperate ice rheology. The differences in the surface velocity are not as large as one would expect from the differences in the basal velocity. This is due to a reduction of the deformational part of the velocity in the BWt sliding simulation, caused by lower ice temperatures and therefore stiffer ice than for the REFt simulation.
The differences between simulated and observed surface velocities for the grounded ice increase for the shown simulations, REFt and BWt, whereas they decrease for the cold ice rheology simulations, REFc and BWc (Table 3).
To illustrate the temporal evolution of the basal thermal regime, maps of the homologous temperatures are shown at selected time-steps (223–226 ka) for the BWt simulation in Figure 7. The main path of the sticky spot is overlaid as a thick dashed line. At 223 ka the sticky spot is not visible in the temperature field at the observed (low surface velocity in Fig. 2) position in the SWIS, but further upstream of the glacier. In subsequent time-steps this cold area first increases and then moves down to the sticky-spot location (224– 225 ka). Later the cold area decreases (226–227.5 ka), another cold ice area starts to develop at the same upstream location (227.5–228 ka; not shown) and the cycle begins again. Similar behaviour can be observed for a colder ice area at CC close to the grounding line, while in all other areas the temperature remains mainly unchanged. All simulations including the spin-up, where we use the temperate ice rheology (SUt, CONt, BWt, QUt, TAt) show this cyclic behaviour.
Maps of the simulated basal homologous temperature of the grounded ice area at the selected time-steps, 223–226 ka, for the BWt simulation. The thick dashed line is the main path of the sticky-spot location over time. The –0.1°C temperature contour is shown in light grey.
We take the TIAF at the base as a proxy for the basal thermal regime and also as a measure of the importance of basal sliding, since sliding is only allowed if the temperature at the base is at pressure-melting point or slightly below (sub-melt sliding). The TIAF, the mean ice velocity (averaged by volume) and the mean basal homologous temperature are shown in Figure 8 for the total time span of 30 ka of the sliding and the control simulations. We discriminate between simulations with temperate (CONt, REFt, BWt, QUt, TAt) and entirely cold (CONc, REFc, BWc, QUc, TAc) ice rheology.
Temporal evolution of the mean basal homologous temperature of the grounded ice, the mean over all ice velocities and the temperate ice area fraction (TIAF) at the base (top to bottom) during the 30 ka for the control (CONt, black), reference (REFt, red) and sliding simulations (dark to light blue) for the different flux-routing methods (BWt, QUt, TAt). Simulations performed with the cold ice rheology (CONc, REFc, BWc, QUc, TAc) are represented by thinner curves for comparison.
Compared with the reference (REFt) and the control (CONt) simulations, all three sliding simulations under wet basal conditions (BWt, QUt and TAt) have a lower TIAF, indicating a colder ice base. Although the sliding area is reduced, the mean ice velocity increases compared with the reference and control simulations.
All temperate ice rheology simulations show strong temporal variations, due to the internal feedback between flow and temperature fields resulting from the microscopic water content, as already seen in the spin-up experiment. The velocities adjust to the new basal conditions in <2000 years, followed by a slight increase over the simulation time. At 222.5 ka the minimal TIAF is reached in the REFt. At that time, the simulated mean velocity without sliding (CONt) is higher than the velocity with sliding (REFt). Up to 213 ka the mean velocity and the TIAF of all experiments are nearly in phase, but have different amplitudes. After that BWt, QUt and TAt are running out of phase, while REFt and CONt are slightly in phase for another ~7 ka before they diverge. Simulations with basal water exhibit larger temporal variability than the REFt and CONt experiments. At least for the first half of the sliding simulations, the mean velocity of the reference and control simulations shows a saw-tooth-like waveform with a repetition frequency of ~5000 a–1 that is modified by another frequency close to ~2000 a–1. For all sliding simulations, the occurrence of the maxima in TIAF corresponds well with maxima in the mean velocity. In the non-sliding control simulation the maxima in the mean velocity occur after the maxima of the TIAF (~300–500 years, e.g. at 214, 218 and 222 ka). The three flux-routing schemes work similarly up to ~226 ka, when they start to diverge. Simulations performed with a cold ice rheology (thinner curves in Fig. 8) result in higher TIAFs but lower mean ice velocities, and show much smaller temporal variations. Although the temporal variations for all quantities shown in the figure are small in the cold ice rheology simulations, they increase considerably for simulations that include basal water in the sliding relation. Simulations with cold and temperate ice rheology show the same general picture: more basal sliding results in higher mean velocities, but smaller TIAF.
Sticky-spot area
In the observed surface velocity field (Fig. 2), as well as in the simulations performed in this study, a relatively slow-moving part in the SWIS has been identified at x ≈ 270 km and y ≈ 200 km. This area shows a large temporal variability in all temperate ice rheology simulations, including the spin-up experiment. In Figure 9 the surface velocity, the basal velocity, the basal friction heat, the basal strain heat, the thickness of the basal temperate ice layer and the basal homologous temperature are shown for CONt, REFt and BWt. For a later comparison CONc is also shown.
Temporal evolution of the ice at the sticky-spot location for the control (CONt, black), the reference (REFt, red) and the sliding experiment including Budd and Warner routing (BWt, blue). The surface velocity, u s, the basal velocity, u b, the basal friction heating, q frict, the strain heat of the lowermost gridpoint, the thickness of the basal temperate layer, H temp, and the basal homologous temperature are shown from top to bottom. For comparison, the values of the cold ice rheology simulation without basal sliding (CONc) are shown in grey.
In the figure, two main classes of simulation need to be distinguished: (1) experiments where basal sliding is suppressed (CONt and CONc) and (2) simulations where basal sliding is allowed (REFt and BWt). Since all sliding experiments that include basal water routing have very similar behaviour, only the BWt simulation is shown in the figure, as a representative. The simulation with cold ice rheology (CONc) shows no temporal variations in any quantity (e.g. a constant surface velocity of 19 m a–1), meaning that the ice is already in steady state after SUc. All other simulations in Figure 9 show strong temporal variations. The temporal behaviour is similar for at least 15 ka after the spin-up, before the simulations start to deviate considerably. As shown above (Fig. 8), the minimum basal temperature decreases with increasing basal sliding (CONt: –7.3°C; REFt: –7.8°C; BWt: –8.1°C). The periods of temperature cycles are slightly less than 5 ka and much longer than the periods identified for the domain average values (Fig. 8). The minima in the basal temperature are missing for the REFt and CONt simulations, if a temperate basal ice layer is available. Both simulations show a maximum thickness of temperate ice of H temp = 9.6 m, which is the thickness of the lowermost ice layer at this position. The missing minima in the strain-heating term are also correlated with the occurrence of temperate ice at the base. The range of strain heating over time is between ~10-4 and 1: 5 × 10-3 W m–3. The temporal means (CONt: 0.8; REFt: 0.9; BWt: 0: 7 × 10-3 W m–3) are clearly above the 0: 5 × 10-3 W m–3 value for the cold ice simulation, CONc. The heat supplied by basal friction is small, with temporal means of 1 × 10-3 W m–2 (REFt) and 4 × 10-3 W m–2 (BWt) compared with 53 × 10–3 W m–2 supplied by the geothermal heat flux. The sliding velocity is very low at this position (REFt: 0.7 m a–1, maximum) even for the enhanced sliding simulation (BWt: 4.4 m a–1, maximum). Highest surface velocities are reached in peaks in the REFt simulation, where the velocity matches the observation of 33 m a–1. In the more general view, the surface velocities in the BWt simulation show higher local maxima, with magnitudes of ~30 m a–1, which are in the range of 10% of the observations. During time intervals with a cold ice base, the simulated surface velocities are much smaller than observed.
Discussion
We run the model into thermal steady state based on the present-day climate (surface temperatures). This overestimates the mean ice temperature and therefore the contribution of ice deformation to the surface velocity. However, it is only possible to study the impact of subglacial water and microscopic water inclusions when perturbation experiments are started from a model in a steady state, so the results are not overlain by the model transient behaviour.
Due to the multilayer approach of the full-Stokes model we were able to run the simulations over long timescales. At the end of the thermocoupled spin-up simulation, performed over 200 ka without basal sliding, the simulated surface velocities already agree well with the observed velocities. We therefore conclude that the occurrence of ice streams in the model domain is mainly controlled by the geometry (Reference Winsborrow, Clark and StokesWinsborrow and others, 2010).
The basal sliding parameter, , used in subsequent simulations over another 30 ka is relatively low, taking into account the already high contribution of deformation to the total velocity. We do not state that our choice of sliding parameters (p, q, v) is the only way to match the observations. However, the aim of this study is to analyse the effect of subglacial water as an additional contributor to the sliding, and hypothesize a similar effect for other parameter choices.
Including the basal water in the sliding relation increases the mean velocity of the ice and therefore ice transport over the grounding line, compared with the reference simulations. The three different flux-routing methods lead to similar results, most likely because the main water flux points along the grid orientation in the y –direction. In this case the differences between the different schemes are smallest (Fig. 4; Reference Le Brocq, Payne, Siegert and AlleyLe Brocq and others, 2009).
The mean differences between all sliding simulations and the observed velocity field are <5 m a–1 for the grounded ice, with a high correlation, ~0.8. The best fit to the observations is given by the Budd and Warner routing in combination with a cold ice rheology, but we hypothesize that other combinations of model parameters in the sliding law or a lower value of the microscopic water content would lead to slightly different results.
Simulations with basal water do not show increased sliding at all locations compared with the reference simulations. Some areas at CC show a slight decrease in the surface velocity because of the changes in the basal thermal regime. This corresponds well with pronounced glacier flow in that area. The simulations show that sliding, especially the enhanced sliding due to basal water, leads to a reduction of the TIAF and the mean basal temperature (Fig. 8). Since both quantities are taken as domain averages, temporal changes in small areas or just spatial differences between different time slices are not resolved well, but are clearly visible between different time-steps (Fig. 7). Although the area of sliding decreases and the viscosity of the ice at the base increases due to lower temperatures, the mean ice velocity increases with increased sliding. We conclude that enhanced flow occurs in smaller areas. The enhanced flow is also visible in the distribution of the basal velocities in Figure 6.
Our simulations also show decreasing basal temperatures with increasing basal sliding. The effect is largest if we use enhanced sliding due to basal water (lower basal temperatures and TIAFs in Fig. 8). In a two-dimensional (2-D) flowline model, Reference PattynPattyn (1996) argued that an increase in basal velocity caused an increase in vertical advection rates transporting colder ice downward to the base, decreasing the area subjected to basal melting (TIAF), decreasing the total area of lubrication and therefore also decreasing the area of enhanced sliding and the velocity itself. It was found that whenever basal velocity plays a major role, its effect on the ice-sheet dynamics is largely controlled by the ice-sheet temperature. It was further noted that the disintegration of an ice sheet is partly counteracted by colder basal temperatures due to enhanced sliding. In contrast to the findings of Reference PattynPattyn (1996) the overall mean velocity increases if we allow for sliding.
In our 3-D simulations the flow is more concentrated in smaller areas (the flow pattern changes horizontally), which was not possible in Pattyn’s 2-D flowline model. The simulations BWc, QUc and TAc are very similar to the type 1 sliding relation used by Reference PattynPattyn (1996), but show only minor temporal variations that are not comparable with the cyclic behaviour found by Pattyn, because we use a fixed geometry. If we account for basal water in the sliding, the sliding velocities increase and therefore the advection of colder ice to the base increases, explaining the larger temporal variations of the TIAF, the mean basal temperature and the mean velocity in the BWc, QUc and TAc simulations, compared with the REFc simulation.
Reference Van Pelt and OerlemansVan Pelt and Oerlemans (2012) recently described a basal water feedback based on the assumption that higher sliding velocities induce more heat production by frictional heating and consequently more basal meltwater lubricating the bed. Our results do not show this feedback. This is probably due to lower sliding velocities in our simulations which do not produce enough friction heat at the base to change the temperature regime significantly, and due to the missing feedback between the flow and geometry evolution of the ice in our model.
Although we run our model with constant forcing and a fixed geometry, all simulations performed with a temperate ice rheology do not run into steady state, independent of whether we allow for basal sliding or not (e.g. CONt, REFt, BWt, QUt and TAt in Fig. 8). The domain-averaged variations of the TIAF, the velocity and the basal temperature in this figure appear small, but the temperate ice rheology simulations lead to significant temporal variations at several areas in the model domain, such as the sticky spot, where the surface velocity varies by a factor of 5 over the 30 ka simulation time (Fig. 9).
As suggested by Reference RobinRobin (1955) and discussed in more detail by Reference Clarke, Nitsan and PatersonClarke and others (1977), the ice is subjected to creep instability, based on nonlinear coupling between the temperature and the flow of ice via the creep relation (Eqn (12)). The parameter values of Reference Paterson and BuddPaterson and Budd (1982) used in the rate factor (Eqn (13)) vary by three orders of magnitude over the temperature range found in the ice. Thus a small increase in deformation rate, by increasing the strain heating and therefore the temperature in the ice, will result in a further increase in deformation rate.
Reference Payne and DongelmansPayne and Dongelmans (1997) applied a 3-D, thermo-mechanical ice-sheet model to an idealized ice-sheet geometry. The work they reported showed that localized creep instabilities can be responsible for the onset of streaming, based only on internal feedbacks within the ice. We infer that creep instability is responsible for the temporal variations of the ice flow in our model. If a creep instability occurs locally, the resulting enhanced flow at the base causes increasing horizontal and vertical advection rates. In particular, the vertical advection of cold ice to the base will decrease the temperature in the basal ice layers and therefore the ice deformation. The whole process gives rise to a cyclic behaviour (cf. non-sliding simulation CONt). The process is further amplified if we also take basal sliding into account. This can be seen in Figure 8, where the simulations with the highest sliding velocities (BWt, QUt and TAt) also show the largest temporal variability for the TIAF, the mean velocity and mean basal temperature. Due to sliding at the base, shearing in the basal ice layers becomes less effective, thus reducing the amount of strain heating in these layers, and, in addition, more cold ice is advected from the surface to the base. The reduction of strain heating at the base caused by sliding cannot be compensated by the additional friction heat and prevents the development of a temperate ice layer at the sticky-spot location, at least for the simulated time period in the BWt simulation (Fig. 9). This is in good agreement with the findings of Reference Blatter and HutterBlatter and Hutter (1991), who concluded that strain heating is the necessary mechanism to produce a temperate ice layer at the base. The parameterization of the temperate ice rheology with a constant water content if the ice reaches the pressure-melting point supports the occur-rence of a creep instability, since ice at the pressure-melting point becomes even softer and enhances the positive feedback of the creep instability. In particular for the relative high water content of 1% (upper limit for which Eqn (14) holds true) the softening effect is largest. For future studies our model would therefore benefit from an explicit solution of the enthalpy and subsequently temperature and the water content, as in Reference Aschwanden, Bueler, Khroulev and BlatterAschwanden and others (2012).
Conclusions
Motivated by the importance of water for the rheology of temperate ice (microscopic water content) and for the basal lubrication (macroscopic water sheet), we examined their long-timescale effects using the 3-D, thermocoupled, full-Stokes ice-flow model TIM-FD3 on a 2.5 km horizontal grid, for an area covering major ice streams in western Dronning Maud Land adjacent to the Brunt and Riiser-Larsen Ice Shelves. We implemented three different flux-routing algorithms for the subglacial meltwater and modified a Weertman-type sliding relation to account for higher sliding velocities under wet basal conditions. Subsequent to full-Stokes spin-up simulations over 200 ka, different sliding simulations, each covering a model time span of 30 ka, were performed. We find the occurrence of the major ice streams in the study area to be mainly controlled by the ice and bedrock geometry. Smaller glaciers draining into the Brunt Ice Shelf only appear as pronounced individual glaciers when lubrication coupled to subglacial water layer thickness is taken into account. We observed in the simulations that sliding with subglacial water routing leads to increased velocities, although the effective sliding area decreases because of the temperature change at the base. Simulations performed with a temperate ice rheology result in cyclic behaviour of the ice flow on millennial timescales at distinct locations in the model domain, which is even stronger if the water layer thickness is considered in the sliding law. Further investigations of the cyclic behaviour are planned for the future and require the solution of the enthalpy evolution to improve model estimates of the water content in the ice.
Acknowledgements
We thank M. Rückamp and F. Ziemen who read drafts of the text and contributed useful ideas and discussions.
Open access
