Numerical Ice Sheet Models Research Articles

Investigating the stability of subglacial lakes with a full Stokes ice-sheet model

AbstractDespite the large amount of subglacial lakes present underneath the East Antarctic ice sheet and the melt processes involved, the hydrology beneath the ice sheet is poorly understood. Changes in subglacial potential gradients may lead to subglacial lake outbursts, discharging excess water through a subglacial drainage system underneath the ice sheet. Such processes can eventually lead to an increase in ice flow. In this paper, a full Stokes numerical ice-sheet model was employed which takes into account the ice flow over subglacial water bodies in hydrostatic equilibrium with the overlying ice. Sensitivity experiments were carried out for small perturbations in ice flow and basal melt rate as a function of ice thickness, general surface slope, ice viscosity and lake size, in order to investigate their influence on the subglacial potential gradient and the impact on subglacial lake drainage. Experiments clearly demonstrate that small changes in surface slope are sufficient to start and sustain episodic subglacial drainage events. Lake drainage can therefore be regarded as a common feature of the subglacial hydrological system and may influence, to a large extent, the present and future behavior of large ice sheets.

An improved numerical scheme to compute horizontal gradients at the ice-sheet margin: its effect on the simulated ice thickness and temperature

AbstractIn three-dimensional numerical ice-sheet models that use finite-difference schemes, the position of ice margins is poorly represented with a regular quadratic grid. As a result, in a centered difference scheme, the surface gradient term and the flux divergence term computed for the gridpoints next to the ice margin may be inaccurate. In this paper, an improved scheme is presented that computes the horizontal gradients at the ice-sheet margin using an asymmetric (upstream) second-order difference scheme in order to avoid using information from the zero-thickness gridpoints. The model is applied to an idealized synthetic geometry to obtain a steady-state ice-sheet topography. The improved model shows a realistically smooth thickness distribution near the margin. Thermomechanical coupling is found to enhance the error near the margin. The error in simulated thicknesses with the centered-difference method was significantly reduced with the new upstream scheme.

Towards a GIS assessment of numerical ice-sheet model performance using geomorphological data

AbstractA major difficulty in assimilating geomorphological information with ice-sheet models is the lack of a consistent methodology to systematically compare model output and field data. As an initial step in establishing a quantitative comparison methodology, automated proximity and conformity analysis (APCA) and automated flow direction analysis (AFDA) have been developed to assess the level of correspondence between modelled ice extent and ice-marginal features such as end moraines, as well as between modelled basal flow directions and palaeo-flow direction indicators, such as glacial lineations. To illustrate the potential of such an approach, an ensemble suite of 40 numerical simulations of the Fennoscandian ice sheet were compared to end moraines of the Last Glacial Maximum and the Younger Dryas and to glacial lineations in northern Sweden using APCA and AFDA. Model experiments evaluated in this manner were ranked according to level of correspondence. Such an approach holds considerable promise for optimizing the parameter space and coherence of ice-flow models by automated, quantitative assessment of multiple ensemble experiments against a database of geological or glaciological evidence.

Glaciology of the British Isles Ice Sheet during the last glacial cycle: form, flow, streams and lobes

The last ice sheet over the British Isles, together with other mid-latitude Pleistocene ice sheets, and in contrast to the modern ice sheets of Greenland and Antarctica, had a relatively low profile, low summit elevation and extensive, elongated lobes at its margin. A thermo-mechanically coupled numerical ice sheet model, driven by a proxy climate, has been used to explore the properties that would have permitted these characteristics to develop. The approach, the key to quantitative palaeoglaciology, is to determine the boundary conditions that permit the simulated ice sheet to mimic the evolution of the real ice sheet through the last glacial cycle. Simulations show how a British ice sheet may have been confluent with a Scandinavian ice sheet during some parts of its history and how unforced periodic and asynchronous oscillations could occur in different parts of its margins. Marginal lobes are a reflection of streaming within the ice sheet. Such streams can be ephemeral, dynamic streams located because of ice sheet properties, or fixed streams whose location is determined by bed properties. The simulations that best satisfy constraints of extent, elevation and relative sea levels are those with major fixed streams that strongly draw down the ice sheet surface. In these, the core upland areas of the ice sheet were cold based at the Last Glacial Maximum, basal streaming velocities were between 500 and 1000 ma −1 compared with surface velocities of 10–50 ma −1 in inter-stream zones, shear stresses were as low as 15–25 kPa in streams compared with 70–110 kPa in upland areas and 60–84% of the ice flux was delivered to the margin via streams.

