Three-dimensional Thermomechanical Ice-sheet Model Research Articles

The North American Cordillera—An Impediment to Growing the Continent-Wide Laurentide Ice Sheet

Abstract This study examines the evolution of a continental-scale ice sheet on a triangular representation of North America, with and without the influence of the Cordilleran region. Simulations are conducted using a comprehensive atmospheric general circulation model asynchronously coupled to a three-dimensional thermomechanical ice-sheet model. The atmospheric state is updated for every 2 × 106 km3 increase in ice volume, and the coupled model is integrated to steady state. In the first experiment a flat continent with no background topography is used. The ice sheet evolves fairly zonally symmetric, and the equilibrium state is continent-wide and has the highest point in the center of the continent. This equilibrium ice sheet forces an anticyclonic circulation that results in relatively warmer (cooler) summer surface temperatures in the northwest (southeast), owing to warm (cold) air advection and radiative heating due to reduced cloudiness. The second experiment includes a simplified representation of the Cordilleran region. The ice sheet’s equilibrium state is here structurally different from the flat continent case; the center of mass is strongly shifted to the east and the interior of the continent remains ice free—an outline broadly resembling the geologically determined ice margin in Marine Isotope Stage 4. The limited glaciation in the continental interior is the result of warm summer surface temperatures primarily due to stationary waves and radiative feedbacks.

The Greenland Ice Sheet at the peak of warming during the previous Interglacial

The Last Interglacial (LIG or the Eemian) between ca. 130 and 115 kyr BP is probably the best analogue for future climate warming for which increasingly better proxy data are becoming available. The volume of the Greenland Ice Sheet (GrIS) during this period is of particular interest to better assess how much and how fast sea-level can rise in a future Earth undergoing gradual climatic warming. Sea-level during the LIG is inferred to have been up to 9 meter higher than today, but contribution of the GrIS into this rise remains unclear. Various ice-sheet modeling studies have come up with a very broad range of the LIG volume loss by the GrIS to between 60 cm and 6 m of equivalent sea-level rise. This wide range is explained by the sensitivity of GrIS models to the imposed climatic conditions and to poor knowledge of the LIG climate itself in terms of the magnitude and precise timing of the maximum warming, as well as in terms of spatial and annual patterns. To partially circumvent these uncertainties we made use of the newest temperature record over the Central Greenland reconstructed from the isotopic composition of the recently obtained NEEM ice core containing undisturbed LIG segment to build the climatic forcing of the model. The NEEM record unequivocally indicates times of the start and of the end of the LIG warming in Greenland as well as the duration of the warmest time period within the Eemian. Using a three-dimensional thermomechanical ice-sheet model, we produced an ensemble of possible LIG configurations by varying only four key parameters for temperature, precipitation rate, surface melting magnitude and melting pattern within realistic bounds. The outcome of a series of the numerical experiments is a variety of glaciologically consistent GrIS geometries corresponding to a wide range of possible «climates». To constrain the ensemble of GrIS geometries, we used data inferred from 5 Greenland ice cores such as the presence or absence of LIG ice, borehole temperature and isotopic composition. Lagrangian backtracing of particles was used to calculate non-climatic biases in isotopic records introduced by horizontal advection, systematic latitudinal contrast and local elevation changes. Comparison of model-generated ice-core characteristics with the observed data enabled to narrow down the ensemble to a bound on the GrIS contribution to the LIG sea-level rise of between 1.3 and 2.9 m with the best choice of 1.8–2.2 m. This conclusion in general supports the point of view about the modest GrIS contribution to global sea level rise during the LIG.

Interactions between topographically and thermally forced stationary waves: implications for ice-sheet evolution

