2.1. Mechanical model: no back pressure

Recent years have seen a reappraisal of the “back-pressure” concept of marine ice-sheet dynamics (Reference PatersonPaterson, 1994; Reference BentleyBentley, 1997, Reference Bentley1998). A number of observations (Reference VaughanVaughan, 1993; Reference Whillans and van der VeenWhillans and Van der Veen, 1993; Reference MacAyeal, Rignot and HulbeMacAyeal and others, 1998) indicate very limited coupling between grounded ice and shelves in a number of cases. These observations are consistent with theoretical considerations discussed by Hindmarsh (Reference Hindmarsh and Peltier1993, Reference Hindmarsh1996), who pointed out that grounded ice masses with sufficiently high basal traction have their mechanics described by the shallow-ice approximation (Reference HutterHutter, 1983), and that significant longitudinal and other non-shallow stress effects were in consequence confined to a narrow boundary layer close to the grounding line. A large number of theoretical studies of the grounding-line area (Reference Herterich, Van der Veen and OerlemansHerterich, 1987; Reference Barcilon and MacAyealBarcilon and MacAyeal, 1993; Reference LestringantLestringant, 1994) support this idea. This means that in modelling marine ice sheets, there are large areas of grounded ice where the modelling of back pressure is not important; this could include all grounded areas excepting the ice aprons at the mouths of ice streams. This means that the sheet model does not need a coupled shelf model to compute the grounded-ice dynamics. A consequence of this is that rather than a limited number of steady solutions for ice-sheet geometry existing, which is true for non-marine ice sheets, an infinite number of steady geometries are possible (Reference Hindmarsh and PeltierHindmarsh, 1993, Reference Hindmarsh1996); in other words, the system is in neutral equilibrium. This has significant implications for modellers, as ensuring dynamical consistency between physical and numerical models of marine ice sheets is somewhat delicate (Reference HindmarshHindmarsh, 1996).

The fact that the back-pressure concept has been reappraised means that there are now two distinct schools of marine ice-sheet numerical modelling: (i) models where shelf and sheet are coupled across the grounding line with a computation of longitudinal stresses in this area (Reference Lingle and ClarkLingle and Clark, 1985; Reference Van der VeenVan der Veen, 1985; Reference HuybrechtsHuybrechts, 1992; Reference PattynPattyn, 1996), and (ii) uncoupled models where grounding-line advance is computed from the kinematics (Reference Hindmarsh and PeltierHindmarsh, 1993, Reference Hindmarsh1996). Models following method (i) can be criticized because they have not yet attempted to demonstrate that the numerical models are not adding spurious elements to the dynamics of the marine ice-sheet system; method (ii) can be criticized because it ignores the subtlety of grounding-line dynamics.

We follow method (ii) and reply to the latter criticism by saying that so far as the large-scale dynamics of marine ice sheets are concerned, complex stress fields in the grounding line are just a subtlety; there is no evidence or theoretical argument that they affect the flow of sheet/inland ice where basal traction is high. While these stress fields may affect streams, they are only of significance for areas of streams with slopes similar to those of ice shelves (Reference Hindmarsh and PeltierHindmarsh, 1993); such areas are of limited areal extent. In any case, the dynamical complexities of ice rises or marine ice sheets without streams are sufficiently complicated that they warrant separate investigation. We do represent sliding in this paper such that a very high proportion of the total velocity is due to slip, and check that the shallow-ice approximation is still valid. This is done by computing the “traction number” (Reference Hindmarsh and PeltierHindmarsh, 1993), which is the ratio of the surface longitudinal stress to the basal shear stress. If this number is around unity, non-shallow-ice approximation effects are unimportant (see Table 3, shown later). The computed traction numbers are never significantly greater than unity.

We do not represent lateral variation in flux (“streaming”), and it is possible that in streams where longitudinal stress gradients play a significant role, the neutral equilibrium concept on which the numerical schemes used in this paper are based does not hold. This issue remains open; thus we assert that the analysis is valid for slow-flowing inland ice, ice streams with mechanics described by the shallow-ice approximation and possibly for all grounded ice. Since we do not model a shelf, the model cannot represent shelf grounding caused by sea-level drop.

2.2. The ice-sheet equation in a spherical system

The continuity equation is

(1)

