Sea ice feedbacks influence the isotopic signature of Greenland ice sheet elevation changes: last interglacial HadCM3 simulations

Changes in the Greenland ice sheet (GIS) affect global sea level. Greenland stable water isotope (δ¹⁸O) records from ice cores offer information on past changes in the surface of the GIS. Here, we use the isotope-enabled Hadley Centre Coupled Model version 3 (HadCM3) climate model to simulate a set of last interglacial (LIG) idealised GIS surface elevation change scenarios focusing on GIS ice core sites. We investigate how δ¹⁸O depends on the magnitude and sign of GIS elevation change and evaluate how the response is altered by sea ice changes. We find that modifying GIS elevation induces changes in Northern Hemisphere atmospheric circulation, sea ice and precipitation patterns. These climate feedbacks lead to ice-core-averaged isotopic lapse rates of 0.49 ‰ (100 m)⁻¹ for the lowered GIS states and 0.29 ‰ (100 m)⁻¹ for the enlarged GIS states. This is lower than the spatially derived Greenland lapse rates of 0.62–0.72 ‰ (100 m)⁻¹. These results thus suggest non-linearities in the isotope–elevation relationship and have consequences for the interpretation of past elevation and climate changes across Greenland. In particular, our results suggest that winter sea ice changes may significantly influence isotope–elevation gradients: winter sea ice effect can decrease (increase) modelled core-averaged isotopic lapse rate values by about −19 % (and +28 %) for the lowered (enlarged) GIS states, respectively. The largest influence of sea ice on δ¹⁸O changes is found in coastal regions like the Camp Century site.

1 Introduction

Ice core records of stable water isotopes (δ¹⁸O) yield useful information on past climate change over the last several glacial–interglacial cycles (e.g. Sime et al., 2009). Alongside site elevation, ice core δ¹⁸O is affected by variations in site temperature, sea ice, evaporation conditions and transport pathway effects (Sime et al., 2013; Werner et al., 2018). Improving our understanding of the elevation signal captured in Greenland ice core isotopic records means we also need to consider these impacts.

The last interglacial period (LIG – between around 130 and 115 ka) was the last time when the volume of the Greenland ice sheet (GIS) is believed to have been considerably reduced (Robinson et al., 2011; Stone et al., 2013). This period was characterised by warmer-than-present-day conditions in the high latitudes and stronger summer time insolation (Hoffman et al., 2017; Capron et al., 2017). Over Arctic land areas, LIG summer temperatures are estimated to have been around 4–5 ^∘C above present-day values (e.g. CAPE Last Interglacial Project Members, 2006) and the North Greenland Eemian Ice Drilling (NEEM) ice core record suggests surface temperatures were 8 ± 4 ^∘C warmer compared to the last millennium (when accounting for elevation changes in the GIS) (NEEM community members, 2013).

LIG global mean sea level is believed to have risen by between 6 and 9 m compared to present-day levels (Kopp et al., 2009; Dutton et al., 2015) which likely indicates both reduced Antarctic (DeConto and Pollard, 2016; Sutter et al., 2016) and Greenland ice sheets (Tarasov and Peltier, 2003; Lhomme et al., 2005; Helsen et al., 2013; Calov et al., 2015). The contribution of the GIS to this LIG high stand remains unclear: previous studies suggest a possible contribution anywhere between +0.3 and +5.5 m to global mean sea level (Cuffey and Marshall, 2000; Robinson et al., 2011; Quiquet et al., 2013; Stone et al., 2013; Plach et al., 2019). Interestingly though, some total air content measurements have been interpreted as indicative of the elevation over central Greenland remaining nearly unchanged (only few hundred metres lower than today) (Raynaud et al., 1997), and the NEEM air content data have also been interpreted as indicative of a lowering of the surface elevation of only 130 ± 300 m relative to the present (NEEM community members, 2013).

Seven deep ice cores, that likely contain some LIG ice, have been recovered from the GIS with stable water isotope records (δ¹⁸O and δD): NEEM, North Greenland Ice Core Project (NGRIP), Greenland Ice Sheet Project 2 (GISP2), GRIP, DYE3, Camp Century and Renland (Johnsen and Vinther, 2007; NEEM community members, 2013). Isotope–elevation slopes derived from spatial data from central and northwest Greenland suggest that a change in elevation of 100 m may provide a 0.62 ‰ and 0.72 ‰ change in δ¹⁸O, respectively (Dansgaard, 1973; Johnsen et al., 1989; Poage and Chamberlain, 2001). Though the suggested global average isotope lapse rate is 0.3 ‰ (100 m)⁻¹ (Blisniuk and Stern, 2005), it is possible that isotope–elevation relationships vary more widely at high latitudes because of the higher temperature variability (Rowley and Garzione, 2007) or sea ice effects (Holloway et al., 2016a, 2017; Malmierca-Vallet et al., 2018).

While the LIG global average warming is in line with projections for the end of the century (Clark and Huybers, 2009; Hoffman et al., 2017), there is considerable uncertainty on the timing of the sea level high stand during the LIG (Kopp et al., 2009; Düsterhus et al., 2016; Barlow et al., 2018) and the magnitude considerably surpasses near-future (2100–2200) projections (e.g. Fischer et al., 2018; Irvali et al., 2020). Thus, the LIG represents a relevant period, when the implications of changes in ice sheet elevation are highly pertinent for midterm to long-term future projections (DeConto and Pollard, 2016). An improved understanding of the isotopic response to GIS elevation changes may therefore help improve the interpretation of LIG Greenland isotope data and help constrain the GIS response to future sea level and temperature scenarios.

In this study, we first investigate the isolated impact of GIS elevation changes on Greenland δ¹⁸O and the underlying processes (Sect. 3.1). We perform a suite of idealised elevation change simulations with the isotope-enabled Hadley Centre Coupled Model version 3 (HadCM3) climate model to analyse the response of Greenland temperature (Sect. 3.1.2) and precipitation (Sect. 3.1.4), Arctic sea ice (Sect. 3.1.5) and atmospheric circulation (Sect. 3.1.3) to these GIS elevation changes. The second part of this study focuses on testing to which extent variations in the background climate state (Arctic sea ice extent) may influence the isotopic lapse rate values at different Greenland ice core sites (Sect. 3.2). For this analysis, we additionally use a second set of simulations that investigate the joint impact of Arctic sea ice change and GIS changes.