Role of transition zones in marine ice sheet dynamics

It has been suggested that rapid thinning of the West Antarctic ice sheet is due to increased melting under ice shelves caused by a gradual ocean warming (Shepherd et al., 2004). Payne et al. (2004) showed that such melting could lead to an acceleration of grounded ice flow. In this paper, we analyze the response of a marine ice sheet to different perturbations near the grounding line using a numerical ice sheet model that takes into account longitudinal stress coupling and grounding line migration at subgrid precision, based on a novel technique. Results show that stress transmission or longitudinal coupling across the grounding line plays a decisive role. The grounding line migration is a function of the length scale over which the basal conditions change from frozen to the bed to floating, the “transition zone.” We demonstrate that thinning of the ice shelf due to bottom melting has a negligible effect on the grounded ice mass. Only perturbations at the grounding line or reduction in buttressing of the ice shelf substantially thins the grounded ice sheet. Inclusion of lateral drag does not alter these results qualitatively. Marine ice sheets with large transition zones, such as ice streams, seem highly sensitive to such perturbations compared to ice sheets with small transition zones, such as an abrupt ice sheet/ice shelf junction.

Modelling iceberg trajectories, sedimentation rates and meltwater input to the ocean from the Eurasian Ice Sheet at the Last Glacial Maximum

We model iceberg flow paths from the Eurasian Ice Sheet, and the associated meltwater production and sedimentation rates within the Norwegian–Greenland Sea during the last glaciation. Results from a numerical ice sheet model, an atmospheric general circulation model and an ocean general circulation model are collated and used to provide iceberg production rate, wind field and surface current forcings to an iceberg trajectory model. The iceberg model then determines how icebergs issuing from the Eurasian Ice Sheet travel across the ocean and eventually melt. In addition the release of iceberg sediments is also predicted. The results show that iceberg trajectories are complex and that common features of iceberg movement are clustering in zones of convergence and exit into the North Atlantic through the Iceland–Faeroes Channel. Eurasian icebergs do not penetrate into the interior of the Arctic Ocean. The gathering of icebergs produces a complex meltwater pattern that does not follow the conceptual idea of decreasing meltwater production with distance from the ice sheet margin. Sedimentation results are compared with the meltwater configuration and are found to be a poor indicator of past zones of iceberg melt, with zones of sedimentation extending significantly less far.

GRANTISM: An Excel TM model for Greenland and Antarctic ice-sheet response to climate changes

Over the last decades, the response of large ice sheets on Earth, such as the Greenland and Antarctic ice sheets, to changes in climate has been successfully simulated with large-scale numerical ice-sheet models. Since these models are highly sophisticated, they are only applicable on the scientific level as they demand a large amount of CPU time. Based on similar physics, a computationally fast flowline model of the Greenland and Antarctic ice sheet is presented here, primarily designed for educational purposes. Using an over-implicit numerical scheme, the model runs fast and behaves in a similar way to changes in background temperature forcing as major ice-sheet models do. A user-friendly interface and the implementation within a common spreadsheet program (Excel TM) make the model suitable for the classroom.

Comparing predicted and observed spatial boundaries of geologic phenomena: Automated Proximity and Conformity Analysis applied to ice sheet reconstructions