ABSTRACTThis study examines mutual interactions between stationary waves and ice sheets using a dry atmospheric primitive-equation model coupled to a three-dimensional thermomechanical ice-sheet model. The emphasis is on how non-linear interactions between thermal and topographical forcing of the stationary waves influence the ice-sheet evolution by changing the ablation. Simulations are conducted in which a small ice cap, on an idealised Northern Hemisphere continent, evolves to an equilibrium continental-scale ice sheet. In the absence of stationary waves, the equilibrium ice sheet arrives at symmetric shape with a zonal equatorward margin. In isolation, the topographically induced stationary waves have essentially no impact on the equilibrium features of the ice sheet. The reason is that the temperature anomalies are located far from the equatorward ice margin. When forcing due to thermal cooling is added to the topographical forcing, thermally induced perturbation winds amplify the topographically induced stationary-wave response, which that serves to increase both the equatorward extent and the volume of the ice sheet. Roughly, a 10% increase in the ice volume is reported here. Hence, the present study suggests that the topographically induced stationary-wave response can be substantially enhanced by the high albedo of ice sheets.

Greenland ice-sheet volume sensitivity to basal, surface and initial conditions derived from an adjoint model

AbstractWe extend the application of control methods to a comprehensive three-dimensional thermomechanical ice-sheet model, SICOPOLIS (SImulation COde for POLythermal Ice Sheets). Lagrange multipliers, i.e. sensitivities, are computed with an exact, efficient adjoint model that has been generated from SICOPOLIS by rigorous application of automatic differentiation. The case study uses the adjoint model to determine the sensitivity of the total Greenland ice volume to various control variables over a 100 year period. The control space has of the order 1.2 × 106 elements, consisting of spatial fields of basal flow parameters, surface and basal forcings and initial conditions. Reliability of the adjoint model was tested through finite-difference perturbation calculations for various control variables and perturbation regions, ascertaining quantitative inferences of the adjoint model. As well as confirming qualitative aspects of ice-sheet sensitivities (e.g. expected regional variations), we detect regions where model sensitivities are seemingly unexpected or counter-intuitive, albeit ‘real’ in the sense of actual model behavior. An example is inferred regions where sensitivities of ice-sheet volume to basal sliding coefficient are positive, i.e. where a local increase in basal sliding parameter increases the ice-sheet volume. Similarly, positive (generally negative) ice temperature sensitivities in certain parts of the ice sheet are found, the detection of which seems highly unlikely if only conventional perturbation experiments had been used. The object of this paper is largely a proof of concept. Available adjoint-code generation tools now open up a variety of novel model applications, notably with regard to sensitivity and uncertainty analyses and ice-sheet state estimation or data assimilation.

High-resolution numerical simulation of Younger Dryas glaciation in Scotland

We use a 500 m resolution three-dimensional thermomechanical ice-sheet model forced by a scaled GRIP temperature pattern to retrodict the extent of glaciers during the Younger Dryas episode in Scotland. Using empirical data from sources spanning half a century we systematically perturb temperature depression, precipitation distribution, and the amount of basal sliding to identify the parameter space that most closely reproduces the glacier margins identified from field investigations. Arithmetic comparison of predicted ice cover with empirical glacier extent enables mismatch to be quantified and an ‘optimum-fit’ timeslice to be identified. This ‘best-fit’ scenario occurs with a maximum mean annual temperature depression from present of 10 ∘ C and steep eastward and northward gradients imposed on a modern precipitation distribution, coupled with restricted basal sliding. Even small deviations around these values produce considerably less well-fitting glacier configurations, suggesting that only a narrow range of optimal parameterisation exists. Mismatch between modelled and empirically reconstructed glacier extents occurs as a consequence of local conditions not accommodated by the model, such as wind-blown snow accumulation and lake (loch) bathymetry, and where geological evidence is equivocal. At the domain scale, however, our simulation suggests that inception of Scottish glaciers occurs rapidly, and leads to a coherent ice cap within 400 years of initial climatic cooling. According to our model, the ice cap begins to decay relatively early in the stadial, probably as a result of increasing aridity leading to net thinning and recession of many of its margins, with final and catastrophic collapse of the ice cap occurring largely within a century.

Numerical modeling investigations of the subglacial conditions of the southern Laurentide ice sheet