2 Methods

2.1 Model description

We use the isotope-enabled HadCM3 general circulation model (GCM) to simulate the isotopic response to idealised variations in the elevation of the GIS. This GCM has been widely used to examine present, past and future climates (Tindall et al., 2009, 2010; Sime et al., 2013; Holloway et al., 2016b) and consists of a coupled ocean, atmosphere and sea ice model. The atmosphere component has a horizontal grid spacing of 2.5^∘ (latitude) by 3.75^∘ (longitude) and has 19 vertical levels (Gordon et al., 2000). The horizontal grid resolution of the ocean component is 1.25^∘ by 1.25^∘ with 20 vertical levels (Gordon et al., 2000). Tindall et al. (2009) presents the implementation of the water isotope code in HadCM3.

The HadCM3 sea ice output over the Arctic Ocean has been previously validated against observational sea ice data (Meier et al., 2017; Peng et al., 2013) by Malmierca-Vallet et al. (2018). Under PI conditions, HadCM3 simulates rather little summer sea ice over the central Arctic Ocean and too much winter sea ice over the Norwegian, Barents, Labrador and Bering seas. For a full validation of the sea ice model, interested readers are referred to Malmierca-Vallet et al. (2018) and Gordon et al. (2000).

2.2 Experimental setup

2.2.1 GIS elevation change simulations

We run a first ensemble of 16 idealised elevation change HadCM3 simulations with greenhouse-gas and orbital forcing centred at 125 000 years BP (125 ka) (see Table A1). A 125 ka control experiment (hereafter 125Control) is performed including a present-day GIS configuration (IceBridge BedMachine Greenland, version 3; Morlighem et al., 2017a, b). To generate the idealised elevation changes, we scale up and down the GIS height. In particular, we scale elevations relative to the elevation at the NEEM ice core site, following

$$\begin{array}{}\text{(1)}& \mathit{\beta }=\mathrm{\Delta }z/Z_{\mathrm{NEEM}},\end{array}$$

where Δz is the elevation change prescribed, Z_NEEM is the elevation at the NEEM ice core site in the present-day GIS configuration, and β is the scaling percentage. GIS elevations are then decreased/increased by β (ranging between ± 2 % and ± 48 %; Table A1):

$$\begin{array}{}\text{(2)}& Z_{\mathrm{new}}=Z_{\mathrm{ini}}\pm \left(\mathit{\beta }\times Z_{\mathrm{ini}}\right),\end{array}$$

where Z_ini is the two-dimensional array for the present-day GIS and Z_new is a new two-dimensional array with modified elevations. Our simulations use Δz equal to ± 50, 100, 300, 500, 700, 900, 1000 and 1300 m (see Table A1). The 16 idealised simulations are identical except for the GIS elevation, which is decreased/increased by ± 2 % up to ± 48 % (β scaling percentage; Table A1). This decrease/increase is applied evenly over the entire GIS. This simple method is used because it shows a well-controlled and more comprehensible idealised framework for sensitivity studies about the dependence of δ¹⁸O and climate on the magnitude and sign of GIS elevation changes.

To help isolate the impact of elevation changes, the present-day GIS extent is unmodified. Each elevation change simulation is time integrated for a total of 475 years, which ensures appropriately spun-up atmosphere and upper ocean fields. In each case, the final 50 years of the simulations are considered for averaging.

2.2.2 Joint impact of changes in sea ice and GIS morphology

A second ensemble of 32 simulations with different GIS configurations and sea ice retreat scenarios is used to help explore the joint impacts of sea ice change and GIS change on Greenland (see Figs. A1 and A2 and Table A1). This second set of simulations help us explore both changes in the extent (land-ice fraction) (Fig. A2) and elevation of the GIS (Fig. A1).

The methodology outlined in Domingo et al. (2020) is used to generate these simulations. The parameterisation of the set of GIS morphologies and sea ice retreat scenarios is performed by means of a principal component analysis (PCA) approach. Due to the spherical geometry of the Earth, the application of a classical PCA to our data would be inappropriate. We thus apply a particular case of generalised PCA analysis (weighted PCA), described in Jolliffe (2002) and Salter et al. (2019). In summary, (1) GIS morphology modes are calculated from an initial ensemble of 14 LIG GIS reconstructions (Robinson et al., 2011; Born and Nisancioglu, 2012; Helsen et al., 2013; Quiquet et al., 2013; Stone et al., 2013; Calov et al., 2015; Langebroek and Nisancioglu, 2016), and (2) associated relevant sea ice retreat scenarios are calculated from 22 Arctic sea ice change experiments performed by Malmierca-Vallet et al. (2018) (see Fig. A1 and Table A1). In particular, we generate a set of 32 nine-dimensional random vectors. The first eight components of each vector are independently normally distributed and are used to generate a new GIS morphology via linear combination of the PCs (same procedure as in Domingo et al., 2020). The ninth component represents instead heat flux and is uniformly distributed between 0 and 120 W m⁻². We follow the methodology of Holloway et al. (2016a, 2017) and Malmierca-Vallet et al. (2018) on sea ice forcing and implement these heat fluxes to the bottom of the Arctic sea ice (see Table A1). This sea ice forcing is kept constant through the complete annual cycle, and thus the model still calculates the seasonal cycle of sea ice growth and decay. Sea ice varies over time with the coupled model, and both the oceanic and atmospheric components of HadCM3 respond to variations in sea ice (for more details on the methodology, see Malmierca-Vallet et al., 2018).

The resulting 32 LIG GIS morphologies show strong variation in terms of both height and ice extent (Figs. A1 and A2). Some morphologies show a rather small retreat of the GIS, and others show a possible division of the GIS into two domes; some display strong ice loss in the south, while others show substantial ice retreat in the north.

After the design GIS morphologies and associated sea ice forcing are generated, HadCM3 is used to model the isotopic response to these modified GIS morphologies and sea ice retreat scenarios at the ice core sites. All LIG experiments are forced with orbital parameters and greenhouse-gas (GHG) values appropriate for 125 ka and time integrated for a total of 475 years.

3 Isotopic simulation results