Comparing predicted with observed geologic data is a central element of many aspects of research in the geosciences, e.g., comparing numerical ice sheet models with geomorphic data to test ice sheet model parameters and accuracy. However, the ability to verify predictions using empirical data has been limited by the lack of objective techniques that provide systematic comparison and statistical assessment of the goodness of correspondence between predictions of spatial and temporal patterns of geologic phenomena and the field evidence. Much of this problem arises from the inability to quantify the level of agreement between straight or curvilinear features, such as between the modeled extent of some geologic phenomenon and the field evidence for the extent of the phenomenon. Automated Proximity and Conformity Analysis (APCA) addresses this challenge using a system of Geographic Information System-based buffering that determines the general proximity and parallel conformity between linear features. APCA results indicate which modeled output fits empirical data, based on the distance and angle between features. As a result, various model outputs can be sorted according to overall level of agreement by comparison with one or multiple features from field evidence, based on proximity and conformity values. In an example application drawn from glacial geomorphology, APCA is integrated into an overall model verification process that includes matching modeled ice sheets to known marginal positions and ice flow directions, among other parameters. APCA is not limited to ice sheet or glacier models, but can be applied to many geoscience areas where the extent or geometry of modeled results need to be compared against field observations, such as debris flows, tsunami run-out, lava flows, or flood extents.

Assessing the ability of numerical ice sheet models to simulate grounding line migration

Grounding line migration is a key process affecting the stability of marine ice sheets such as the West Antarctic ice sheet (WAIS). Grounding line motion is often included in numerical models simulating the past and future evolution of the WAIS; however, little attention has been paid to the numerical consistency of these models. The aim of this paper is to assess the ability of simple versions of existing marine ice sheet models to simulate grounding line migration. In particular, we investigate the response of the grounding line to external forcing and the sensitivity of the models' predictions to their numerics and the mechanical coupling between ice sheet and shelf. From the model comparison, there is no consensus on how the grounding line should react to changes in boundary conditions. A crucial finding is the strong dependency of models using a fixed grid on numerical details such as the horizontal grid size. This implies that we should be very wary about grounding line predictions from such models. Including mechanical coupling at the grounding line does not seem to change the qualitative behavior of the models. This suggests that the way the grounding line is treated in marine ice sheet models dominantly determines the grounding line dynamics. We find that models that employ a moving grid to explicitly track the grounding line do not share many of the deficiencies of the fixed grid models. We conclude that at present, no reliable model of the grounding line is available, and further model development is urgently needed.

Spectral roughness of subglacial topography and implications for former ice-sheet dynamics in East Antarctica

Ice-sheet basal erosion is controlled by ice dynamics (including the basal thermal regime) and the lithology of the substrate. Spatial variation in subglacial roughness is therefore likely to be a function of these controls. In Antarctica, very little is known about former ice dynamics and sub-ice geology. Here, we calculate the spectral roughness of subglacial East Antarctica from an analysis of radio-echo sounding data. As the modern glaciological setting is understood reasonably well in East Antarctica from numerical ice-sheet modelling, we are able to compare the roughness calculations with contemporary ice-sheet dynamics. We show that ice divides at Ridge B and Dome A are underlain by rough terrain and a cold thermal regime. At the Dome C ice divide, however, the bed is noticeably smooth at the shorter wavelengths (∼5–10 km) but rougher at longer wavelengths. One interpretation is that subglacial morphology has been eroded, leaving a smoothed landscape that contains valleys and hills. Regardless of the origin of the valleys, their erosion is only possible if the ice sheet were more dynamic than at present. Hence, whereas Dome A and Ridge B have a bed morphology consistent with the present ice sheet, the morphology at Dome C is likely to predate the present ice-sheet configuration.

Fast-flowing outlet glaciers of the Last Glacial Maximum Patagonian Icefield

Glacial geomorphology around the Northern Patagonian Icefield indicates that a number of fast-flowing outlet glaciers (the continuation of ice streams further upglacier) drained the icefield during the Last Glacial Maximum. These topographically controlled fast-flowing glaciers may have dictated the overall pattern of Last Glacial Maximum ice discharge, lowered the ice-surface profile, and forced the ice-divide westward. The influence of the fast-flowing outlet glaciers on icefield behavior also helps to explain why the configuration of the Patagonian Icefield at the Last Glacial Maximum is not accurately represented in existing numerical ice-sheet models. Fast-flowing outlet glaciers would have strongly influenced ice discharge patterns and therefore partially decoupled the icefield from climatically induced changes in thickness and extent.