AbstractSub- and proglacial bed conditions influence advance and retreat of an ice sheet. The existence and distribution of frozen ground is of major importance for better understanding of ice-flow dynamics and landform formation. The southern margin of the Laurentide ice sheet (LIS) was dominated by the presence of relatively thin ice lobes that seem to have been very sensitive to external and internal physical conditions. Their extent and dynamics were highly influenced by the interaction of subglacial and proglacial conditions. A three-dimensional thermomechanical ice-sheet model was coupled with a model for the thermal regime in the upper Earth crust. The model has been applied to the LIS in order to investigate the spatial distribution of thermal conditions at the bed. The evolution of the whole LIS was modeled for the last glacial cycle, with primary attention on correct reconstruction of the southern margin. Our results show extensive temporal and spatial frozen ground conditions. Only a slow degradation of permafrost under the ice was found. We conclude that there are significant interactions between the ice sheet and the underlying frozen ground and that these influence both ice dynamics and landform development.

Subglacial erosion and englacial sediment transport modelled for North American ice sheets

Reconstructions of North American glacial history have either been based on interpretation of the geologic record or on numerical models of continental ice dynamics. Because such models make no explicit reference to the geological products of glaciation, the full potential of a combined approach has not been realized. We present the results of a process-based model of sediment production, entrainment, transport and deposition that represents a first step toward uniting these two lines of investigation. In this model basal-ice processes are dependent on ice sheet geometry, hydrology and dynamics. A three-dimensional thermomechanical ice sheet model yields the ice geometry, flow and thermal conditions that are used by a model of subglacial hydrology; together these drive the basal process model. Comparisons of simulation results with mapped debris distributions of the Hudson Bay Paleozoic carbonates and of the distinctive Dubawnt Group are used to test the validity of the transport model. Dispersion of debris radially out of Hudson Bay is greatly overestimated by the model, and dispersion toward the bay equally underestimated. This is because the modelled ice sheet has a single dome over Hudson Bay at glacial maxima, including Last Glacial Maximum, rather than a saddle separating ice domes over Quebec-Labrador and Keewatin. Estimated rates of glacial erosion for the Cordilleran and Laurentide ice sheets are compared against those predicted by the model. The model underestimates these rates because it underestimates the vigour of the quarrying process.

Simulations of the Northern Hemisphere through the last Glacial-interglacial cycle with a vertically integrated and a three-dimensional thermomechanical ice-sheet model coupled to a climate model

We present simulations of the Northern Hemisphere land ice through the last glacial-interglacial cycle with a vertically integrated ice-sheet model and a three-dimensional thermomechanical ice-sheet model. Both models are coupled asynchronously to the zonally averaged Louvain-la-Neuve climate model, which includes simplified treatments of the atmosphere, ocean and sea ice. The two-dimensional vertically integrated ice-sheet model, which contains no thermomechanical coupling, was developed in spherical coordinates (Marsiat, 1994). The three-dimensional thermomechanical ice-sheet model was developed using the two-dimensional vertically integrated model as source. We compare results of the vertically integrated with those of the thermomechanical ice-sheet model. in the thermomechanical model the deformation properties of ice depend on the temperature within the ice and the enhancement factor; the latter is introduced to model, in a simplified approach, the different flow properties of Pleistocene and Holocene ice due to varying dust content. The computations with the thermomechanical model show that the growth and decay of the Northern Hemisphere ice sheets can be modelled with a common enhancement factor for all ice sheets. It is shown that there are model set-ups for the thermomechanical model yielding temporal developments of the total ice volume comparable to those of the vertically integrated model. Furthermore, we demonstrate that for the coupled climate/cryosphere system the total ice volume depends considerably on the enhancement factor.

Modelling of Last Glacial Maximum ice sheets using different accumulation parameterizations

We use a three-dimensional thermomechanical ice-sheet model, previously tested on the Greenland ice sheet, to reconstruct Last Glacial Maximum (LGM) ice sheets. We compare the effects on the results of the ice-sheet model of three different accumulation parameterization schemes. In the first and second schemes, LGM precipitation is computed from the present precipitation, taking and not taking into account moisture transport. In the third scheme, LGM precipitation and surface temperatures are computed using outputs of an atmospheric global circulation model (AGCM), treated in anomaly mode. Results are compared to the last reconstruction of the Northern Hemisphere ice sheets (Peltier, 1994), computed using global rebound rates in a visco-elastic model of the Earth’s crust. The first two accumulation parameterizations do not give satisfactory reconstructions of the LGM ice sheets, since they are unable to compute realistic LGM climatic conditions. The third method gives very satisfactory results, which leads us to conclude that the best way to obtain realistic LGM climatic conditions is to use AGCM outputs.