We present results from (1) 16 GIS elevation-change-only scenarios and, (2) 32 experiments with combined Arctic sea ice forcing and modified GIS morphology (see Table A1). A two-sided Student's t test is utilised to estimate the statistical significance of changes (von Storch and Zwiers, 2001). Hereafter, lapse rates are defined to be positive, if the analysed atmospheric variable decreases with elevation.

3.1 GIS elevation change scenarios

3.1.1 Mean annual δ¹⁸O changes at ice core sites

At the NEEM deposition site (around 205 ± 20 km upstream of the NEEM drill site due to ice flow; NEEM community members, 2013), the 125Control experiment simulates a precipitation-weighted δ¹⁸O (hereafter δ¹⁸Op) anomaly of 1.9 ‰ compared to PI. The lowered GIS experiments have δ¹⁸Op anomalies which vary between 0.6 ‰ and 6.4 ‰ whilst the increased GIS elevation experiments act to decrease δ¹⁸Op anomalies by as much as −3.9 ‰ in the most extreme scenario (p1300) relative to 125Control (Fig. 1a).

Figure 1 Change in δ¹⁸Op, temperature, precipitation flux and winter (March) sea ice extent as a function of GIS elevation changes (m). Changes are calculated as anomalies compared to 125 ka control experiment. Ice core sites displayed: (a–d) NEEM, (e–h) NGRIP, (i–l) GRIP, (m–p) GISP2, (q–t) DYE3 and (u–x) Camp Century. Results for each of the 16 GIS elevation change scenarios are represented by black dots. Solid blue lines signify best fit lines (y=ax). Also shown are ± 3 SD (blue-shaded envelopes) on the best fit lines. Additionally, a second fit of the data using an exponential function ($y=\pm a{e}^{-bx}+c$) is included (red lines with dashes). Note that in the last column of plots, winter sea ice extent vs. GIS elevation change differs for the various ice core sites because of the different elevation changes at each ice core site compared to the 125 ka control.

The 125Control experiment shows δ¹⁸Op anomalies of 1.5 ‰ at GRIP and GISP2 and 1.6 ‰ at NGRIP compared to PI. Depending on the prescribed reduction of the ice sheet elevation, δ¹⁸Op anomalies compared to 125Control vary between 0.5 ‰ and 6.5 ‰ at NGRIP (Fig. 1e), between 0.3 ‰ and 6.8 ‰ at GISP2 (Fig. 1m) and between 0.2 ‰ and 6.6 ‰ at GRIP (Fig. 1i). In contrast, the increased elevation scenarios show a decline in δ¹⁸Op anomalies of up to −3.6 ‰ at NGRIP and GRIP and −4.0 ‰ at GISP2 in the most extreme scenario (p1300) relative to 125Control (Fig. 1e, i and m).

With respect to the PI simulation, the 125Control simulation shows a δ¹⁸Op rise of 1.0 ‰ at DYE3 and 1.8 ‰ at Camp Century. The decreased elevation scenarios present δ¹⁸Op anomalies ranging from 0.4 ‰ to 5.7 ‰ at DYE3 and from 0.7 ‰ to 4.6 ‰ at Camp Century compared to 125Control (Fig. 1q and u). δ¹⁸Op anomalies decrease up to −3.9 ‰ at DYE3 and −3.3 ‰ at Camp Century in p1300 compared to 125Control (Fig. 1q and u).

δ¹⁸Op anomalies are weaker for increases in GIS elevation than for decreases (Fig. 1). Core-averaged δ¹⁸Op anomalies are −3.7 ‰ compared to 6.1 ‰ for a 1300 m increase/decrease in relative elevation, respectively. This results in a non-linear isotopic lapse rate across the elevation change scenarios (Fig. 1), which will be discussed in the following sections.

For the rest of the results section, for clarity, we focus particularly on two example scenarios which depict medium–high GIS elevation changes (experiments marked in red in Table A1).

3.1.2 Surface air temperatures

The orbital forcing dominates the climate in the 125Control simulation. In Greenland, summer local temperature increases exceed 3.5 ^∘C due to the large increase in summertime insolation (Fig. 2c).

Figure 2 Modelled annual (ANN), winter (DJF) and summer (JJA) surface air temperature anomalies for m900 (d–f) and p900 (g–i) compared to the 125Control simulation. Also shown are temperature anomalies for 125Control compared to the PI simulation (a–c). Only the anomalies statistically significant at the 95 % confidence level are displayed.

The local surface climate over Greenland is noticeably affected by local changes in GIS surface elevation. Decreases in GIS elevation act to increase surface air temperatures (SATs) across Greenland, and vice versa (Figs. 1 and 2). In Greenland, the scenario with decreased elevation (m900) simulates positive SAT anomalies all year round compared with 125Control (Fig. 2d–f). Annual local temperature increases exceed 4.5 ^∘C in m900 relative to 125Control. As expected, the increased elevation scenario (p900) shows negative SAT anomalies throughout the year relative to the 125Control experiment (Fig. 2g–i). In central regions over Greenland, local temperature decreases exceed −4.5 ^∘C during both summer and winter seasons.

A previous study by Merz et al. (2014a) found that the sensitivity of Greenland's LIG climate to GIS topography changes is seasonally diverse. For example, in winter, strong cooling conditions are found over some areas that become ice free and flat, while the remaining ice domes show warmer conditions (Merz et al., 2014a). Areas that become ice free are characterised by weak surface winds and turbulence, barring an efficient sensible heat flux and leading to very strong temperature inversion (Merz et al., 2014a). Our idealised elevation change simulations do not show this temperature inversion mechanism; this is most likely linked to the unchanged land-ice distribution in our experiments.

Averaging across six ice core sites (Camp Century, NEEM, NGRIP, GRIP, GISP2 and DYE3), temperature lapse rates vary slightly from 0.47 ^∘C (100 m)⁻¹ for the lowered GIS states to 0.44 ^∘C (100 m)⁻¹ for the enlarged GIS states (Fig. 1, solid blue lines); however, these changes are not statistically significant.

3.1.3 Atmospheric circulation