Late Quaternary ice sheet history of northern Eurasia

The maximum limits of the Eurasian ice sheets during four glaciations have been reconstructed: (1) the Late Saalian (>140 ka), (2) the Early Weichselian (100–80 ka), (3) the Middle Weichselian (60–50 ka) and (4) the Late Weichselian (25–15 ka). The reconstructed ice limits are based on satellite data and aerial photographs combined with geological field investigations in Russia and Siberia, and with marine seismic- and sediment core data. The Barents-Kara Ice Sheet got progressively smaller during each glaciation, whereas the dimensions of the Scandinavian Ice Sheet increased. During the last Ice Age the Barents-Kara Ice Sheet attained its maximum size as early as 90–80,000 years ago when the ice front reached far onto the continent. A regrowth of the ice sheets occurred during the early Middle Weichselian, culminating about 60–50,000 years ago. During the Late Weichselian the Barents-Kara Ice Sheet did not reach the mainland east of the Kanin Peninsula, with the exception of the NW fringe of Taimyr. A numerical ice-sheet model, forced by global sea level and solar changes, was run through the full Weichselian glacial cycle. The modeling results are roughly compatible with the geological record of ice growth, but the model underpredicts the glaciations in the Eurasian Arctic during the Early and Middle Weichselian. One reason for this is that the climate in the Eurasian Arctic was not as dry then as during the Late Weichselian glacial maximum.

Late Quaternary ice sheet history of northern Eurasia

The maximum limits of the Eurasian ice sheets during four glaciations have been reconstructed: (1) the Late Saalian (>140 ka), (2) the Early Weichselian (100–80 ka), (3) the Middle Weichselian (60–50 ka) and (4) the Late Weichselian (25–15 ka). The reconstructed ice limits are based on satellite data and aerial photographs combined with geological field investigations in Russia and Siberia, and with marine seismic- and sediment core data. The Barents-Kara Ice Sheet got progressively smaller during each glaciation, whereas the dimensions of the Scandinavian Ice Sheet increased. During the last Ice Age the Barents-Kara Ice Sheet attained its maximum size as early as 90–80,000 years ago when the ice front reached far onto the continent. A regrowth of the ice sheets occurred during the early Middle Weichselian, culminating about 60–50,000 years ago. During the Late Weichselian the Barents-Kara Ice Sheet did not reach the mainland east of the Kanin Peninsula, with the exception of the NW fringe of Taimyr. A numerical ice-sheet model, forced by global sea level and solar changes, was run through the full Weichselian glacial cycle. The modeling results are roughly compatible with the geological record of ice growth, but the model underpredicts the glaciations in the Eurasian Arctic during the Early and Middle Weichselian. One reason for this is that the climate in the Eurasian Arctic was not as dry then as during the Late Weichselian glacial maximum.

The effect of basal topography on ice sheet dynamics

Classical shallow-ice theory assumes that bed topography under ice sheets has slopes comparable to the surface slope of the ice sheet. A modification of the classical steady-state theory which allows for significant bed topography on shorter length scales has recently been developed by Morland (Proc. R. Soc L. Ser. A., 456, 1711-1739), but his theory requires explicit integration of the ice-flow equations over the topography length scale, which may be below the grid size of typical numerical ice sheet models. Here we present a method for parameterising the effect of basal topography of wavelengths much greater than ice thickness but much smaller than the horizontal extent of the ice sheet on the bulk flow of the ice sheet. In particular, we are able to show through the use of a multiple-scales expansion technique that the effect of such topography is described by a simple correction factor applied to the classical expression for ice flux. This correction factor dispenses with the need to integrate explicitly over the topography length scale and could allow the effect of such topography to be included in numerical models with limited grid size. Examples are given for the practical implementation of this 'correction factor method' in calculations of the steady-state shape of ice sheets.