Modelling of Last Glacial Maximum ice sheets using different accumulation parameterizations

We use a three-dimensional thermomechanical ice-sheet model, previously tested on the Greenland ice sheet, to reconstruct Last Glacial Maximum (LGM) ice sheets. We compare the effects on the results of the ice-sheet model of three different accumulation parameterization schemes. In the first and second schemes, LGM precipitation is computed from the present precipitation, taking and not taking into account moisture transport. In the third scheme, LGM precipitation and surface temperatures are computed using outputs of an atmospheric global circulation model (AGCM), treated in anomaly mode.Results are compared to the last reconstruction of the Northern Hemisphere ice sheets (Peltier, 1994), computed using global rebound rates in a visco-elastic model of the Earth’s crust. The first two accumulation parameterizations do not give satisfactory reconstructions of the LGM ice sheets, since they are unable to compute realistic LGM climatic conditions. The third method gives very satisfactory results, which leads us to conclude that the best way to obtain realistic LGM climatic conditions is to use AGCM outputs.

A three-dimensional climate—ice-sheet model applied to the Last Glacial Maximum

A quasi-three-dimensional (3-D) climate model (Sellers, 1983) was used to simulate the climate of the Last Glacial Maximum (LGM) in order to provide climatic input for the modelling of the Northern Hemisphere ice sheets. The climate model is basically a coarse-gridded general circulation (GCM) with simplified dynamics, and was subject to appropriate boundary conditions for ice-sheet elevation, atmospheric CO2concentration and orbital parameters. When compared with the present-daysimulation, the simulated climate at the Last Glacial Maximum is characterized by a global annual cooling of 3.5°C and a reduction in global annualprecipitation of 7.5%, which agrees well with results from other, more complex GCMs. Also the patterns of temperature change compare fairly with mostother GCM results, except for a smaller cooling over the North Atlantic and the larger cooling predicted for the summer rather than for the winter over Eurasia.The climate model is able to simulate changes in Northern Hemisphere tropospheric circulation, yielding enhanced westerlies in the vicinity of the Laurentide and Eurasian ice sheets. However, the simulated precipitation patterns are less convincing, and show a distinct mean precipitation increase over the Laurentide ice sheet. Nevertheless, when using the mean-monthly fields of LGM minus present-day anomalies of temperature and precipitation rate to drive a three-dimensional thermomechanical ice-sheet model, it was demonstrated that within realistic bounds of the ice-flow and mass-balance parameters, veryreasonable reconstructions of the Last Glacial Maximum ice sheets could be obtained.

A three-dimensional climate—ice-sheet model applied to the Last Glacial Maximum

A quasi-three-dimensional (3-D) climate model (Sellers, 1983) was used to simulate the climate of the Last Glacial Maximum (LGM) in order to provide climatic input for the modelling of the Northern Hemisphere ice sheets. The climate model is basically a coarse-gridded general circulation (GCM) with simplified dynamics, and was subject to appropriate boundary conditions for ice-sheet elevation, atmospheric CO2 concentration and orbital parameters. When compared with the present-daysimulation, the simulated climate at the Last Glacial Maximum is characterized by a global annual cooling of 3.5°C and a reduction in global annualprecipitation of 7.5%, which agrees well with results from other, more complex GCMs. Also the patterns of temperature change compare fairly with mostother GCM results, except for a smaller cooling over the North Atlantic and the larger cooling predicted for the summer rather than for the winter over Eurasia.The climate model is able to simulate changes in Northern Hemisphere tropospheric circulation, yielding enhanced westerlies in the vicinity of the Laurentide and Eurasian ice sheets. However, the simulated precipitation patterns are less convincing, and show a distinct mean precipitation increase over the Laurentide ice sheet. Nevertheless, when using the mean-monthly fields of LGM minus present-day anomalies of temperature and precipitation rate to drive a three-dimensional thermomechanical ice-sheet model, it was demonstrated that within realistic bounds of the ice-flow and mass-balance parameters, veryreasonable reconstructions of the Last Glacial Maximum ice sheets could be obtained.