To better understand the variations in atmospheric circulation that occur in response to changes in surface elevation, we show changes in the low-level wind pattern (at 850 hPa) and mean sea level pressure (MSLP) field. The 125Control simulation exhibits a widespread decrease in summer (JJA) MSLP compared to PI (Fig. 3c); the warmer SATs and Arctic sea ice loss in 125Control result in a warmer and less stable atmosphere at northern high latitudes during summer. During both winter and summer, the 125Control experiment shows no major differences in wind direction or strength relative to PI (Figs. 4a–d and A3a–b).

Figure 3 ANN, DJF and JJA mean sea level pressure anomalies (Pa) for m900 (d–f) and p900 (g–i) compared to the 125Control simulation. Also shown are sea level pressure anomalies for 125Control compared to the PI simulation (a–c). Only the anomalies statistically significant at the 95 % confidence level are displayed.

Figure 4 Absolute DJF and JJA low-level winds (at 850 hPa) for PI (a–b), 125Control (c–d), m900 (e–f) and p900 (g–h). Shading displays wind speed (m s⁻¹).

Over the Norwegian Sea, there is an increase in winter MSLP (local increases exceed +50 Pa) in m900 relative to 125Control (Fig. 3e). This increase is coincident with a sea ice increase (Fig. 5d) and cooler SATs (Fig. 2e) over the same region compared to 125Control. Around northern Greenland, scenario m900 shows a decrease in summer MSLP relative to 125Control (local decreases exceed −50 Pa; Fig. 3f); this is coincident with an decline in sea ice concentration (Fig. 5a). Scenario p900 shows an increase in annual MSLP over the central Arctic Ocean (Fig. 3g).

Figure 5 Sea ice concentration anomalies (%) for summer (September) and winter (March) for scenarios m900 (a, d) and p900 (c, f) compared to the 125Control simulation. Also shown is absolute sea ice concentration for the 125Control simulation for summer (b) and winter (e). Only the anomalies statistically significant at the 95 % confidence level are displayed.

Over Greenland, the surface winds respond to variations in the GIS surface elevation. Strong anticyclonic flow centred over Greenland is characteristic of the PI and 125Control simulations (Figs. 4 and A3). In the decreased elevation experiment (m900), the Greenland anticyclone becomes smaller (especially during the winter months compared to 125Control; Figs. 4c–f and A3c) and, over northeast Greenland, local wind vectors suggest air mass inflow from the Arctic Ocean contrary to the common outflow observed in the PI and 125Control (Figs. 4 and A3c). In contrast, the scenario with increases in the GIS elevation (p900) displays an enhanced anticyclonic flow particularly during winter (Figs. 4g–h and A3e).

3.1.4 Changes in precipitation pattern

During summer, 125Control shows an enhanced precipitation rate compared to PI mainly across southwestern and central Greenland (Fig. 6c). This is in line with results from other climate models (e.g. Otto-Bliesner et al., 2006; Merz et al., 2014b).

Figure 6 Annual (ANN), winter (DJF) and summer (JJA) precipitation anomalies for m900 (d–f) and p900 (g–i) compared to the 125Control simulation. Also shown are precipitation anomalies for 125Control compared to the PI simulation (a–c). Only the anomalies statistically significant at the 95 % confidence level are displayed.

There is a rise in precipitation rate over much of Greenland throughout the year in m900 compared to 125Control (Fig. 6d to f). This is expected as the lowering of the orography leads to a wider spread of precipitation across Greenland from the east and west which is blocked by the higher and steeper elevation of the present-day GIS. Local increases over southeast Greenland exceed 0.8 mm d⁻¹ in m900 during winter (Fig. 6e). This increase in precipitation is in accordance with a reduction in winter sea ice concentration along the east coast of Greenland relative to 125Control (Fig. 5d).

For the increased elevation scenarios, local changes in precipitation rate relative to 125Control are less widespread and smaller than for the decreased elevation scenarios during both seasons (Fig. 6d–i). Over southeast Greenland, p900 is up to 0.6 mm d⁻¹ drier than the 125Control simulation during winter (Fig. 6h).

Precipitation increases linked to elevation decreases are much larger than the drying linked to elevation increases, implying non-linearities in the climate response to GIS elevation change (Figs. 1 and 6). The core-average (averaging across the six core sites: NEEM, NGRIP, GRIP, GISP2, DYE3 and Camp Century) precipitation lapse rate varies from 0.097 mm yr⁻¹ (100 m)⁻¹ for the decreased GIS elevation states to 0.009 mm yr⁻¹ (100 m)⁻¹ for the increased GIS elevation states (Fig. 1, solid blue lines). Nevertheless, this is considerably influenced by the DYE3 and Camp Century ice core sites (Fig. 1s and w, solid blue lines). The DYE3 ice core site shows a much steeper relationship than the other ice core sites: 0.26 mm yr⁻¹ (100 m)⁻¹ for the decreased elevation scenarios and 0.134 mm yr⁻¹ (100 m)⁻¹ for the increased elevation scenarios (Fig. 1s, solid blue lines). At the Camp Century site, precipitation tends to increase with increasing surface elevation (−0.091 mm yr⁻¹ (100 m)⁻¹ elevation increase) (Fig. 1w, solid blue lines). Removing the Camp Century site increases the core-average precipitation gradient to 0.029 mm yr⁻¹ (100 m)⁻¹ elevation increase. The different behaviour found at the Camp Century site is likely linked to the orographic enhancement of precipitation (Johnson and Hanson, 1995; Frei and Schär, 1998; Petersen et al., 2004; Roe and Baker, 2006) among the enlarged GIS states. Reductions in Camp Century height also result in marginal increases in precipitation rate (0.05 mm yr⁻¹ (100 m)⁻¹) which are probably related to the weakening of the Greenland anticyclone and smaller barrier effect.

3.1.5 Changes in sea ice

For the PI simulation, the September mean sea ice extent is 5.8×10⁶ km². The 125Control simulation shows a reduced September mean of 4.4×10⁶ km² relative to PI; larger seasonal and latitudinal insolation variations (linked to the orbital forcing) lead to Arctic sea ice loss during summer/spring (e.g. Otto-Bliesner et al., 2006).

GIS elevation reductions lead to an increase in winter sea ice extent, whereas increases in the GIS elevation result in winter sea ice retreat (Fig. 1). Thus, the loss/increase of winter sea ice extent acts as a positive/negative feedback on δ¹⁸O. In contrast to δ¹⁸O and SAT, variations in winter sea ice extent are smaller for decreases in GIS elevation compared to increased elevation scenarios. For example, the March sea ice extent is reduced by −4.2 % in p900 and increased by +1.7 % in m900 compared to the 125Control simulation.