Modelling the southern extent of the last Icelandic ice‐sheet

AbstractA three‐dimensional, time‐dependent numerical ice‐sheet model is used to investigate the southern extent of the last Icelandic ice‐sheet since the Last Glacial Maximum (LGM, ca. 20 000 14C yr BP). Ice development over southern Iceland is forced using a linear relationship between mass balance and altitude based on observations over Sólheimajökull, a southern outlet glacier of the Mýrdalsjökull ice‐cap. A continentality factor is introduced that raises inland equilibrium‐line altitudes (ELAs) and also slackens the mass‐balance/altitude gradient driving the model. With this calibration the present‐day ice distribution can be reconstructed. Growth from ice‐free to total ice cover is assessed and the sensitivity of ice extent to ELA change is shown to be non‐linear. The model indicates that an ELA lowering of 500 m, consistent with a ca. 5°C temperature depression for southern Iceland at the LGM, would enable glacier ice to cover the whole land surface within the modelled area and inundate any putative ecological refugia. Modelling with a 300 m ELA depression, consistent with previous reconstructions of Younger Dryas (YD, 10 000–11 000 14C yr BP) ice extent, indicates that the principal outcrops of the Sólheimar ignimbrite were glaciated at this time, suggesting that an origin during this stadial (and a correlation with North Atlantic Ash Zone I / Vedde Ash) is problematic. Copyright © 2003 John Wiley & Sons, Ltd.

A comparison of Eulerian and Lagrangian methods for dating in numerical ice-sheet models

AbstractAccurate dating in ice sheets is required for a correct interpretation of palaeoclimatic records and for incorporation of material characteristics in the flow law which depend on ice age. In this paper, we make a comparison between a Lagrangian and Eulerian approach to the ice advection problem in numerical ice-sheet models. This comparison is first performed for a schematic two-dimensional ice sheet of Nye–Vialov type with a prescribed stationary velocity field. Severalcases are examined which incorporate basal melting, basal sliding and an undulating bed. A further comparison is made with an analytical solution for the ice divide. Both methods are also applied in a thermomechanical model of the Antarctic ice sheet for steady-state present-day conditions. Our main conclusion is that, for similar discretization parameters, the Lagrangian method produces less error than an Eulerian approach using a second-order upwinding finite-difference scheme, though the difference is small (<1%) for the largest part of the model domain. However, problems with the Lagrangian approach are introduced by the dispersion of tracers, necessitating the use of interpolation procedures that are a main source of additional error. It is also shown that a cubic-spline approximation of Lagrangian trajectories improves accuracy, but such a method is computationally hardly applicable in large-scale ice-sheet models.

Streaming flow in an ice sheet through a glacial cycle

AbstractGeological evidence indicates that the flow of the last European ice sheet was dominated by numerous large ice streams. Although some were ephemeral, most were sustained along well-defined axes at least during the period of retreat after the Last Glacial Maximum. A thermomechanically coupled three-dimensional numerical ice-sheet model has been used to simulate the ice sheet through the whole of the last glacial cycle, but with a spatial resolution that is high enough to capture streaming behaviour. An experiment with a smoothed bed is used to explore the self-organizing behaviour of streams when they are not forced by bed topography. On such a bed, streams typically have a width of 1–10 km, much narrower than the inferred European ice streams. An experiment using a realistic topography suggests that widths of ice streams are strongly influenced by topography, and tend to be of order 100 km. Moreover, even where the topography is muted, it stabilizes the locations of ice streams which, once formed, tend to be sustained along pre-existing axes. The model creates patterns of streaming that are similar to inferred patterns, suggesting strong topographic forcing. In a simulation using a realistic bed in which the ice was very cold and basal melting rarely occurred, streams were again very narrow. Widespread streaming under low driving stresses tends to reduce ice-sheet thicknesses compared with weak streaming or models that do not produce streaming. Consequently, ice thicknesses are smaller and tend to be consistent with the results of sea-level inversions based on geophysical Earth models.