where H is the ice thickness, a is the mass balance and q is the mass flux. Assuming an axi-symmetric ice sheet resting at the surface (r = R) of the planet, functions are a position of latitude θ, and no dependence on longitude (φ) exists. Figure 1 illustrates the spherical coordinate system (r, θ, φ) in terms of a Cartesian coordinate system (x, y, z) for an arbitrary scalar field M. The divergence of the horizontal flux q in the coordinate system (θ, φ) on a sphere becomes

(2)

where q_θ is obtained from the vertically integrated ice horizontal-velocity profile. Figure 2 illustrates the dependence of thickness field H, surface elevation s and basal elevation b on latitude. Free water depth is denoted f, and the thickness and free water depth at the grounding line are denoted H _m, f _m, respectively. The velocity v has two components (v _r, v_θ). For the moment, we assume ice motion to be exclusively due to internal deformation (no basal sliding), with a non-linearly viscous flow law, meaning that the meridional (outward) flux component q_θ can be written as a function of local ice-sheet thickness H and upper surface radial gradient ∂s/∂θ (Reference MarshallMarshall, 1996),

Fig. 2. Ice sheet resting upon a spherical Earth. Owing to axial symmetry, and the fact that on the Earth’s surface r is equal to the earth radius R, the spatial variables reduce to the co-latitude θ. The span θ _m determines the ice thickness at the grounding line H _m as a function of the local free water depth f _m via the flotation criterion. The ice flow q_θ(θ) is computed using the shallow-ice approximation.

(3)

where C = 2A(T)(ρ_ig) ⁿ /(n + 2), with A(T) a rate factor for the ice, g the gravitational acceleration, ρ _i the ice density and n the exponent for the creep power law (Reference GlenGlen, 1955)

(4)

Here, e_ij is the deformation-rate tensor, and τ_ij is the deviator stress tensor which has second invariant

(5)

By construction, the surface s(θ) and bedrock elevations b(θ) are related by

(6)

(Fig. 2).

Sliding can also be incorporated into the description as follows. By allowing sliding at the bedrock contact, the velocity component at the base of the ice column (V _h(r, θ); see Fig. 2) is no longer zero and is given by an appropriate sliding law. We chose a Weertman-type sliding which relates basal velocity to the basal shear stress σ _b according to

(7)

Integrating this velocity term all along the ice-column thickness H gives , the extra flux due to the basal sliding. By expressing the basal shear stress as

(8)

this extra flux due to basal sliding can be written in a form similar to that resulting from internal deformation in Equation (3)

(9)

Taking ℓ equal to n (3 in the computations), the total flux now becomes

(10)

In our experiments, the value for A _s has been tuned so that the flux due to basal sliding is a multiple α of that due to internal deformation at the grounding line and in the initial steady state. Given Equation (10) and the value for C in the previous paragraph, this finally gives for A _s

(11)

with A(T)the rate factor for the ice. Basal sliding is omitted in this study except where it is specifically mentioned as being included.

2.3. Computation of steady profiles

All of the dynamic computations reported in this paper have steady initial conditions, and examining the steady profile itself is also instructive. The steady profile is obtained from Equation (1) by setting ∂H/∂t = 0. Then, from Equation (2) and the zero-flux condition at the divide, we obtain

(12)

(13)

which can be integrated to give

(14)

We discretize over latitude θ by constructing N + 1 nodes with index k ∊ (0, N). By expressing the flux in terms of a finite-difference formula as a function of the ice thickness H and upper surface elevation s at the two successive nodes (k and k + 1) of the grid, we obtain from Equation (3)

(15)

where the subscript k + 1/2 stands for the flux at the grid centre between the nodes k and k + 1, and Δ _θ represents the difference in co-latitude between two gridpoints. Computation of the steady profile consists of a backward integration from the grounding line up to the ice divide.

The ice thickness at the grounding line H _m obeys the flotation criterion (Reference Van der VeenVan der Veen, 1985; Reference HindmarshHindmarsh, 1996)

(16)

where ρ _w is the sea-water density (1028 kg m⁻³), θ _m is the prescribed span in degrees of latitude and f _m is the free water depth at the grounding line (see Fig. 2). Equation (15) is solved for S_k which, using the flux expression (14), gives

(17)