The decreased elevation scenario (m900) displays an increase of winter sea ice concentration on the Norwegian Sea and on the southeastern coast of Greenland compared to the 125Control simulation (Fig. 5d). The reduced cyclogenesis off the southeast coast of Greenland (Fig. A5) results in growth of winter sea ice over these regions (Fig. 5d). This is probably associated with a decrease in wind-driven ocean heat transport (e.g. Stone and Lunt, 2013; Davini et al., 2015). The increased elevation scenario (p900) experiences the same forcing but in the opposite direction (Figs. 5f and A5).

We also find some local changes in summer sea ice concentration; while p900 shows decreases of summer sea ice over the Beaufort Sea, it shows increases over the Fram Strait area. Similar patterns are found in m900 but in the opposite direction and of lower magnitude (Fig. 5a, c).

We ascribe these changes in summer sea ice to variations in ocean salinity caused by anomalous downwelling or upwelling, induced by anomalously low or high sea level pressure over the Arctic (Jackson and Vellinga, 2012). In HadCM3, the geostrophic balance of the Beaufort gyre can be altered ageostrophically by wind stresses linked to low-frequency sea level pressure variability (Jackson and Vellinga, 2012). Our increased elevation scenario (p900) shows high sea level pressure anomalies over the Arctic basin (Fig. 3) which lead to downwelling in the centre of the Arctic basin and upwelling along the coasts, respectively (Fig. A4). Since the surface water is fresher and colder than the subsurface water, this results in salinification near the coasts and freshening in the centre of the basin (Fig. A4). The same mechanisms apply to the decreased elevation scenario (m900) but in the opposite direction (Fig. A4).

The increase in wind speed along the Fram Strait in p900 compared to 125Control, and vice versa for m900 (Fig. A5), also affects the advection of sea ice from the Arctic to the Atlantic Ocean (Davini et al., 2015).

3.2 The response of the isotopic lapse rate to changes in the background climate state

Malmierca-Vallet et al. (2018) demonstrate the importance of Arctic sea ice changes as a control on LIG Greenland ice core δ¹⁸O because of its impact on both the regional temperature increase and the moisture source. Thus, we also study 32 simulations that examine the joint impact of modified Arctic sea ice retreat and modified GIS morphology (considering both changes in the extent and elevation of the GIS; see Table A1 and Sect. 2.1).

We use the sea ice retreat simulations of Malmierca-Vallet et al. (2018) to isolate the impacts of δ¹⁸O due to sea ice variation. This allows us to test to which extent Arctic sea ice changes may influence isotopic lapse rate values. Figure A6 shows the change in δ¹⁸O as a function of winter (March) sea ice retreat. We (1) remove the orbital forcing effect by calculating δ¹⁸O anomalies compared to the 125 ka control simulation, and (2) we only analyse scenarios with winter sea ice retreat lower than 55 %, due to the lack of sensitivity of Greenland δ¹⁸O to sea ice losses greater than 50 % (for more details, see Malmierca-Vallet et al., 2018).

Figure 7 δ¹⁸Op anomalies as a function of ice core site elevation change (m) relative to PI (first column) and 125Control (second column). Also shown are sea-ice-corrected δ¹⁸Op anomalies compared to 125Control (third row). Ice core sites displayed: NEEM (first row), NGRIP (second row), GRIP (third row), GISP2 (fourth row), DYE3 (fifth row) and Camp Century (sixth row). Triangles represent results for the 16 elevation change experiments. Dots represent results for the 32 simulations that examine the joined impact of Arctic sea ice retreat and modified GIS shape. Solid lines signify best fit curves (first column, $y=a+bx$; second and third columns, y=ax) and shade envelopes represent ± 3 SD uncertainty on the best fit lines.

To calculate the sea-ice-corrected δ¹⁸O anomalies, we deduct the sea-ice-associated δ¹⁸O effect from the total δ¹⁸O anomalies. Figure 7 shows the resulting sea-ice-corrected δ¹⁸O anomalies (Fig. 7, purple curve fits) as well as total δ¹⁸Op anomalies (compared to both the PI and 125Control simulations) not corrected for sea ice changes (Fig. 7, red and blue curve fits, respectively) for both sets of simulations (Fig. 7, triangles for elevation change simulations and dots for simulations looking at the combined impact of sea ice retreat and GIS changes).

When considering total δ¹⁸Op anomalies (relative to 125Control) not corrected for sea ice changes, a non-linear δ¹⁸O lapse rate is observed over Greenland (Fig. 7, second column). The core-average δ¹⁸O lapse rate varies from 0.29 ‰ (100 m)⁻¹ for the enlarged GIS states to 0.49 ‰ (100 m)⁻¹ for the lowered GIS states (Fig. 7, second column). These results thus strongly suggest a non-linearity in the isotope–elevation relationship, with higher δ¹⁸O–elevation gradients for lowered GIS states, and vice versa.

When further deducting the winter sea ice effect, we find an almost stationary core-average δ¹⁸O lapse rate, slightly varying from 0.38 ‰ (100 m)⁻¹ for the enlarged GIS scenarios to 0.39 ‰ (100 m)⁻¹ for the lowered GIS scenarios (Fig. 7, third column). The sea ice effect increases δ¹⁸O–elevation gradients by 28 % in the enlarged GIS states and decreases δ¹⁸O–elevation gradients by −19 % in the lowered GIS states. Indeed, this suggests that sea ice changes may strongly influence linearity in the isotope–elevation relationship.

4 Discussion

4.1 Response of Arctic sea ice and atmospheric circulation to GIS elevation changes