Late Weichselian iceberg, surface-melt and sediment production from the Eurasian Ice Sheet: results from numerical ice-sheet modelling

Results from a numerical ice-sheet model were matched with geological evidence detailing the extent and timing of the Late Weichselian ice sheet in the Eurasian Arctic, through simple adjustments in the model’s environmental inputs. The numerical model predicts the spatial and temporal variation in (1) subglacial sediment transport and deposition, (2) rates of iceberg calving and (3) rates of supraglacial melting over the past 30 000 yr. The ice sheet is characterised by a series of ice streams occupying bathymetric troughs in the Barents Sea and west of Norway, which act as sources for glacial sediment build-up and iceberg production. Our model indicates that significant volumes of sediment are deposited over the Bear Island fan (3000 km 3 in 12 000 yr), Franz Victoria Trough fan (500 km 3 in 8000 yr) and a series of smaller fans west of Norway (which combine to yield 4000 km 3 in 12 000 yr). The Norwegian Channel ice stream operated for only 4000 yr, during which time 400 km 3 of sediments were deposited on the adjacent fan. The ice sheet experienced two major periods of iceberg production where the rate of calving increased by ∼50% at about 15 000 and 12 500 yr ago. However, the time-dependent response of iceberg calving to sea-level rise for individual ice streams is shown to be more complicated than the generalised ice-sheet response. The marine portions of the ice sheet decayed after 15 000 yr ago, resulting in several embayments at the mouths of bathymetric troughs. A second pulse of enhanced iceberg calving at 12 500 yr ago caused further decay of the marine ice sheet and the ice margin retreated back to the shorelines of island archipelagos in the northern Barents Sea. Ice loss through surface melting was restricted mostly to the southern margin of the ice sheet.

The timing of initiation of fast-flowing ice streams during a glacial cycle inferred from glacimarine sedimentation

Fast-flowing ice streams drain huge basins within modern ice sheets and are a mechanism for rapid iceberg production and mass loss. Some ice streams are known to be unstable, and to switch between fast flow and stagnation. Numerical ice-sheet models and geophysical observations of large-scale streamlined glacial landforms and large glacier-derived sedimentary fans along high-latitude continental margins indicate the presence of ice streams in Quaternary ice sheets. High-resolution AMS dated sedimentary records from an inter-fan area on the Svalbard continental margin yield evidence on the timing of initiation of fast glacier flow during a glacial cycle of ice-sheet growth and decay. Sedimentation rates and iceberg-rafted debris influx show that debris delivery to the Svalbard margin was relatively high during the period of Late Weichselian ice-sheet growth from about 30 000–18 000 yr BP. Both parameters fell dramatically in inter-fan areas when full-glacial conditions, with ice at the continental shelf edge, were established at about 18 000 yr BP. This shift is interpreted to mark the initiation of a fast-flowing ice stream, the convergence of flow from the surrounding ice sheet into this ice stream, and the consequent slowing of ice flow and sediment delivery to the inter-fan area beyond the ice-stream margins. This provides a constraint on the timing of ice-stream initiation during the cycle of ice-sheet growth and decay.

On glacial inception and orography

There is no doubt that orography plays a crucial role in the initiation of ice sheets. Mass balance observations on glaciers in high-mountain and polar climates show that the specific balance varies first of all with altitude. The height–mass balance feedback (HMB-feedback) is the most effective process to turn local ice caps into big ice sheets. Here, a study is presented on the processes that determine the strength of the HMB-feedback. This is done by carrying out a set of integrations with a numerical ice-sheet model for a bounded domain and a sinusoidally varying bed shape. Results are analysed for different wavenumbers and amplitudes. It turns out that the relative importance of the geometric effect (higher mountains allow earlier glacierisation) and the ice-mechanical effect (dependence of ice flow on the bed slope) varies significantly with wavenumber and amplitude. With regard to full glaciation of the domain, this implies that for a given wavenumber there is a preferred amplitude (that is, the drop in the equilibrium-line altitude needed to initiate full glaciation has a minimum).