A comparison of different ways of dealing with isostasy: examples from modelling the Antarctic ice sheet during the last glacial cycle

The bedrock isostatic response exerts a strong control on ice-sheet dynamics and is therefore always taken into account in ice-sheet models. This paper reviews the various methods normally used in the ice-sheet modelling community to deal with the bedrock response and compares these with a more sophisticated full-Earth model. Each of these bedrock treatments, five in total, is coupled with a three-dimensional thermomechanical ice-sheet model under the same forcing conditions to simulate the Antarctic ice sheet during the last glacial cycle. The outputs of the simulations are compared on the basis of the time-dependent behaviour for the total ice volume and the mean bedrock elevation during the cycle and of the present rate of uplift over Antarctica. This comparison confirms the necessity of accounting for the elastic bending of the lithosphere in order to yield realistic bedrock patterns. It furthermore demonstrates the deficiencies inherent to the diffusion equation in modelling the complex deformation within the mantle. Nevertheless, when characteristic parameters are varied within their range of uncertainty, differences within one single method are often of the same order as those between the various methods. This overview finally attempts to point out the main advantages and drawbacks of each of these methods and to determine which one is most appropriate depending on the specific modelling requirements.

A comparison of different ways of dealing with isostasy: examples from modelling the Antarctic ice sheet during the last glacial cycle

The bedrock isostatic response exerts a strong control on ice-sheet dynamics and is therefore always taken into account in ice-sheet models. This paper reviews the various methods normally used in the ice-sheet modelling community to deal with the bedrock response and compares these with a more sophisticated full-Earth model. Each of these bedrock treatments, five in total, is coupled with a three-dimensional thermomechanical ice-sheet model under the same forcing conditions to simulate the Antarctic ice sheet during the last glacial cycle. The outputs of the simulations are compared on the basis of the time-dependent behaviour for the total ice volume and the mean bedrock elevation during the cycle and of the present rate of uplift over Antarctica. This comparison confirms the necessity of accounting for the elastic bending of the lithosphere in order to yield realistic bedrock patterns. It furthermore demonstrates the deficiencies inherent to the diffusion equation in modelling the complex deformation within the mantle. Nevertheless, when characteristic parameters are varied within their range of uncertainty, differences within one single method are often of the same order as those between the various methods. This overview finally attempts to point out the main advantages and drawbacks of each of these methods and to determine which one is most appropriate depending on the specific modelling requirements.

Basal temperature conditions of the Greenland ice sheet during the glacial cycles

A high-resolution, three-dimensional thermomechanical ice-sheet model, which includes isostasy, the possibility of ice-sheet expansion on the continental shelf and refined climatic parameterizations, was used to investigate the basal thermal regime of the Greenland ice sheet. The thermodynamic calculations take into account the usual terms of heat flow within the ice, a thermally active bedrock layer and all of the effects associated with changes in ice thickness and flow pattern. Basal temperature conditions are documented with respect to glacial–interracial shifts in climatic boundary conditions, both in steady state as during simulations over the last two glacial cycles using the GRIP δ180 record. It is found that the basal temperature field shows a large sensitivity in steady-state experiments but that, during a glacial cycle, basal temperature variations are strongly damped, in particular in central areas. A comparison has been made with measured data from deep ice cores and the implications are discussed.

Basal temperature conditions of the Greenland ice sheet during the glacial cycles

A high-resolution, three-dimensional thermomechanical ice-sheet model, which includes isostasy, the possibility of ice-sheet expansion on the continental shelf and refined climatic parameterizations, was used to investigate the basal thermal regime of the Greenland ice sheet. The thermodynamic calculations take into account the usual terms of heat flow within the ice, a thermally active bedrock layer and all of the effects associated with changes in ice thickness and flow pattern. Basal temperature conditions are documented with respect to glacial–interracial shifts in climatic boundary conditions, both in steady state as during simulations over the last two glacial cycles using the GRIP δ 180 record. It is found that the basal temperature field shows a large sensitivity in steady-state experiments but that, during a glacial cycle, basal temperature variations are strongly damped, in particular in central areas. A comparison has been made with measured data from deep ice cores and the implications are discussed.