Our lowered GIS experiments show similar climate behaviour to previous studies where GIS is removed (Toniazzo et al., 2004; Stone and Lunt, 2013; Davini et al., 2015). During summer, warming over the GIS is enhanced due to its lower elevation and Arctic sea ice retreat (Lunt et al., 2004). During winter, warm anomalies extend over the Arctic Ocean. These anomalies are related to a smaller anticyclone over Greenland (Fig. 4h) (Stone and Lunt, 2013; Davini et al., 2015). Support also comes from Merz et al. (2014a), who show, for their perturbed LIG experiments with reduced GIS, a smaller anticyclone as well as decreased wind velocities. The lowering of the elevation also leads to a weakened ice sheet barrier effect, permitting cyclonic systems to get into more central and northern areas of Greenland. These modelled results are in agreement with the findings of Merz et al. (2014a) and Hakuba et al. (2012), who show that a decrease in the height and size of the GIS weakens the barrier effect, permitting more moisture to be advected to the plateau. The reduction in cyclogenesis over the Norwegian Sea and off the southeast coast of Greenland due to lowering of the GIS elevation leads to the growth of further sea ice, especially during the winter months (Fig. 5g and h). Reduced ocean heat transport due to weakened wind-driven currents (Stone and Lunt, 2013) may also contribute to decreased surface temperatures over the Norwegian and Barents seas and on the southeast coast of Greenland (Fig. 2e and h).

When elevation is increased, colder GIS temperatures occur, but a compensating warming occurs around Greenland (Fig. 2). Similar mechanisms causing this surface temperature pattern were discussed in Merz et al. (2014a), who investigated the sensitivity of Greenland's LIG climate to GIS topography changes. In particular, Merz et al. (2014a) found cooling over areas of higher elevation (eastern Greenland) but warming on the periphery of the ice sheet. GIS topography changes influence the Greenland's surface energy balance through changes in surface winds and turbulent heat fluxes (Merz et al., 2014a). In our experiments, the glacial anticyclone over Greenland intensifies as elevation increases, especially during the winter months (Fig. 4i). This leads to enhanced cyclogenesis over the Barents Sea and southeast coast of Greenland which result in winter sea ice retreat over the same regions (Fig. 5f).

4.2 Lapse rates in response to GIS elevation changes

Averaging across the six ice core sites, we find temperature lapse rates of 0.47 ^∘C and 0.44 ^∘C (100 m)⁻¹ for the decreased/increased elevation change scenarios, respectively (Fig. 1). We also find that the wetting related to decreases in GIS elevation is higher than the drying related to increases in elevation.

Our temperature lapse rates compare well with previous estimates which show that the near-surface temperature lapse rate can generally differ from the free-air lapse rate (gradient of moist adiabatic cooling of 0.65 ^∘C (100 m)⁻¹) (Marshall et al., 2007; Fausto et al., 2009; Gardner et al., 2009; Erokhina et al., 2017). Marshall et al. (2007) monitored 25 sites with altitudes between 130 and 2010 m across the Prince of Wales Icefield and observed a mean daily temperature lapse rate of 0.41 ^∘C (100 m)⁻¹. Furthermore, Gardner et al. (2009) showed temperature lapse rates near four glacier surfaces, in the Canadian High Arctic, of 0.49 ^∘C (100 m)⁻¹ (ablation season mean). Gardner et al. (2009) also suggested that lower temperature lapse rates are expected under a warming climate and linked this to the negative relationship found between lower-troposphere temperatures and slope lapse rates.

4.3 The response of the isotopic lapse rate to the background climate state

Isotope–elevation gradients have tended to be calculated from modern surface data (e.g. Dansgaard, 1973): a present-day spatial relationship is presumed to apply to temporal changes. This disregards any impact that variations in the ice sheet elevation may have on the atmospheric circulation, precipitation patterns and eventually the isotopic composition.

Our idealised elevation change simulations with HadCM3 allow a fuller calculation. We find a smaller core-average δ¹⁸O lapse rate for enlarged GIS states (0.29 ‰ (100 m)⁻¹) than for the lowered GIS states (0.49 ‰ (100 m)⁻¹) (Fig. 7). Hence, δ¹⁸O–elevation gradients do not remain constant across the parameter space of elevation changes. This strongly suggests non-linearities in the isotopic response to Greenland elevation change.

We also find that winter sea ice variations can increase/decrease modelled core-averaged isotopic lapse rate values by about +28 % and −19 % for the enlarged/lowered GIS states, respectively (Fig. 7). These results thus suggest that sea ice variations may have a strong influence on δ¹⁸O–elevation gradients, especially at coastal areas such as the Camp Century ice core site (Fig. 9). In particular, at this location, we find that the winter sea ice effect decreases the δ¹⁸O–elevation gradient by −24 % in the lowered GIS states and increases the δ¹⁸O–elevation gradient by as much as 92 % in the enlarged GIS states (Fig. 9e). While the largest influence of sea ice on δ¹⁸O changes is found at the Camp Century site, the DYE3 site shows the smallest (Fig. 9e–f). These results point to elevation changes as a likely driver (together with GHGs and orbital forcing) on LIG δ¹⁸O changes at the DYE3 ice core site. This is in agreement with previous LIG GIS modelling studies which propose a significant LIG lowering around the DYE3 area, even the total loss of ice (Robinson et al., 2011; Helsen et al., 2013).

Interestingly, LIG isotopic lapse rates and the PI spatially derived isotopic lapse (0.37 ‰ (100 m)⁻¹) modelled with HadCM3 are lower than the modern spatially derived gradients of 0.62 and 0.72 ‰ (100 m)⁻¹ in central and northwest Greenland, respectively (Dansgaard, 1973; Johnsen et al., 1989; Vinther et al., 2009). Furthermore, our LIG core-average δ¹⁸O lapse rates are also somewhat lower than the lapse rate of 0.56 ‰ (100 m)⁻¹ modelled (with the isotope-enabled version of the European Centre Hamburg Model version 4) over Greenland for the LIG period by Sjolte et al. (2014). Note that our modelled isotopic lapse rates contemplate the dynamical response of atmospheric circulation to GIS elevation changes and Arctic sea ice variations, whereas previous studies disregard these effects.

These elevation change simulation results thus suggest possible non-linearity in isotope–elevation gradients. It would be useful if this was checked with other models to assess model dependence in the results. The HadCM3 resolution does not permit it to represent the steep GIS margins; this may be behind some of the model–data mismatches (Toniazzo et al., 2004).

4.4 Implications for NEEM δ¹⁸O and elevation reconstructions