where θ_{k +1/2} represents the position (k + 1/2)Δ _θ . In Equation (17), (H_k + H _{k+ 1})/2 is first approximated by H_{k +1}, leading to a first estimate of S_k which is then used to yield an estimate of H_k . This latter value is then used for a more accurate estimate of (H_k + H_{k +1})/2 in Equation (17) which is solved again. The iteration is repeated until convergence is obtained. The convergence is extremely fast since a relative variation in the solution less than 10 is reached after two or three iterations. The particular design of this scheme ensures consistency with the numerical schemes used for the time-dependent equations (Reference HindmarshHindmarsh, 1996). When sliding is incorporated, Equation (17) is modified by use of Equations (10) and (11) so that the term [(H_k + H_{k +1})/2 ^{n +2} in the denominator is replaced by [(H_k + H_{k +1})/2] ^{n +1} {[(H_k + H_{k +1})/2] + αH _m}.

Throughout this study, we only consider isothermal ice with a uniform flow law coefficient A(T) varying between 1.5 × 10⁻¹⁸ and 1.5 × 10⁻¹⁷ Pa⁻³ a⁻¹, which, according to Reference PatersonPaterson (1994), corresponds to ice temperatures between about −30° and −10°C, respectively (without accounting for any enhancement factor resulting from fabric). The exponent for the flow law n has been set to 3, and the ice density is supposed to be uniform and equal to 917 kg m⁻³ Different steady-state profiles are described in section 3.

2.4. Computing the response to forcing

2.4.1. Grounding-line migration and moving grid

Following Reference HindmarshHindmarsh (1996), the time-dependent evolution for the ice sheet is solved approximately using finite differences to discretize the continuity equation on a normalized horizontal domain. Thus, we transform the spatial variable by constructing a normalized co-latitude ξ = θ/θ _m such that the ice divide and the grounding line correspond to ξ = θ and ξ = 1, respectively, implying solution over the domain ξ ∊ [0, 1]. It is known that such a transform permits (but does not ensure) consistency of the stability properties of the discretized equations with those of the associated differential equation. The grounding-line (angular) velocity is computed from a total differentiation of the flotation criterion (16) at the grounding line with respect to time; such a form is equivalent to a statement of mass conservation (Reference SalamatinSalamatin, 1989; Reference HindmarshHindmarsh, 1996). Then, by replacing ∂H _m/∂t by its expression in the continuity equation (1), and ∂/∂θ by (1/θ _m) (∂/∂ξ), we can follow the procedure in Reference HindmarshHindmarsh (1996), but now in the different ξ-coordinate system used in this paper, to obtain

(18)

where ∂f _m/∂t is the time derivative of the free-water depth at the grounding line. Since we neglect local variations in sea level arising from variations in the elevation of the geoid, the sea surface is assumed to be uniformly flat throughout time so that (∂f/∂ξ)| _{ξ =1} = (−∂b/∂ξ)| _{ξ =1}. The quantity ∂f _m/∂t is the difference between the eustatic rate of sea-level change and the rate of change of bedrock elevation. Owing to grid motion consequent upon changes in the ice-sheet span , the transformation to the normalized domain introduces an advection term into the continuity equation, leading to the expression

(19)

The discretization using an explicit time-stepping scheme for Equation (19) as well as that for Equation (18) are given in detail in the Appendix.

In the numerical experiments reported below, the forcing consists mainly of sea-level changes which influence the dynamic behaviour via the grounding-line migration (see Equation (18)). Temperature changes also exert a strong control on the rheological properties of ice through the temperature-dependent rate factor A(T) (see, e.g., Reference PatersonPaterson, 1994). However, temperature changes within the ice sheet in response to climatic forcing are neither uniform nor instantaneous, and these complex thermomechanical coupling effects are considered in the companion paper (Reference Hindmarsh and MeurHindmarsh and Le Meur, 2001). Our main concern in this paper is understanding the coupling between ice flow and bedrock deflection.

3.1. Steady profiles

Steady profiles were computed using the iteration (17) for different combinations of the parameters defining ice-sheet and bedrock geometry and the accumulation rate. We first tried to reproduce two major typical sizes with initial ice-sheet semi-spans of 500 and 1000 km, representative of the modern West and East Antarctic ice sheets, respectively. The rate factor A was also adjusted depending on the size of the ice sheet. The larger ice sheet is assumed to consist of relatively cold ice (≃ −30°C with a corresponding rate factor A = 1.5 × 10⁻¹⁸ Pa⁻³ a⁻¹), while the smaller ice sheet had a rate factor of 1.5 × 10⁻¹⁷ Pa⁻³ a⁻¹ roughly corresponding to −10°C. Basal sliding was also included in some of the calculations for the warmer ice sheet. This requires a slightly modified computation as described in section 2.3. The effects of the depth of the base (in particular the depth at the grounding line) were also investigated for the smaller ice sheet with values of both 500 and 1000 m below sea level used.

In addition, the consequences of different accumulation distributions were considered. We first considered uniform accumulation rates all over the ice-sheet surface, choosing 0. 03 and 0.05 m a⁻¹ as accumulation rates for the larger and smaller configurations, respectively. We also introduced spatial variation into the patterns by constructing them to increase exponentially from the ice divide to the coast. Galling a ₀ the accumulation at the divide, and assuming the accumulation rate to be a ₀ k _a at the margin (grounding line, θ = θ _m), the general expression for accumulation rate as a function of co-latitude becomes

(20)

In that case, the steady continuity equation (12) becomes

(21)

Integrating by parts twice yields

(22)

and application of the zero-flux condition at θ = 0 yields

(23)

We then obtain

(24)

This new expression is used in Equation (17) in exactly the same way as previously described.

The different parameter combinations used here are summarized in Table 1. All the different accumulation parameterizations are labelled with the corresponding (a ₀, k _a) combination (note that k _a = 1 gives a uniform accumulation rate). Where sliding is specified, the ratio α = 1. Corresponding computed profiles are shown in Figure 3.

Table 1. Parameter sets used in Figure 3

Fig. 3. Initial steady-state ice-sheet profiles corresponding to parameter sets defined in Table 1 and computed according to Equation (17).

All the profiles exhibit the characteristic parabolic shape for steady-state ice sheets. Since the marginal accumulation is greater, the configuration where the accumulation rate varies in space leads to an increased total ice input and therefore a thicker ice sheet (compare runs 1 and 2 as well as runs 4 and 6), despite the fact that the same accumulation rate obtains at θ = 0. Allowing for basal sliding for ξ ∊ [2/3, 1] increases the overall ice velocity which, under the same accumulation pattern, gives a thinner steady ice sheet (compare runs 2 and 3). Run 5 uses a rate factor A = 1.5 × 10⁻¹⁸ Pa⁻³ a⁻¹ one order of magnitude smaller than run 4. It indicates that the resulting colder ice, by being more viscous, deforms less readily and leads to a thicker ice sheet. Finally, run 7 shows the effect of a shallower depth for the bedrock. In comparison with run 6, we see that the effect of increasing the bedrock elevation by 500 m increases the ice surface elevation by almost the same amount close to the divide. Since the two systems are in steady state, the necessity of maintaining the same flux through the decreased thickness required by the flotation condition implies that the margin slope is much greater for the shallower sea (as already shown in Reference HindmarshHindmarsh, 1996). This can in principle have an influence on the evolving ice sheet, since rate of expansion depends on the ice surface gradient (see the term (∂H/∂ξ)| _{ξ= 1} in Equation (18).

3.2. Time-dependent response

3.2.2. Sea-levelforcing

We now consider transient effects induced by sea-level changes. Sea-level fall and rise at a constant rate of 13 mm a⁻¹ are considered over a period of 10 kyr. This replicates the 130 m sea-level change between glacial and interglacial epochs; the rate of rise is comparable with the average rate during deglaciation. The computed profiles as well as the initial steady profile for the two ice-sheet sizes and bed configurations are shown in Figure 4. For the larger ice sheet, the other parameters are as for run 2 (Table 1), whereas the configuration for the smaller ice sheet is as for run 6, but basal sliding is not permitted. When sea level increases, the outermost part of the ice sheet is floated off, thereby leading to a grounding-line retreat. Conversely, a sea-level drop permits a grounding-line advance.

Fig. 4. Grounding-line advance or retreat after 10 kyr in response to sea-level fall or rise for both the small (500 km; b, d) and the large (1000 km; a, c) ice sheet with the two different bedrock depths (−1000 m (a, b) and −500 m (c, d)). The initial steady profile is plotted with a dashed line, as are the final sea-level stands at the end of the simulation (±130 m). No bedrock deflection is permitted in these calculations.

At this stage, it is worth pointing out the difference induced by uniform and spatially varying accumulation parameterizations on the ice-sheet evolutions driven by sea-level change. Owing to the local character of the flux equations (3) and the ice-thickness change at a given point in space (Equation (19) without the advection correction term , changes in the grounding-line position cannot modify the inland ice-sheet profile, provided the accumulation is spatially uniform. On the other hand, the exponential dependence of the accumulation term on normalized position leads to a time-varying accumulation term at a given location as the span changes, which in turn will perturb the ice-sheet profile (Fig. 4). In other words, the accumulation rate varies as a function of normalized position, which means the accumulation rate at a given position in physical space changes as the grounding line moves. Ultimately, this causes the ice-sheet elevation to change in the grounded area. Owing to the relatively small change of the ice-sheet span, this change in thickness is not very pronounced in Figure 4.

Some cases where this effect is more obvious are shown in Figure 5, which shows the effect of spatially uniform and non-uniform accumulation rates upon ice-sheet retreat in response to sea-level rise (the configuration here corresponds to Figures 4b and 6b). Accumulation rates are the same (50 mm a) at the divide, but are 10 times larger at the margin for the non-uniform case. This means that the margin flux for this case is larger, meaning that the slope is larger and the whole ice sheet is thicker. The larger slope also means a lesser retreat. In the uniform-accumulation-rate case (case 1, Fig. 5), ice-sheet elevation does not change in grounded areas owing to the local form of the flux expression arising from the local form of the driving stress given by the shallow-ice approximation. In the non-uniform case (case 2, Fig. 5), as a result of the construction of the accumulation distribution, changes in the grounding-line position cause changes in the accumulation distribution within the ice sheet, which means that the ice-sheet elevation in the grounded area changes as well. In order to produce a substantial retreat, a very high sea-level change was specified (total rise of 530 m). The fact that the uniform-case profile does not change shows that the numerical method is working successfully. These experiments were repeated with twice as many gridpoints (200), and very similar results were obtained (see the corresponding grey curves, which superpose almost exactly on the 100 gridpoint black curves in the zoom box).

Fig. 5. Effect of spatially uniform (case 1) and non-uniform accumulation rates (case 2) on ice-sheet retreat in response to sea-level rise (the parameter set here corresponds to Figures 4b and 6b). These experiments were repeated with twice as many gridpoints (200) (corresponding grey curves, which superpose almost exactly on the 100-gridpoint black curves in the zoom box).

Reference HindmarshHindmarsh (1996) observed a significant change of ice thickness at the ice divide by exclusively forcing his model by sea level, even though a uniform accumulation rate was specified, which is an erroneous result. This problem was solved by using higher-order finite-difference methods (see Appendix) which can be shown to produce dynamically consistent numerical schemes using a suitable modification of the techniques described by Reference HindmarshHindmarsh (1996, appendices 1 and 2).

Figure 4 also shows the strong effect of the free water depth at the grounding line where a shallower bedrock induces a greater ice-surface gradient. This effect tends to reduce the grounding-line migration arising from sea-level forcing. To see why, we note that since the flux is locally determined and its gradient is approximately in equilibrium with the local accumulation rate, if forcing is due to sea-level change only, the grounding-line formula (Equation (18)) reduces to

(25)

when the bedrock is flat. (In fact, for display purposes, the accumulation rate in Figure 4 is exponentially varying in space in order to prevent the different profiles from being coincident, but the flux gradient continues to more or less balance the accumulation rate, and the profiles are very close to steady-state profiles.) It is therefore clear from Equation (25) that for a given rate of sea-level change ∂f _m/∂t (recall there is no bedrock change being considered at the moment), the grounding-line motion will be inversely proportional to the marginal ice-thickness gradient (1/θ _m) (∂H/∂ξ)| _{ξ= 1}. The computed grounding-line migration rates are extremely small; they can also be shown to be consistent with the initial margin-thickness gradient in Table 2 and Equation (25). The small retreats are a consequence of the fact that the thickness gradients are quite high owing to the ice being cold and therefore stiff. However, the range of slopes we compute for both steady-state and time-dependent results (see Tables 2 and 3) are consistent with observations of slope, which show a range of 0.1% to several per cent (Reference McIntyreMcIntyre, 1983; Reference Jackson and KambJackson and Kamb, 1997). For some context, notice that a thickness gradient of only 0.1% will produce a retreat of 112 km for a 100 m sea-level rise; such a retreat is still small compared with that observed for the Siple Coast (Reference Licht, Jennings, Andrews and WilliamsLicht and others, 1996; Reference Conway, Hall, Denton, Gades and WaddingtonConway and others, 1999).

Table 2. Steady ice-sheet thickness gradient at the margin and corresponding grounding-line retreat computed from Equation (19) integrated over the 10 kyr

Table 3. Effects of including basal sliding for case B in Figure 4 and in Table 2 in terms of total retreat after 10 kyr, ice surface slope at the grounding line and traction number

We now consider cases with the same free water depth, but for the two different ice-sheet spans (compare Fig. 4a with Fig. 4b, or Fig. 4c with Fig. 4d). When the span is larger, the grounding-line flux is also larger, leading to a greater slope at the grounding line. This effect is compounded by the fact that we have specified the experiment such that the larger ice sheet also has a smaller rate factor. This slope difference has predictable effects on the total displacement of the grounding line, as well as the respective grounding-line velocity and ice-sheet span for the same cases (shown in Fig. 6).

Fig. 6. Same experiments as in Figure 4, showing grounding-line velocity (solid line) and ice-sheet span (dashed line) during the 10 kyr of simulation. Retreats (negative velocities) and advance (positive velocities) computations are shown for each case.

Table 3 shows the effects of including basal sliding for the configuration corresponding to Figure 4b. One can see that basal sliding increases the total retreat after 10 kyr, and it also reduces the slope such that the geometrical considerations described above still apply. Although sliding makes the ice sheet more sensitive to sea-level rise, the total retreat after a sea-level rise of 130 m over 10 kyr still remains limited.

Also noticeable from the total displacement after 10 kyr in Figure 4 is the asymmetric grounding-line motion in response to sea-level rise and fall, respectively. The ice sheet shrinks slightly more easily than it expands. This is readily interpretable once one realizes that the ice-sheet expansion is limited by the necessity of providing the extra ice volume through accumulation, whereas nothing similar can prevent the ice sheet from retreating as fast as sea level dictates. The features are demonstrated in Figure 6, where the grounding-line velocity and the ice-sheet span are plotted against time. When considering retreat, the grounding-line speed steadily increases, whereas with grounding-line advance, the velocity reaches a maximum generally after about 1 kyr and subsequently decreases. This asymmetric behaviour appears to depend upon the relative rates of accumulation and sea-level change. For example, in case A, the discrepancy between the retreat and advance magnitudes which is 13.5% with a sea-level rate of change of 13 mm a⁻¹, is only 6.8% when the rate of sea-level change is halved.

A final point is that given that the speed of retreat increases in all cases, it may be of interest to address the longer-term evolution of the ice sheet. However, if one wants to reproduce a realistic deglaciation forced by sea level, either a non-zero rate of change of sea level has to be specified for the first 10 kyr only (at the same rate of 13 mm a⁻¹), or it needs be specified to last for longer but at a lower rate. In the first case, after 10 kyr of sea-level rise (13 mm a⁻¹), the speed of retreat rapidly falls to zero within < 200 years, leading to the same total retreat after 50 kyr as after 10 kyr (in fact there is a slight decrease in the total retreat of 40 m which is due to a slight readvance during the last 40 kyr of constant sea level). In the second case, a halved rate of sea-level change over twice as long a period (20 kyr) leads to the same final retreat with rates which are exactly halved. Similar results were also obtained in case of ice-sheet advance in response to sea-level fall.

6.1. Bedrock model

The results above have made it clear that any process modifying depth and slope at the grounding line through isostatic Earth accommodation under surface-varying loads will play a role in the grounding-line dynamics.

The bedrock model used here (Reference Le MeurLe Meur 1996a, b) is based upon previous approaches (Reference LongmanLongman, 1962; Reference FarrellFarrell, 1972; Reference Wu and PeltierWu and Peltier, 1982) which led to a new generation of SGVE spherical Earth models. Self-gravitation is a fundamental aspect, and the perturbation in the gravity field induced by deformation is not only explicitly computed but also fully included as an extra driving force in the stress-balance equations. The spherical approach is effected by a spectral decomposition onto a basis comprising the spherical harmonic system of orthonormal functions. The spherical approach allows proper account to be taken of the radial distribution of the Earth’s mechanical properties, especially when computing the far-field response. The model incorporates the entire mechanical structure of the planet from the outer elastic lithosphere down to the core with a three-layer mantle lying in between, where most of the isostatic accommodation takes place. The Earth’s structure is summarized in Table 4. A significant feature is the use of a Maxwell body model to represent the viscoelastic response of the mantle; this introduces the characteristic time-scales which are mainly responsible for the Earth dynamics. The viscoelastic approach is formulated using the conceptual framework of the “correspondence” principle (Reference BiotBiot, 1954; Reference PeltierPeltier, 1974). The Green’s functions cause the solution to have an instantaneous elastic response and several relaxing viscous modes. The response to any loading history can be obtained by a double space- and time convolution of these Green’s functions with the space and time distribution of both ice and water overburdens.

Table 4. Earth parameters used in the unit bedrock model

For the sake of more accurate results as well as a faster computation, the model has been coded with the symbolic manipulator Mathematica (Reference Spada, Yuen, Sabadini, Morin and GasperiniSpada and others, 1990, Reference Spada, Sabadini, Yuen and Ricard1992). By replacing approximate numerical calculations with exact symbolic operations, many of the propagating round-off errors can be avoided in the numerous matrix products involved. Given an Earth structure, the corresponding set of Love numbers (amplitudes in the spectral domain from which Green functions are summed later) represents the actual output for the unit model. This only needs to be computed once.

6.2. Asynchronous coupling between the two models

Following Reference Le Meur and HuybrechtsLe Meur and Huybrechts (1998) it was found that computing the bedrock evolution every 100 years was a sufficiently fine temporal resolution to reproduce correctly the ice-sheet/Earth dynamics (the time-step for the ice model being about 0.2 year). This time-scale emerges as a consequence of the characteristic time-scales for the viscous mantle response and is a trade-off between computational effort and accuracy.

Following summation in a Legendre series of the elastic and viscous Love numbers computed by the unit bedrock model (see, e.g., Reference Le Meur and HindmarshLe Meur and Hindmarsh, 2000), the impulse Green function response reads

(26)

where G ^E(θ) represents the elastic instantaneous Green’s function, and each of the jth viscous-mode Green’s functions. From there the bedrock response at time t for any of the kth gridpoint along the flowline (and also at sea) can be expressed as

(27)

where the summation represents the spatial convolution. This consists of integrating the contribution to deformation of all the neighbouring points lying within a circle of radius 1000 km. The circle (symbolized as D_k ) is discretized along a local latitude (index j_α ) and a local longitude (index j_β ), and is the angular distance between the gridpoint k and the point whose contribution is considered. The element area is Δ_∑ = R sin is θΔ _θ Δ _φ where Δ _θ ,Δ _φ are the grid dimensions on the surface of the sphere. It should be noted that lithospheric rigidity causes the steady bedrock response to deviate from local compensation, being now partly driven by the neighbouring loads.

The time integration expresses the fact that the current response does not exclusively depend upon the current state of loading, but also includes the contribution from past loads (see H(i_α, j _β , t′) in Equation (27)) because of the delayed viscous response. A memory period (Δ_mem) therefore has to be chosen according to the particular viscosity structure of the model, so that the viscous response is completely accounted for. In our case, with the viscosity values adopted (see Table 3), it appears enough to convolve 15 kyr back in time in order to retain > 98% of the viscous response.

Every 100 years, the loading function H(i _α , j _β , t) is updated from the current ice distribution and sea level, and is used to compute the new bedrock profile according to Equation (27). Bedrock response at sea is computed up to an appreciable distance (1.5 × the ice-sheet span) from the ice-sheet edge, since change of the depth of the seabed induces a loading variation by modifying the water column. In this respect, local sea-level changes induced by geoidal undulations should theoretically also be taken into account. The need for such a level of sophistication is, however, questionable since it would require a global approach for solving the sea-level equation (e.g. Peltier, 1994), with limited gains compared with the accuracy of both the ice-sheet and the bedrock models. Finally, the newly computed deflection rate is used during the succeeding bedrock interval (100 years).

7.1. Specification of initial conditions for the bedrock deflection

The uncoupled SGVE model has the dynamics of a linear system; in other words, it relaxes exponentially towards that isostatic equilibrium corresponding to the current loading. The fact that linear systems admit superposition of solutions gives us some flexibility in formulating initial conditions. One possibility (method a) is to choose an ice-sheet span, compute the profile on a flat bed, determine from this the loading, and compute the steady bedrock deflection induced by this load. With this deflection and with the same span, one computes the new ice-sheet profile, and proceeds iteratively until convergence is obtained. In practice, around 20 iterations are needed for the same level of accuracy as in the computation of the steady profile in section 3. This approach has the merit of conceptual simplicity, as we simply compute the steady deflection of an initially flat bed, and the approach shows clearly what the equilibrated Earth response to an ice sheet means in terms of bedrock displacement. This aspect is considered in section 7.3 where the bedrock response computed according to method a reveals some interesting geometrical features.

In fact, as a consequence of superposability, a flat bed is not a unique choice for a bed in isostatic equilibrium with a spatially uniform load. Superposability means that we can choose any bed profile and define it to be in isostatic equilibrium with the initial ice-sheet profile. Lithospheric deflections due to subsequent loading changes are readily computable using superposition principles. We can justify the approach by pointing out that in general the bedrock profile corresponding to ice-free conditions will not be a flat bed. We term this method b. We shall carry out investigations using these two methods.

7.2. Time-dependent response

The same sensitivity tests on the depth and the slope at the grounding line as for a rigid bed (see section 5) were carried out, but with the bedrock deflection included using method b. The results (darker contours) are displayed along with those computed for the rigid Earth (pale grey) in Figure 8. For the parameter set considered, the inclusion of iso-statically induced bedrock deflection results in an increase in the total computed grounding-line retreat. This increase becomes more pronounced for larger depths and, to a lesser extent in this range of parameters, for greater slopes.

Figure 9 shows some ice-sheet and bedrock profiles arising from these coupled evolutions. The 13 mm a sea-level rate of increase over the simulation leads to a total rise of 130 m which significantly depresses the bedrock beyond the grounding line. Underneath the ice sheet, the bedrock deflection hardly changes. As shown in the lower zoom box and in the upper inset, the bedrock response noticeably increases the final water depth at the grounding line (from 1213 m to 1245 m, i.e. by about 15%), whereas it has almost no effect on the slope (change is less than one-fiftieth). Note that since there are no physical mechanisms represented in the ice-sheet model which allow substantial change of the inland ice sheet, no isostatic recovery is expected and the consequence of deglaciation of distal areas (e.g. in the Northern Hemisphere) is increased bedrock deflection at the grounding line via sea-water loading.

Fig. 9. Simulations of coupled ice sheet and bedrock carried out with the initial state corresponding to case 1 in Figure 7 using method b. This figure shows final geometry after 10 kyr of simulation when the bedrock response is either accounted for (black) or disregarded (grey).

Owing to the fact that changes in the total loading around the grounding line are relatively small (the ice which has been floated off is replaced by a water column of similar weight), the bedrock response there is mainly the result of the regional bending of the lithosphere in response to the sea-level change and to changes in the loading beyond the grounding line (Fig. 9). A consequence of this is that the local slope is hardly modified, whereas the free water depth at the grounding line increases rather more (see the enlargement and the time-dependent evolution of these two fields for both cases in the upper part of Figure 9). The effect of the change in bedrock deflection is to increase the rate of change of free water depth, which has a similar effect to the sea-level forcing. This effect is directly detectable by inspection of the grounding-line migration-rate formula (Equation (18)) where the rate of change of bedrock elevation ∂b/∂t associated with the rate of sea-level change contributes to the rate of change of the free water depth ∂f _m/∂t.

7.3. The steady-state bedrock response

Finally we investigate the coupled response of ice sheet and the Earth using method a, iteratively computing a steady deflection corresponding to the ice-sheet loading; we now assume that the ice-free bedrock profile is a flat bed (Fig. 10a) or with pre-existing topography as in Figure 9 (see Fig. 10b). Note that ice-free conditions imply that the bedrock is now submerged, meaning that the additional load we actually consider is that due to the ice freeboard. Both panels in Figure 10 show the resulting steady configuration when bedrock accommodation is either neglected (grey) or accounted for (black) according to method a as previously described. In particular, they show the effects of the rigidity of the lithosphere on the overall geometry of the bedrock response under an ice sheet.

Fig. 10. Computed steady bedrock deflections according to method a. This figure shows the regional nature of the deflection, which occurs as a consequence of lithospheric rigidity. Rigid Earth (no isostatic response) and corresponding ice-sheet profiles (with or without isostasy) are also shown. (a) Ice-free equilibrated flat bed; (b) Ice-free bedrock with the same relief as in Figure 9.

The response is regional in the sense that it not only depends on the local load but also includes the contribution from all neighbouring points lying within a certain distance. In other words, since the sea level is kept constant, the bedrock deflection due to the presence of inland ice extends beyond the ice-sheet margin. A corollary is that the rigidity of the lithosphere reduces the maximum deflection compared with local compensation according to Archimedes’ principle. For example, with an ice/asthenosphere density ratio of 917/3350, local compensation would give a maximum depression below the ice divide of about 515 m instead of the 368 m we obtained (see Fig. 10a).

The more rigid the plate (generally the thicker), the smoother the overall profile (less pronounced amplitude and more widely spread features). In particular, this means that the slope at the grounding line is strongly controlled by the lithospheric thickness, a parameter which is unfortunately poorly constrained and which probably undergoes strong lateral variations.