Introduction
The volume of marine ice sheets depends on the ice thickness H at the grounding line. Taking H as the boundary value, inland ice thickness evolves such that the surface slope can drive an ice flow sufficient to drain the total accumulation. This view holds for a simple steady state, an isotropic rheology and a frozen bed, in particular. In general, basal melting due to a geothermal heat flux and a transient occurrence of water at the ice-sheet base makes the system more difficult to understand. Ice shelves contribute a buttressing effect as a result of lateral shear forces at their margins. By exerting a back-stress, the ice shelf’s geometry and shape codetermine the thickness and position of the grounding line (e.g. Reference Dupont and AlleyDupont and Alley, 2005).
As early as 1978 Mercer drew attention to the ‘threat of disaster’, a possible collapse of the West Antarctic ice sheet following a climate-induced disintegration of the hitherto buttressing ice shelves (Reference MercerMercer, 1978). A first hint was measured after the disintegration of the Larsen B ice shelf. Reference Scambos, Bohlander, Shuman and SkvarcaScambos and others (2004) reported further retreat of the grounding line after the break-up in 2002. Similarly, Pine Island Glacier (PIG) showed retreat of the grounding line together with upstream thinning and increasing balance velocities (Reference Payne, Vieli, Shepherd, Wingham and RignotPayne and others, 2004). PIG did not yet disintegrate, but was exposed to high ocean-induced melt rates causing the thinning of its floating part. Subsequent heave from a sill led to a reduction in basal stresses and decreased buttressing.
Hydrographic conditions in the Amundsen Sea reveal relatively warm water masses (Reference Hellmer, Jacobs, Jenkins, Jacobs and WeissHellmer and others, 1998) which are held to favour further extensive melting in the PIG area (Reference Thoma, Grosfeld and LangeThoma and others, 2006). Anticipated melt rates of 33ma–1 (Reference Jacobs, Jenkins, Giulivi and DutrieuxJacobs and others, 2011) would lead to further retreat of the grounding line and, due to an inland-sloping bedrock, could initiate a substantial loss of grounded ice (Reference Payne, Vieli, Shepherd, Wingham and RignotPayne and others, 2004).
To quantify ice-volume changes and grounding-line migrations, we apply the dynamic three-dimensional (3-D) ice-sheet/ice-shelf model RIMBAY (Revised Ice sheet Model Based on frAnk pattYn; Reference Thoma, Grosfeld, Smith and MayerThoma and others, 2010a) to a regular model domain, in the first instance. The model is forced by surface accumulation and by basal melting in the floating part. We perform time-dependent simulations in order to investigate volume changes relevant to sea-level changes. Starting at steady state, the system’s reaction after 1000 years of constant melt rates is monitored. Thus, we prescribe evenly distributed melting as well as spatially varying melt rates derived from our ocean circulation model ROMBAX (Revisited Ocean Model based on Bryan And CoX; Reference Thoma, Grosfeld and LangeThoma and others, 2006). The latter is not fully coupled over the whole model period, but is invoked for 3 years to simulate basal melting according to the initial ice-shelf draft. Here our focus is on the effect of distributed melting on the 3-D figure of our ice compartment. The mode of operation using RIMBAY and ROMBAX as a coupled model is described by Reference Thoma, Grosfeld, Smith and MayerThoma and others (2010a).
The Ice Model
Our 3-D ice-sheet/ice-shelf model RIMBAY is based on the dynamic–thermodynamic model constructed by Reference PattynPattyn (2003). A later version of RIMBAY is described by Reference Thoma, Grosfeld, Smith and MayerThoma and others (2010a); several changes, especially concerning the numerical treatment, have been made (Thoma, unpublished information). In its most complex version, RIMBAY can solve the momentum balance by means of the full-Stokes equations (Reference PattynPattyn, 2008). By solving the temperature field within the ice, fully coupled thermomechanical simulations are possible.
Instead of using the model in its highest complexity, we employ the commonly used approximations, namely the shallow-ice approximation (SIA) and the shallow-shelf approximation (SSA). Given the small ratio of ice thickness H to horizontal extent L in the presented experiment, this introduces only negligible errors. A comprehensive description of the SIA and SSA is given by Reference Greve and BlatterGreve and Blatter (2009).
The basal boundary condition for ice sheets and ice shelves determines the flow condition and hence the flow field of the ice body. In the case of a frozen bed, the entire flow is driven by internal shearing. At this time, no attempt is made to incorporate sliding into the simulations. Both SIA and SSA contain the depth-averaged temperature-dependent coefficient A and the exponent n from Glen’s flow law (Reference PatersonPaterson, 1994):
where 
denote the components of the strain tensor, τ′
ij
are the stress-deviator tensor components and τ is the effective stress. Because our perturbation studies span timescales of the order of 103 years, we can omit the solution for the temperature field. For ice sheets the latter would attain a stationary state after 104–105 years. This assumption allows use of a constant A, still in the perspective of the experiments described in this paper. We chose a constant A = 3.0×10–17 kPa–3 s–1, corresponding to an ice temperature of 263 K, and the exponent n=3.
The SIA and SSA yield the velocity field from a given ice-thickness distribution. While the SIA offers a local solution of the flow components u and v, the SSA can be solved as a boundary-value problem. At the ice-shelf front, one can specify the velocity gradient according to (Reference MayerMayer, 1996)
where ε is the effective strain rate and H the ice thickness. We use an ice density ρ ice of 915 kgm–3; the sea-water density ρ w is taken as 1028 kgm–3. At the ice shelf’s inland boundary, the grounding line, the velocities from the last grounded ice-sheet node are taken as the boundary value. Hence, no attempt is made to resolve the exact grounding-line position within a gridcell. The error that could thus be introduced affects the value of the horizontal strain rate in the first SSA gridcell downstream of the grounding line. Subsequent thinning due to an overestimated velocity would cause retreat of the grounding line, to the next gridpoint at most. In this case, the mass flux from the ice sheet could be correctly transferred into the ice shelf.
In time-dependent simulations, the continuity equation yields temporal ice-thickness changes according to
where a denotes surface accumulation, m is the basal melt rate and U the depth-averaged flow velocity. In the case of basal freezing, m becomes negative.
All experiments are performed on a rectangular model domain consisting of an ice sheet, ice shelf and an open-ocean compartment. The bedrock topography resembles a depression opening to one side of the domain. Water depths range from 200 to ~600m at the open-ocean side. The model domain has 130×100 evenly spaced nodes with a horizontal grid size of 5 km ×5 km. Initially, it is covered by a 350m thick slab of ice except for the outermost 30 gridlines covering the area with higher water depths. Due to a flotation condition, we find ice-shelf nodes where H<H bed ρw/ρ ice (Fig. 1). During time-dependent simulations, Eqn (3) is solved at different boundary conditions for accumulation a and melting m. We apply a constant accumulation rate of 0.25mw.e. a–1, while basal melting varies in the experiments (Table 1). The time-step for all experiments is 0.25 years.
Table 1. Comparison of model results. Given melt rates have been active from T = 16.000 to T = 17.000. Numbers for OceanA and OceanB represent average values
Fig. 1. Bottom topography (a) and initial ice surface elevation (b) at time T = 0. Values are given in metres. The ice-sheet/ice-shelf system gains mass from flux across the inland boundary (3000m2 a–1) and from surface precipitation (0.25mw.e. a–1).
The model does not handle the possibility of an advancing ice-shelf front. Instead, all ice passing its initial position (x = 500 km over most of the domain except along the edges) would calve as iceberg. Contrary to this constraint, the ice shelf may retreat if the continuity equation (Eqn (3)) yields a negative ice thickness. All perturbation experiments start from a steady state, which was taken from a 16.000 years spin-up run. The corresponding ice volumes evolve as shown in Figure 2. Figure 3 shows the ice thickness distribution and ice flow at that time.
Fig. 2. Time-dependent ice volume. Spin-up of the ice model and perturbation test. Ice volume reaches a quasi-steady state at T = 16.000 years. An imposed melt rate of 0.5 ma–1 from T = 16.000 to T = 17.000 years causes a volume loss of ~3.3%. The system needs ~5.000 years to recover after melting is switched off.
Fig. 3. INITIAL. Surface elevation at t = 16.000 years in metres (a). The grounding-line position coincides with the onset of small surface gradients. Velocities reach 1150ma–1 at the ice-shelf front. Arrows indicate flow direction (b).
The Ocean Model
The ocean model ROMBAX is a primitive-equation Ocean General Circulation Model (OGCM), which has been applied to sub-ice-shelf cavities and adjacent oceans since the early 1990s in a precursor version. A detailed description is given by Reference Grosfeld, Gerdes and DetermannGrosfeld and others (1997) and Reference Thoma, Grosfeld and LangeThoma and others (2006). Successful simulations used ROMBAX to model melting and freezing under prominent Antarctic ice shelves, as well as the circulation in sub-ice lakes (Reference Thoma, Grosfeld and MayerThoma and others, 2007, Reference Thoma, Grosfeld, Mayer and Pattyn2010b). The latest version simulates a thermohaline circulation driven by buoyancy forces, but forcing by wind or sea-ice formation is also taken into account.
A major driving mechanism for ocean currents within an ice-shelf cavity is the heat exchange across the ice/ocean boundary. Due to the pressure dependence of the freezing point of sea water according to Reference Foldvik and KvingeFoldvik and Kvinge’s (1974) formulation,
a 400 m deep ice-shelf base would melt at ocean temperatures of about –2.2˚C. T 0 = –1.89˚C is the freezing point of sea water (at a salinity of 34.7) at the ocean surface, and p (z) is the pressure at depth z. Relatively warmer currents with respect to the in situ freezing point provide heat for melting, become fresher, and tend to rise, due to mixing with fresh water, to shallower areas. Here, being supercooled with respect to a shallower boundary, the water has the potential to form ice platelets, which will accrete at the ice-shelf base. In ROMBAX, the process of ice/ocean interaction is implemented according to the three-equation formulation of Reference Holland and JenkinsHolland and Jenkins (1999), applying heat- and salt-conserving formulations with respect to the empirical function for the freezing point. Here, frazil-ice formation is not explicitly modelled, but expressed as a bulk formula.
The growth of ice platelets in conjunction with the release of latent heat from this phase transition makes the upper ocean layer warmer and more saline. These water masses could then descend and recirculate or could leave the cavity as Ice-Shelf Water (ISW). At the continental-shelf break, ISW contributes to the formation of Bottom Water when mixing with branches of Circumpolar Deep Water.
All water nodes, either of the open ocean or of the sub-ice-shelf cavity, belong to the ocean model. Unlike RIMBAY, which works on a Cartesian grid, ROMBAX uses geographic coordinates. The horizontal resolution of 0.05˚ in latitude ×0.2˚ in longitude corresponds to gridcells of approximately 5 km ×5.7 km at 75˚ S. In the vertical direction the water column is divided into 14 sigma levels in the open ocean, while 10 levels extend into the ice-shelf cavity.
Using two grid types means that every coupling step of both models requires a coordinate transformation. We use a transformation such that the open-ocean boundary faces the eastern end of the model domain which extends to 76˚ S.
At this boundary, a prescribed mass transport stream function Ψ= 4 Sv (1 Sv = 106m3 s–1) forces a cyclonic circulation within the ice-free ocean area. Respectively, vertical profiles of the potential temperature Θ, and salinity S, kept constant along the open-ocean boundary, serve as boundary values for tracers. In order to keep the experiments simple at this stage, we do not allow temporal changes at the ocean’s eastern boundary. Although ROMBAX would allow for incorporating wind action as well as sea-ice formation and flow, we chose to force the ocean just by prescribing boundary values of Ψ, Θ and S. However, within the model domain, ice-shelf/ocean interaction will feed back on the circulation until a steady state evolves.
Experiments
Our main objective is to investigate the response of the ice-sheet/ice-shelf system to changing melt rates at the ice-shelf base. We focus especially on the movement of the grounding line and the change in grounded ice volume. Thus, we specify two experiments using uniformly distributed melting, and two tests with melt rates from the ocean model. The latter are chosen to investigate the effect of spatially varying melting on the 3-D-ice distribution. One experiment is made to test whether the system readjusts to a pre-perturbation state. If this is the case, coupling across the grounding line seems well defined in our model. Table 1 gives an overview of the experiments concerning ice volume, extent and thickness.
Initial
All experiments start from the result of the control run INITIAL at a quasi-steady state which is attained after 16.000 model years. At this time the model has spun up to an ice volume of 173.800 km3 (Fig. 2). The corresponding thickness distribution and the vertically integrated flow can be inferred from Figure 3. All contributions of ice across the boundary from a virtual catchment area and from precipitation are balanced by calving into the ocean. Ice thickness varies from >1000m at the inland boundary to <200m at the ice-shelf front. Here the magnitude of the ice flow reaches 1150ma–1 as a maximum.
Melt05
To test whether the model will readjust after a perturbation terminates, all floating points are exposed to a uniform melt rate of 0.5ma–1. After 1000 years, the ice volume has decreased by ~3.3%. Accordingly, the grounding-line position has retreated by three grid nodes (15 km), and the ice-sheet area decreased by almost 5.000 km2. When basal melting is switched off again, the volume needs longer than 1000 years to readjust. This is the case after ~5000 years (Fig. 2). After this time the grounding line also reaches its former position.
Melt20
Unlike Melt05 where moderate melt rates are of the order of the surface accumulation, a much higher rate of 2.0 ma–1 is chosen to impose a stronger perturbation. After 1000 years this case shows a strong effect on ice-volume loss. The total ice volume decreases by 25%. Maximum ice thickness in the ice sheet reduces by >5%. The ice-shelf front has retreated by 150 km in the central area while the grounding line has moved inland up to 90 km. Due to balance fluxes from the ice sheet, the ice shelf remains as a 125km wide rim along the grounding line (Fig. 4). Compared to INITIAL, the ice sheet has also lost 20% of its area, which accounts for a loss of ~37.000km2 at this stage (Table 1).
Fig. 4. Ice-surface height for Melt20. Due to a melt rate of 2.0 ma–1, the grounding line has retreated by 90 km within 1000 years. About 900 ice-shelf nodes (>22.000km2) have disappeared beyond the 0-contour in the white area.
A comparison of ice-thickness changes in Melt05 and Melt20 is provided in Figure 5. As expected, the impact of Melt20 on the thickness profile is apparent. Although high melt rates exceeding 2.0 ma–1 have been reported in some areas (Reference HellmerHellmer, 2004), we try to quantify these values independently using ROMBAX. Therefore, we designed two ocean scenarios, yielding melt rates in time and space, driven only by an initial θ,S-distribution in our ocean model.
Fig 5. Resulting ice-thickness profiles at T = 17.000 compared with the initial state at T = 16.000 (blue). Melt05 melt rate is 0.5 ma–1 for the entire ice shelf (Melt20 = 2.0ma–1). The grounding line for Melt20 is situated at x = 45 (=225 km, red arrow) compared with x = 63 (=315 km, blue arrow) for INITIAL. This corresponds to a retreat of 90 km.
OceanA and OceanB
For the hydrography, we specified two model scenarios, OceanA and OceanB. In both, the potential temperature Θ increases with depth, while an increase in salinity S, ranging from 34.52 to 34.82, keeps the stratification stable with respect to the potential density. For OceanA, Θ linearly increases from –1.9˚C at the surface to –1.5˚C at the lowermost level. OceanB has a surface temperature of –1.8˚C and reaches –1.2˚C at the sea floor. This configuration prescribes a moderate ocean warming of 0.1–0.3˚C compared to OceanA. The vertical Θ,S-profiles are kept fixed at the open boundary, while inside the domain the ice-shelf/ocean interaction will change the 3-D Θ,S-pattern. Figure 6a shows a temperature transect along the ice-shelf front for OceanB. In the outflow region, cold water masses from the cavity deflect the isotherms downwards.
Fig. 6. ROMBAX, OceanB. (a) Temperature transect along the ice-shelf front from south to north after 3 model years, looking from the ocean into the ice-shelf cavity. Water masses of higher temperature at the left flowing into the cavity provide heat for basal melting. In the outflow region, colder water deflects the isotherms downward. (b) The stream function represents the vertically integrated mass transport. A maximum of 4 Sv (1 Sv = 106m3 s–1) is prescribed at the open-ocean boundary. Within the ice-shelf cavity, the transport is <0.5 Sv. The ice-shelf front is situated at ~5˚ W.
OceanA and OceanB melt rates calculated after 3 years of integration are shown in Figure 7. In both cases, the highest melt rates occur near the ice-shelf front where the main fraction of the ocean current is situated. This is due to the dynamic principle that ocean currents follow f/H w (f is the Coriolis parameter, H w the water column thickness) contours, meaning that, in our model set-up, only ~10% of the water masses enter the cavity. This value could rise if the relationship of ice thickness to water-column height decreased to a condition similar to that found, for instance, in the Filchner Trough. Compared to the real bathymetry in that area, our model geometry imposes an extra barrier for an exchange of water masses (Fig. 6b).
Fig. 7. Melt rates as derived by ROMBAX for (a) OceanA and (b) OceanB. Average melt rates are 0.35 and 0.55 ma–1, respectively. Highest values occur at the ice-shelf front, for OceanB.
Due to the cyclonic circulation, inflow into the cavity occurs in the southern region. Although the main current evolves along the ice-shelf front, a portion of the water masses propagates along the grounding line. Here the potential for melting is highest, but the small transport carries only enough heat to cause melt rates of 0.5 ma–1 near the grounding line compared to 6.0ma–1 at the ice-shelf front. OceanB yields respective values of 1.0 and 12 ma–1 at some front regions. Average melt rates amount to 0.35 ma–1 for OceanA and 0.55 ma–1 for OceanB.
As in experiments Melt05 and Melt20, we use the output from the ocean model to solve the continuity Eqn (3) in OceanA and OceanB. We simulate the ice thickness evolution for 1000 years, again starting from the control run INITIAL at 16.000 years. Melt rates from Figure 7 last throughout the perturbation time-span. In the case of a retreating ice sheet, a simple extrapolation of melt rates towards the new grounding line handles geometric changes of the cavity. No attempt is made to apply a full coupling of ice and ocean models including repeated reruns of ROMBAX.
As would be expected, H decreases more in the ocean inflow area (y = 15 (=75 km)) than in the outflow area (y=85 (=425 km)). This is the case in both experiments, but the higher melt rates of OceanB make the effect more obvious. The maximum decrease in ice thickness here is in the order of 100 m. Figure 8b provides a 2-D view of the changes as a comparison of Melt20 and OceanB. Decreases exceed 300 m within the ice shelf close to the grounding line in Melt20.
Fig. 8. (a) Resulting thickness profiles for OceanB at T = 17.000 years. High melt rates in the inflow region (green) compared with the profile at T = 16.000 years. High melt rates at the ice-shelf front cause a thinning by up to 30% in the inflow region. Thinning stays moderate, where a colder outflow would allow even freezing (red). Due to symmetry, INITIAL is identical on both sides. (b, c) 2-D views of ice-thickness changes in metres: Melt20– INITIAL (b); OceanB – INITIAL, with different scale (c). The red line marks the boundary of the ice shelf at the end of the simulation.
Discussion and Outlook
With a view to quantifying the stability of the Antarctic ice sheet, we applied a coupled model to a synthetic ice-sheet/ ice-shelf ocean system. Special emphasis is given to the ocean and its contribution to the ice-sheet mass balance. High basal melt rates thereby codetermine the steady-state profile of our idealized model ice-sheet configuration.
From the hysteresis experiment (Fig. 2) we conclude that the loss of ice is faster after a grounding-line retreat than the rebuild following a grounding-line advance. Once ice becomes afloat, the corresponding grid nodes obey the ice-shelf flow, causing higher strain rates with respect to a given surface gradient than compared to an ice-sheet flow. This should also happen in nature, implying that the ice shelf would drain inland ice faster than before. On the other hand, if grounding leads to a change from ice-shelf to ice-sheet flow, smaller velocities (fluxes) in the ice sheet need longer to build up an ice-sheet profile at the new margin.
In our study we investigated a coupled 3-D ice-sheet/ice-shelf ocean model. Although each model part underlies a different timescale to attain a steady state, our objectives concentrate on the time-dependent behaviour of ice shelves. Being forced by perturbations of the basal mass balance, ice shelves react within 103 years or less, depending on size and flow velocity. According to the results in Table 1, our system loses on average about 40 Gt ice a–1, in experiment Melt20. This value, which is even higher at the onset of the melt event, even illustrates a possible impact on a centennial timescale.
Running ROMBAX only once at the beginning of the simulations does not take into account changes in the ice-shelf geometry and their implications for the sub-ice circulation. If we expect a balance of ice dynamics and basal mass flux, erosion of the ice-shelf bottom will make the interface with the ocean steeper. Simultaneously, if the draft at the grounding line stays almost constant, the total magnitude of depth variation will increase. This results in a higher potential for circulation within the cavity, in line with increased melting. In the case of an inland sloping bedrock, a grounding-line retreat would foster melting even more. Therefore we conclude that our restriction to one ROMBAX run at the beginning of the perturbation time period will underestimate the potential for melting. Nevertheless, future experiments will have to use a perennial coupling of ice and ocean models, especially if climate scenarios are taken into account (e.g. Reference Grosfeld and SandhägerGrosfeld and Sandhäger, 2004).
Concerning the coupling to the grounded ice sheet, the main simplification in our model is the use of the SIA and SSA. Since these approximations neglect higher-order stress components, the transition from grounding to floating takes place at a single point, where the flotation criterion holds. No attempt is made to incorporate partly grounded areas like ice streams whose wet bases could explain their accelerated flow.
Being aware of our simplifications, we are confident that our results underestimate the effect of melting ice on the volume change in the grounded parts. By omitting an ice-stream process, our model restrains the discharge of inland ice when a shrinking ice shelf leads to higher thickness gradients at the grounding line. Nevertheless, experiment Melt20 shows that a substantial ice-volume loss could be possible if our ice shelf were exposed to a sufficiently high melt rate.
Still, OceanA and OceanB do not simulate the case of an intrusion of warm water masses onto the continental shelf like those reported for the eastern Amundsen Sea by Reference JenkinsJenkins and others (2010). Here water as warm as +1˚C comes into contact with the floating glacier, causing melt rates higher than 30ma–1 near the grounding line (Reference Payne, Holland, Shepherd, Rutt, Jenkins and JoughinPayne and others, 2007; Reference Jacobs, Jenkins, Giulivi and DutrieuxJacobs and others, 2011).
On the other side of the Antarctic Peninsula, the coastal current reaches temperatures of +0.75˚C at the continental-shelf break off Kapp Norvegia (Reference Fahrbach, Hoppema, Rohardt, Schröder and WisotzkiFahrbach and others, 2004). If these water masses at intermediate depth flushed the Filchner Trough and partially the cavity of the Filchner Ice Shelf, they could also cause high melting.
Our ocean model experiment OceanA yields a mean melt rate of 0.35ma–1 which is in good agreement with 0.32 ma–1 from Reference HellmerHellmer (2004) for the Filchner Ice Shelf. Rates of this magnitude seem not to harm our model ice sheet, but if water masses of higher temperature entered the cavity, average melt rates as high as in scenario Melt20 would be possible. Reference HellmerHellmer (2004) presented such present-day values for the Getz Ice Shelf, as well as for the ice shelves of the eastern Weddell Sea. Our results imply that if this is the case for large ice shelves like Filchner and Ronne, the anticipated ice-volume loss in grounded areas would contribute substantially to sea-level rises. This might happen even on timescales less than 1000 years.
Melt-rate outcomes of ROMBAX show a large variability in magnitude and horizontal distribution. This even holds for our regularly shaped model domain. If applied as a boundary condition to an ice-shelf base, they will influence the ice thickness in both horizontal directions. This implies that simple parameterizations of basal melt rates (e.g. with respect to offshore ocean temperatures) could lead to misinterpreted ice-thickness distributions.
The model results shown demonstrate the possible suitability of the coupled model for application with real topography. Hence, favoured regions of the Antarctic continent can be simulated now, with a special emphasis on their vulnerability to increasing ocean temperatures. To quantify such scenarios in further studies, we aim to apply our coupled model to the Filchner Ice Shelf in the Weddell Sea, and Pine Island Glacier in the Amundsen Sea, including their inland drainage basins.
Acknowledgements
This work was supported by the ice2sea programme, funded through the European Union 7th Framework Programme, grant No. 226375. This is ice2sea contribution No. 044. We gratefully acknowledge the two anonymous reviewers.