Considering the NEEM elevation reconstruction, which indicates NEEM elevation differences of +45 ± 350 m at 126 ka relative to the present day (NEEM community members, 2013), we find a most likely increase in δ¹⁸O values of between +0.8 ‰ and +3.5 ‰ relative to PI (Fig. 8). This relatively falls within the lower end of the uncertainty range of the reconstruction by Domingo et al. (2020): most likely an LIG δ¹⁸O peak of +3.6 ‰ and uncertainty range between +2.7 ‰ and +4 ‰. The relatively small overlap between the δ¹⁸O record and the elevation reconstruction has already been discussed in Domingo et al. (2020) and could possibly reflect uncertainties attached to the air content NEEM elevation reconstruction method. The methodology depends on making corrections to air content measurements related to insolation and temperature in conjunction with secular variations in surface pressure and winds (Martinerie et al., 1994; Krinner et al., 2000; Raynaud et al., 2007; Eicher et al., 2016). In addition, NEEM air content measurements between 127 and 118.3 ka are known to be affected by surface melting (NEEM community members, 2013).

Figure 8 δ¹⁸Op anomalies as a function of the NEEM ice core site elevation change (m) relative to PI. Solid lines signify best fit curves ($y=a+bx$) and shade envelopes represent ± 3 SD on the best fit lines. Triangles represent results for the 16 elevation change experiments. Dots represent results for the 32 simulations that examine the joined impact of Arctic sea ice retreat (sir) and modified GIS shape.

4.5 Relative influence of sea ice and GIS changes on ice core δ¹⁸O

Considering the maximum reduction in NEEM's surface elevation proposed by the NEEM community members (2013) of −305 m at 126 ka, we find that the impact of LIG orbital sea ice changes appears to be the dominant factor determining δ¹⁸O changes (explaining 60 % of the δ¹⁸O anomaly), followed by GIS-driven sea ice changes (Fig. 9a). This is in agreement with previous studies that show the importance of changes in GIS topography and sea ice retreat to explain the LIG warming at the NEEM ice core site (Merz et al., 2014a, 2016; Guarino et al., 2020).

Figure 9 δ¹⁸Op anomalies as a function of the ice core site elevation change (m) relative to PI (red) and 125Control (blue). Also shown are sea-ice-corrected δ¹⁸Op anomalies relative to 125Control (purple). Ice core sites shown: (a) NEEM, (b) NGRIP, (c) GRIP, (d) GISP2, (e) Camp Century, (f) DYE3. Solid lines signify best fit curves and shade envelopes represent ± 3 SD on the best fit lines. Arrows indicate the influence of insolation and sea ice changes on δ¹⁸O changes. In addition, the purple line (sea-ice-corrected δ¹⁸Op anomalies) represents the influence of GIS changes and associated atmospheric circulation changes on δ¹⁸O changes.

To make a comparable analysis at the other ice core sites, of the relative influence of each factor on determining δ¹⁸O changes, we consider the same reduction in surface elevation of −305 m at the other locations (Fig. 9). We find that LIG orbital sea ice changes are the dominant factor determining δ¹⁸O changes at NGRIP, GRIP and GISP2 (accounting for 55 %–58 % of the δ¹⁸O changes) (Fig. 9b–d), while ice sheet changes appear to have the largest impact on δ¹⁸O changes at the DYE3 site (accounting for 48 % of the δ¹⁸O anomaly) (Fig. 9f). The highest sea ice influence is found at Camp Century (explaining 10 % of the δ¹⁸O changes) (Fig. 9e).

Note that the abovementioned relative influence of each parameter on δ¹⁸O changes should be interpreted with caution; these results could change substantially if we were to consider any other possible elevation change scenario. There is no independent gas content information on elevation changes for DYE3 and Camp Century. Moreover, although there are total air content records that were measured on the GRIP (Raynaud et al., 1997) and NGRIP (Eicher et al., 2016) ice cores, the authors show how complex it is to interpret this proxy in terms of elevation changes at the drilling site.

Additional data on elevation changes together with better dated ice, especially at the DYE3 and Camp Century sites, would be particularly valuable to further assess our quantitative elevation change scenarios. In addition, considering sea ice and ice changes in a joint framework following a Gaussian process emulation approach (Domingo et al., 2020) and also accounting for isostatic change may permit a valuable quantitative assessment of how changes in the GIS affected LIG global sea levels.

5 Conclusions

The results of this study are relevant for the interpretation of past climates from Greenland ice core records. Changing GIS elevation in HadCM3 alters the NH atmospheric circulation, sea ice and precipitation patterns over Greenland and further afield. These climate feedbacks result in lower isotope–elevation gradients for enlarged GIS states, and vice versa. Our results thus point to non-linear δ¹⁸O–elevation gradients over Greenland. We further show that isotopic lapse rate values may be significantly influenced by the background climate, in particular, winter sea ice changes.

These model results highlight the importance of the dynamical response of atmospheric circulation to GIS elevation changes when using isotopic measurements to derive past elevation changes: there may be non-linearities in isotope–elevation relationships. Although the underlying mechanisms need further investigation, our finding has important implications for palaeoclimate studies, in which stationary lapse rates are assumed and are normally based on present-day observations. Intermodel comparison studies would be helpful in further developing our understanding of the isotope–elevation gradient over Greenland and how it varies with the background climate state.

Appendix A

Table A1 Full list of scenarios. The simulations highlighted in bold font are the ones primarily discussed in the text. All LIG simulations are performed with greenhouse-gas and orbital forcing centred at 125 ka (eccentricity – 0.04001; obliquity – 23.80^∘; perihelion – 201.3 d of the year; CO₂ – 276 ppmv; CH₄ – 640 ppbv; N₂O – 263 ppbv). The boundary conditions for the PI simulation are the following: eccentricity – 0.0167; obliquity – 23.45^∘; perihelion – 1.7 d of the year; CO₂ – 280 ppmv; CH₄ – 760 ppbv; N₂O – 270 ppbv.

Note (^∗) that to account for the modelled climate variability for a 125 ka control scenario, we use an average of three 125 ka simulations which feature very minor GIS elevation changes between them (maximum scaling percentage of ± 1.8 in Eq. 1) as our 125 ka control.

Figure A1 Surface elevation (metres) of LIG GIS morphologies used in the 32 simulations that examine the joint impact of modified Arctic sea ice retreat and modified GIS morphology. See Table A1 for more details.

Figure A2 Land-ice mask of LIG GIS morphologies used in the 32 simulations that examine the joint impact of modified Arctic sea ice retreat and modified GIS morphology. See Table A1 for more details.

Figure A3 Winter (DJF) and summer (JJA) low-level wind (at 850 hPa) anomalies for m900 (c–d) and p900 (e–f) compared to the 125Control simulation. Also shown are low-level wind anomalies for 125Control compared to the PI simulation (a–b). Only the anomalies statistically significant at the 95 % confidence level are displayed. Shading displays wind speed (m s⁻¹).

Figure A4 Annual salinity (PSU; averaged over the top 535 m) anomalies for m900 (a) and p900 (b) compared to the 125Control simulation. Only the anomalies statistically significant at the 95 % confidence level are displayed. Shading displays wind speed (m s⁻¹).

Figure A5 Change in 10 m winter (DJF) wind speed (in m s⁻¹) for (a) m900 and (b) p900 compared to 125Control. Only the anomalies statistically significant at the 95 % confidence level are displayed.

Figure A6 Simulated δ¹⁸Op anomalies as a function of winter (March) sea ice retreat. The retreat of sea ice is calculated as the percentage change in winter sea ice extent compared to the 125 ka control simulation. Ice core sites shown: (a) NEEM, (b) NGRIP, (c) GRIP, (d) GISP2, (e) Camp Century and (f) DYE3. Triangles represent results from the sea ice sensitivity experiments performed by Malmierca-Vallet et al. (2018). Solid lines signify best fit lines (fit $=b\times \left(\log \right(x)-a)$). Also shown are ± 1 SD (lines with dashes) and ± 3 SD uncertainties (shade envelopes) on the best fit curve. We assume a straight line regression (fit $=a+bx$) for winter sea ice losses between 0 % and 9 % and increases up to −4.5 %. This figure is an adaptation of Fig. 10 shown in Malmierca-Vallet et al. (2018).

Appendix B: The method of isolating the impacts of δ¹⁸O due to sea ice changes

To calculate the sea-ice-corrected δ¹⁸O anomalies, we deduct the sea-ice-associated δ¹⁸O effect from the total δ¹⁸O anomalies (see Sect. 3.2). In particular, we use the sea ice retreat simulations of Malmierca-Vallet et al. (2018) to isolate the impacts of δ¹⁸O due to sea ice variation (Fig. A6). We acknowledge that this approach has its limitations as the interaction among GIS shape and sea ice factors could lead to smaller/larger effects than predicted from the sum of single parameter effects. In order to test the robustness of our sea-ice-correction method, we run eight additional LIG simulations with the purpose of separating the effect of sea ice changes versus GIS shape changes. From the ensemble of 32 simulations, which explore the joint impacts of sea ice change and GIS change over Greenland, we select four simulations (GIS1-SIE-11.49, GIS2-SIE-11.52, GIS13-SIE-14.98 and GIS31-SIE-19) and rerun them with (1) only the sea ice forcing implemented and (2) only the modified GIS shape implemented (see Table A1).

Table B1 Simulated δ¹⁸Op anomalies compared to 125 ka control on six Greenland deep ice cores: NEEM, NGRIP, GRIP, DYE3, GISP2 and Camp Century. For each pair of simulations, δ¹⁸Op anomalies due to (1) GIS shape changes, (2) sea ice changes and (3) sum of single parameter effects are shown (numbers in bold font).

Table B2 Simulated δ¹⁸Op anomalies compared to 125 ka control on six Greenland deep ice cores: NEEM, NGRIP, GRIP, DYE3, GISP2 and Camp Century. δ¹⁸Op anomalies due to the joint impact of GIS shape changes and sea ice changes are shown.

We find that the four simulations that explore the joint impact of GIS shape changes and sea ice changes result in smaller/larger δ¹⁸Op anomalies (compared to the 125 ka control) than those predicted from the sum of single parameter effects (Tables B1 and B2). Nevertheless, differences are not higher than around ± 1 ‰. This is within the model uncertainty of annual mean δ¹⁸Op (Malmierca-Vallet et al., 2018).

Appendix C: Model evaluation

The performance of the atmospheric component of HadCM3 (HadAM3; Pope et al., 2000) over Greenland has been validated against the European Centre for Medium-Range Weather Forecasts (ECMWF) reanalysis (Murphy et al., 2002). They show reasonable agreement with temperature, precipitation and wind observations except for a small cold winter bias (associated with excessive longwave cooling), a warm summer bias (excessive shortwave heating at the surface) and a wet bias (related to inefficient orographic blocking) (Murphy et al., 2002). The performance of the HadCM3 coupled mode suffers from similar errors. In particular, the model's performance over the Greenland region and when coupled to the land surface scheme MOSES 2.1 has been published in Stone and Lunt (2013) and Valdes et al. (2017). There is generally good agreement with observed temperatures (derived from Hanna et al., 2005) and observational annual precipitation (derived from Uppala et al., 2005), except for a summer warm bias (1.9 ^∘C) and an annual wet bias (1.4 mm d⁻¹) in southeast Greenland for both predicted and prescribed vegetation control runs (Stone and Lunt, 2013).

A validation of the isotope output has also been carried out for the atmosphere only (HadAM3; Sime et al., 2013) as well as for the coupled ocean–atmosphere model (Tindall et al., 2009, 2010; Xinping et al., 2012). HadCM3 is able to reproduce the large-scale features of δ¹⁸O in precipitation, covering altitude, latitude, amount and continental effects (Tindall et al., 2009). Moreover, Malmierca-Vallet et al. (2018) provides an evaluation of two control (PI and present-day experiments) HadCM3 isotope simulations over Greenland; HadCM3 shows similar heavy δ¹⁸O biases over Greenland to other models (e.g. Sime et al., 2013; Sjolte et al., 2014).

Similar biases are expected to affect the PI and LIG experiments. To minimise the effect of the model bias over Greenland, and hence any influence on the study results, we use the standard approach and report modelled values as anomalies (LIG minus PI).

In addition, we note that the coarse spatial resolution of HadCM3 complicates the process of reliably modelling δ¹⁸O changes at the coastal margins. Hence, the small Renland ice cap (where LIG ice has been retrieved) is not included in this study, as it is not well captured within the HadCM3 resolution.