Greenland climate simulations show high Eemian surface melt which could explain reduced total air content in ice cores
This study presents simulations of Greenland surface melt for the Eemian interglacial period (∼130 000 to 115 000 years ago) derived from regional climate simulations with a coupled surface energy balance model. Surface melt is of high relevance due to its potential effect on ice core observations, e.g., lowering the preserved total air content (TAC) used to infer past surface elevation. An investigation of surface melt is particularly interesting for warm periods with high surface melt, such as the Eemian interglacial period. Furthermore, Eemian ice is the deepest and most compressed ice preserved on Greenland, resulting in our inability to identify melt layers visually. Therefore, simulating Eemian melt rates and associated melt layers is beneficial to improve the reconstruction of past surface elevation. Estimated TAC, based on simulated melt during the Eemian, could explain the lower TAC observations. The simulations show Eemian surface melt at all deep Greenland ice core locations and an average of up to ∼30 melt days per year at Dye-3, corresponding to more than 600 mm water equivalent (w.e.) of annual melt. For higher ice sheet locations, between 60 and 150 on average are simulated. At the summit of Greenland, this yields a refreezing ratio of more than 25 % of the annual accumulation. As a consequence, high melt rates during warm periods should be considered when interpreting Greenland TAC fluctuations as surface elevation changes. In addition to estimating the influence of melt on past TAC in ice cores, the simulated surface melt could potentially be used to identify coring locations where Greenland ice is best preserved.
The Eemian interglacial period (∼130 000 to 115 000 years ago; hereafter ∼130 to 115 ka) was the last period with a warmer-than-present summer climate on Greenland (CAPE Last Interglacial Project Members, 2006; Otto-Bliesner et al., 2013; Capron et al., 2014). Favorable orbital parameters (higher obliquity and eccentricity compared to today) during the early Eemian period caused a positive northern summer insolation anomaly (and negative winter anomaly) at high latitudes, which led to a stronger seasonality (Yin and Berger, 2010). This stronger seasonality with relatively warm summer seasons is favorable for high melt rates across the Greenland ice sheet.
Unfortunately, the presence of surface melt can influence our ability to interpret ice core records. Measurements of CH4, N2O, and total air content (TAC) can be affected if melt layers are present. Other ice core measurements such as δ18O, δD, and deuterium excess appear to be only marginally affected (NEEM community members, 2013). However, refrozen melt has the potential to form impermeable ice layers (melt layers henceforth) that alter the diffusion of ice core signals.
The observed TAC of ice core records is the only direct proxy for past surface elevation of the interior of an ice sheet; i.e., the TAC is governed by the density of air which mainly decreases with elevation. However, TAC is also affected by low-frequency insolation variations (changing orbital parameters) at both Antarctic and Greenlandic sites (Raynaud et al., 2007; Eicher et al., 2016). Furthermore, Eicher et al. (2016) find a TAC response on millennial timescales (during Dansgaard–Oeschger events), which is hypothesized to be related to rapid changes in accumulation. While TAC can be estimated for each individual ice core without the need for other reference ice cores, another indirect method which has been applied to infer Holocene thinning of the Greenland ice sheet (Vinther et al., 2009) requires several ice cores. Vinther et al. (2009) compare the changes of δ18O at coastal ice caps (stable surface elevation due to confined topography) with Greenland deep ice cores and infer elevation changes. Unfortunately, Eemian ice core records are sparse, and therefore TAC is the only direct method available to estimate surface elevation changes this far back in time. Since the assumed surface elevation also influences the actual Eemian temperature reconstructions and its uncertainty range, an accurate TAC record is of high importance. The following example illustrates this importance: the North Greenland Eemian Ice Drilling project (NEEM)-derived surface temperature anomaly (NEEM community members, 2013) at 126 ka is 7.5±1.8 ∘C (relative to the last 1000 years) without accounting for elevation changes; including the elevation change based on TAC measurements, the temperature estimate becomes 8±4 ∘C. This means that more than half of the uncertainty of this temperature estimate is related to the uncertainty of past surface elevation.
Despite the importance that melt can have for the interpretation of TAC and other variables of ice core records, the number of studies analyzing the frequency of melt layers in Greenland ice cores is limited (Alley and Koci, 1988; Alley and Anandakrishnan, 1995).
This study investigates regional climate simulations and observations at seven deep Greenland ice core sites – Camp Century, Dye-3, East Greenland Ice-Core Project (EGRIP), Greenland Ice Core Project (GRIP), Greenland Ice Sheet Project 2 (GISP2), NEEM, and North Greenland Ice Core Project (NGRIP). Additionally, an ice cap in the vicinity of the ice sheet is examined – the Agassiz ice cap, located in the northern Canadian Arctic. TAC is derived from regional climate and melt simulations at these locations of interest (Sect. 2). Furthermore, the simulated local temperature and melt are evaluated, and the impact on TAC is estimated and compared with ice core observations (Sects. 3 and 4). The results indicate that Greenland ice core records from warm periods, such as the Eemian interglacial period, might be more affected by surface melt than previously considered (Sect. 5).
2.1 Climate and surface mass balance simulations
This study uses climate and surface mass balance (SMB) based on two Eemian time slice simulations with a fast version of the Norwegian Earth System Model (NorESM1-F; Guo et al., 2018) representing (constant) 125 and 115 ka conditions and one pre-industrial (PI; constant 1850 forcing) control simulation. These global simulations are dynamically downscaled over Greenland with the regional climate model Modèle Atmosphérique Régional (MAR v3.7; 25 km×25 km), which was extensively validated over Greenland under present-day climate conditions (Fettweis, 2007; Fettweis et al., 2013a, 2017).
MAR employs a land surface model (SISVAT; Soil Ice Snow Vegetation Atmosphere Transfer) with a detailed snow energy balance (Gallée and Duynkerke, 1997) fully coupled to the model atmosphere. MAR's atmosphere uses the solar radiation scheme of Morcrette et al. (2008) and accounts for the atmospheric hydrological cycle (including cloud microphysics) based on Kessler (1969) and Lin et al. (1983). The snow–ice component of MAR is derived from the Crocus snowpack model (Brun et al., 1992), simulating mass and energy fluxes between snow layers and reproducing snow grain properties as well as their effect on surface albedo. The MAR model has 24 atmospheric layers (up to 16 km above ground) and SISVAT 30 snowpack layers.
The NorESM-F experiments are spun up for 1000 years with constant 1850 forcing (greenhouse gas (GHG) concentrations and orbital parameters) to a quasi-equilibrium state. The PI simulation is run for another 1000 years with constant forcing. The two Eemian time slice simulations are branched off from the initial 1000-year spin-up and run for another 1000 years each with constant 125 and 115 ka forcing, respectively (changed GHG concentrations and orbital parameters compared to PI). For the MAR experiments, NorESM is run for another 30 years for each of the three experiments and the output is saved on a 6-hourly basis. These 30 years are used as boundary forcing for MAR. After disregarding the first 4 years as spin-up, the final 26 years are used for the analysis (thin lines in Figs. 2–6, A1–A3). All climate simulations use a fixed, modern ice sheet geometry, in the absence of a reliable Eemian ice sheet estimate and high computational costs of a coupling with an ice flow model (e.g., Le clec'h et al., 2019).
The MAR SMB is analyzed in a study investigating the influence of climate model resolution and SMB model selection on Eemian SMB simulations (Plach et al., 2018), which amongst other things shows the high importance of solar insolation in Eemian simulations. Additionally, while providing the most complete representation of physical surface processes in the pool of investigated models, MAR shows a less negative SMB than an intermediate complexity model during the warmest Eemian simulations (mainly due to a higher ratio of refreezing).
Furthermore, the discussed SMB is used in a study investigating the Eemian Greenland ice sheet volume with a higher-order ice sheet model (Plach et al., 2019). Plach et al. (2019) showed that different external SMB forcings show a larger influence on the Eemian ice volume minimum than sensitivity experiments performed with internal ice dynamics (like changed basal friction). The ice sheet simulations with the MAR SMB show a moderately smaller Eemian ice sheet with the difference equivalent to ∼0.5 m of sea level rise (with respect to the modern ice sheet).
In this study, the MAR SMB simulations are analyzed at seven deep Greenland ice core locations – Camp Century, Dye-3, EGRIP, GRIP, GISP2, NEEM, NGRIP – and an adjacent ice cap – the Agassiz ice cap (Fig. 1). Due to model topography misrepresentation at the ice sheet margins, i.e., the model topography is lower than in reality at the Agassiz ice cap location (model resolution 25 km), a substitute location (Agassiz_sub) in the vicinity of the ice cap, with a model elevation similar to the observed elevations, is chosen (Table 1).
Agassiz_sub refers to a substitute location used due to model topography misrepresentation. For details, see Sect. 2.
2.2 Observed surface melt
The PI climate and SMB simulations are compared to present-day satellite and temperature observations at the locations of interest. The two observational melt day datasets are both derived from satellite-borne passive microwave radiometers – Scanning Multichannel Microwave Radiometer (SMMR), the Special Sensor Microwave/Imager (SSM/I), and the Special Sensor Microwave Imager/Sounder (SSMIS). The first dataset, MEaSUREs (Greenland Surface Melt Daily 25 km EASE-Grid 2.0, version 1), covers the years 1979 to 2012 and is available for the entire Northern Hemisphere. The melt onset is identified by comparing 37 GHz, horizontally polarized (37 GHz H-Pol) brightness temperatures with dynamic thresholds associated with a melting snowpack (Mote, 2014). Unfortunately, the Agassiz ice cap is not covered by this dataset. The second dataset, T19Hmelt, covers the whole MAR grid at 25 km from May to September for most years between 1979 and 2010. It uses data collected at K-band horizontal polarization (T19H) with a constant brightness temperature threshold of 227.5 K (Fettweis et al., 2011). Both satellite datasets are discussed to show their different sensitivities and to illustrate the uncertainty of these satellite-based melt observations.
The seasonal temperature observations at weather stations and 10 m borehole temperatures (representing annual mean temperatures from 1890 to 2014) are taken from a collection of shallow ice core records and weather station data (Faber, 2016). Finally, the bore hole temperatures from the Agassiz ice cap are taken from Vinther et al. (2008).
2.3 Observed total air content
Firstly, the Dye-3 TAC for the ice core depth range of ∼240 to 1920 m was extracted from Herron and Langway (1987, Fig. 4 therein). Since Souchez et al. (1998) indicate that ice from warmer periods (higher δO18 values), likely Eemian, is located below 2000 m at Dye-3, the presented Dye-3 TAC record does not represent Eemian conditions. Secondly, the GRIP TAC dataset (Raynaud, 1999) covers depths from ∼ 120 to 2300 and ∼2780 to 2909 m, while an age mode is only provided for the upper part (oldest ice at 41 ka). For the deeper sections of the core, a published unfolding of the GRIP core (Landais et al., 2003, age bands in Fig. 3 therein) is used to assign an age to the observations. Thirdly, the GISP2 TAC data were extracted from a Supplement table of Yau et al. (2016) and cover the period from 127.6 to 115.4 ka. Fourthly, the NEEM TAC observations (NEEM community members, 2013) cover the deepest section of the NEEM ice core from ∼2200 to 2500 m depth (corresponding to an age of ∼75 to 128 ka; not continuous) and an example for Holocene conditions from depths between ∼100 and 1400 m (no age provided). Finally, the NGRIP TAC record (Eicher et al., 2016) includes the entire core from ∼130 to 3080 m; however, the sampling resolution varies. An age model is provided for the entire dataset with a maximum age of ∼120 ka. Note that only the Eemian sections for GRIP, GISP2, NEEM, and NGRIP are shown in Fig. 7.
2.4 Calculation of the model-derived total air content
The model-derived TAC is calculated with the annual mean surface pressure and the annual mean near-surface temperature from the MAR regional climate simulations at every location of interest (Martinerie et al., 1992; Raynaud et al., 1997):
where Vc is the pore volume at close-off in cm3 g−1 of ice, Pc the mean atmospheric pressure at the elevation of the close-off depth interval in mbar, Tc the firn temperature prevailing at the same depth interval in K, P0 the standard pressure (1013 mbar), and T0 the standard temperature (273 K). Vc is calculated as a function of Tc following an empirical relation (Martinerie et al., 1994; Raynaud et al., 1997):
This theoretical TAC is then reduced (TACred) depending on the percentage of refreezing of the annual accumulation (RZper):
where TACrefrozen is calculated using Henry's solubility law (Sander, 2015) for N2 and O2 (neglecting other atmospheric gases) to account for air that is dissolved in the meltwater before refreezing:
with and being the aqueous-phase concentration of N2 and O2, respectively:
where and are the atmospheric concentration ratios (0.79 and 0.21), and and are Henry's solubility constants ( and ) for N2 and O2, respectively. Henry's law assumes that the meltwater is in equilibrium with the ambient air at a temperature of 273 K and at the local atmospheric pressure (Eqs. 5 and 6). No air is occluded in the form of bubbles in the freezing process.
The simulated PI annual mean (near-surface) temperatures (1850 climate forcing) at the eight locations of interest (Fig. 2; black columns; short bold lines – ensemble means; short thin lines – individual model years) generally fit well with observations from weather stations (Fig. 2; long bold lines in black; standard deviation in gray shading). However, the annual means inferred from 10 m borehole temperatures (Fig. 2; long bold lines in gray; average of the years 1980 to 2014) are consistently colder than the simulated PI means. The lower borehole temperatures represent snow temperatures which are typically cooler than the ambient air temperatures. Only at the Agassiz site, the borehole temperatures are higher. This exception is likely related to the usage of a substitute location (see Sect. 2).
The annual mean temperatures at most locations only vary by 0.5 ∘C between the time slice simulations, i.e., no large difference between PI (Fig. 2; black) and warmest Eemian simulations (Fig. 2; orange). This is to be excepted since the annually integrated solar irradiance is similar in all time slices.
However, the varying Eemian seasonality (Yin and Berger, 2010) results in consistently ∼3–4 ∘C (with respect to PI; black) warmer summer (JJA; June–July–August) temperatures at all locations for mid-Eemian conditions (125 ka: orange) and cooler temperatures for late Eemian conditions (115 ka: blue). The simulated PI summer temperatures (Fig. 3; black columns; short bold lines – ensemble means; short thin lines – individual model years) show good agreement with observations from weather stations (Fig. 3, long bold lines in black).
The precipitation-weighted temperatures (Fig. A1) show a similar pattern to the JJA temperatures (Fig. 3). This is understandable since most precipitation in Greenland falls around the summer months and these temperatures are calculated by multiplying daily temperatures with daily precipitation, summing up the results over the year and then dividing by the sum of the annual precipitation; i.e., precipitation is used as a weight, instead of time, in annual mean temperatures. Precipitation-weighted temperatures are arguably closer to what is recorded in an ice core (temperature at the time of deposition) and these temperatures show a less pronounced warming for mid-Eemian conditions (125 ka: orange), i.e., maximum 3 ∘C warmer compared to PI (black).
3.2 Number of melt days
Passive microwave satellite data show a strong difference in observed melt days per year (presence of surface water) (Fig. 4; first three columns from the left; brown and green) between central ice core locations (GRIP, GISP2, NGRIP, NEEM, EGRIP), where surface melt is sparse, and locations closer to the margins (Camp Century, Dye-3) and ice caps (Agassiz), where melt is much more frequent. Central locations show between 0 and ∼1 melt days per year in the last ∼30 years for which satellite data are available. The exact values vary depending on the location, satellite dataset, and whether the extreme melt event of 2012 is included.
The simulated PI melt day frequency (Fig. 4, black columns) shows good agreement with the observations (Fig. 4; brown and green columns), i.e., low melt frequencies at the central locations and higher melt frequencies at locations at the margins. However, the simulated PI melt frequencies are generally lower than present-day observations (especially at the Agassiz location), with the exception of Dye-3, which shows a higher simulated melt frequency.
3.3 Melt and refreezing
The 125 ka simulations (Fig. 4; orange columns) show a significantly higher melt frequency at all locations (more than 30 melt days per year at Dye-3) compared to the PI simulations (Fig. 4; black columns) and observations (Fig. 4; brown/green columns). The SMB simulations show surface melt at all ice core locations during the warm mid-Eemian with an annual meltwater production (Fig. A2) for warmer locations of ∼300 (Camp Century) and ∼ 600 (Dye-3). However, even modern dry, high-altitude locations show an annual surface melt of ∼ 60 (GRIP, GISP2), 80 (NGRIP), and up to 120 (EGRIP). NEEM shows ∼ 150 for the warmest Eemian simulations.
The mean simulated amount of refreezing exceeds 40 % of the annual accumulation at most ice core locations under warm mid-Eemian conditions (Fig. 5; thick orange lines). Even at the highest locations, GRIP and GISP2 at ∼3200 m elevation, refreezing surpasses 25 % of the annual accumulation under 125 ka conditions. The largest amount of refreezing is simulated at Agassiz_sub, EGRIP, and Dye-3, where refreezing percentages reach 80 % to 90 %.
3.4 Total air content
Theoretical TAC derived from simulated surface pressure and annual mean temperature (Raynaud et al., 1997) and reduced according to the amount of simulated refreezing (Fig. 6 and Sect. 2) shows significantly lower values for the 125 ka simulations. Most of the higher ice core locations (GRIP, GISP2, NGRIP, NEEM, EGRIP, and Camp Century) show simulated TAC values between 45 and 70 mL kg−1 on average, whereas the respective PI values are between 90 and 100 mL kg−1. At Dye-3, the simulated TAC is about 25 mL kg−1 on average for the warm 125 ka Eemian simulations compared to 75 mL kg−1 during PI. Observed Holocene TAC from ice core records (Fig. 6; horizontal gray shading) fits well with the PI simulations, while observed Eemian TAC (Fig. 6; horizontal orange shading) is not as low as the simulated values.
The Eemian ranges in Fig. 6 are calculated as the average (± 2 standard deviations) of the lowest 10 % of observed Eemian TAC (Fig. 7; used observations are indicated in orange) for NEEM and NGRIP. Due to the low number of Eemian observations at GRIP and GISP2, a different threshold of 20 % is used for this core. For the calculation of the late Holocene ranges in Fig. 6, observations younger than 1000, 2000, and 4000 years are used for GRIP, Dye-3, and NGRIP, respectively. The late Holocene range for NEEM is calculated from the entire Holocene example provided in the NEEM community members (2013) data (nine data points; no age provided).
Finally, TAC observations from the deeper ice core sections (i.e., possibly Eemian; Fig. 7; NEEM, GRIP, GISP2, NGRIP; circles; inverted y axes) are compared with mean simulated TAC for 115 ka (Fig. 7; blue line) and 125 ka conditions (Fig. 7; orange line). For Dye-3, the entire TAC record is shown due to the lack of Eemian observations. However, the ice at the bottom of Dye-3 has been shown to contain pre-Eemian ice (Willerslev et al., 2007). Note that NEEM and GRIP are shown against age based on a more robust chronology involving “unfolding the ice” (NEEM community members, 2013; Landais et al., 2003), while NGRIP and Dye-3 are shown against core depth.
The 115 ka simulations generally fit well with the late Eemian (NEEM, GRIP, GISP2, NGRIP) and Holocene (Dye-3) observations, while the 125 ka simulations are lower than the observations. For NEEM, the lowest TAC observations are within the ice core section influenced by melt (gray shading in Fig. 7; NEEM community members, 2013).
The enhanced Eemian seasonality (Yin and Berger, 2010) and warmer Eemian summers (CAPE Last Interglacial Project Members, 2006; Otto-Bliesner et al., 2013; Capron et al., 2014) are indicators of elevated melt during this period. The recent extreme melt event in Greenland in 2012 and a similar event in 1889 (Nghiem et al., 2012) demonstrate that surface melt on the entire Greenland ice sheet, even at the summit of Greenland, is possible under recent climate conditions. Even though these extreme Greenland-wide melt events were caused by a rare large-scale atmospheric pattern (Neff et al., 2014) and were further enhanced by an externally caused albedo lowering (ash deposition from forest fires; Keegan et al., 2014), it is likely that such events are more frequent in a warmer climate such as the Eemian interglacial period.
The simulations discussed in this study (regional climate plus a full surface energy balance) indicate surface melt and refreezing (Figs. 4 and 5) at all deep Greenland ice core locations. Even central Greenland locations close to summit (GRIP, GISP2) show a melt of ∼60 mm yr−1 (Fig. A2). Due to this high surface melt, TAC measurements derived from these simulations are between ∼25 % (GRIP, GISP2) and ∼80 % (Dye-3, EGRIP) lower than modern (PI) values (Fig. 6). Even though the presented climate simulations show such extensive melt, there are several reasons why these simulations can be interpreted as conservative estimates: (1) the simulated PI melt frequency is mostly lower than satellite observations (Fig. 4; black vs. brown/green columns). However, the observation of higher melt frequencies can likely also be related to the effects of recent global warming which are not represented in the PI climate simulations. (2) Processes like ash deposition which were partly responsible for the extreme Greenland melt events of 2012 and 1889 (Keegan et al., 2014) are not simulated. (3) The climate simulations use a fixed, modern ice sheet geometry and including the neglected lowering and retreat of the Eemian ice sheet would likely increase the simulate warming in many regions.
Many studies suggest a substantial Eemian ice volume reduction (e.g., Van de Berg et al., 2011) particularly in the marginal regions – an overview of previous Eemian studies can be found in Plach et al. (2018). The use of a fixed ice sheet undoubtedly adds additional uncertainties to the presented melt simulations – e.g., neglecting modifications of local wind patterns and surface albedo as regions become deglaciated impacting local near-surface temperature (Merz et al., 2014a), local orographic precipitation following the slopes of the ice sheet (Merz et al., 2014b), or increased katabatic winds caused by steeper ice sheet slopes (Gallée and Pettré, 1998; Le clec'h et al., 2019). However, these uncertainties are much stronger in marginal than in high-altitude regions where the ice elevation changes were more limited. After all, a future, more exhaustive evaluation of Eemian melt at the ice cores sites should investigate different possible ice sheet geometries.
Furthermore, the absence of a simulated annual warming, and proxy data showing Eemian peak temperatures as high as ∘C (NEEM community members, 2013, without altitude corrections) and ∘C (Landais et al., 2016) for NEEM (the North Greenland Eemian Ice Drilling project in northwest Greenland) and ∘C (Landais et al., 2016, lower bound as the record only starts after the peak Eemian warming) for NGRIP (North Greenland Ice Core Project) indicate that the climate simulations might include a cold bias. The simulated JJA temperatures (Fig. 3) and the simulated precipitation-weighted temperatures (Fig. A1) show a peak warming of only ∼3–4 and ∼3 ∘C, respectively. However, the fact that NEEM community members (2013) infer an elevation (at the deposition site) of several hundred meters higher than at NEEM today complicates the interpretation of how well the simulated temperatures fit the proxy-derived observations.
Focusing again on the comparison of melt observations and simulations (Fig. 4), a strong underestimation of melt at the Agassiz site in the PI simulations becomes apparent. This strong underestimation is likely related to the use of a substitute location (geographically shifted, with similar model and observed elevation) necessary due to low model topography at the original core site causing unrealistically high melt simulations. Furthermore, the Agassiz site is only covered by the satellite dataset which appears to be less sensitive to melt (T19Hmelt with less melt than MEaSUREs at all sites) and although Eemian ice is absent at the Agassiz site, the simulated Eemian refreezing percentage (Fig. 5) of approximately 80 % is consistent with the Agassiz melt record, which indicates a complete melt of the annual accumulation during the Holocene optimum ∼10–11 ka (Fisher et al., 2012; Lecavalier et al., 2017).
Another important aspect for the melt interpretation is the formation of melt layers and the amount of meltwater needed to form a (visible) melt layer. While the presented TAC calculations assume Henry's solubility law (Sander, 2015) for the air content of the melt layer, the formation of a melt layer in an ice core is a complicated process, e.g., depending on prevailing snow properties. A higher number of melt layers is not just the result of uniformly higher summer temperatures but a combination of an increased contrast between the pre-melt snow pack temperatures (strongly influenced by winter temperature) and the summer melt rate (a function of summer temperature) (Pfeffer and Humphrey, 1998). Therefore, the enhanced Eemian seasonality might have been favorable for the formation of melt layers.
The simulated 125 ka TAC values are consistently lower than the observations (Figs. 6 and 7). However, at NEEM – the ice core with the most complete Eemian record (likely including peak warming) – the simulated 125 ka TAC seems to be most similar to the lowest observations, indicating that the high amount of simulated melt could explain these observations. The variability of the observed NEEM TAC in the suggested melt zone between 127 and 118.3 ka (gray shading; NEEM community members, 2013) is large, likely due to the varying influence of melt layers.
The Eemian TAC measurements at GRIP, GISP2, and NGRIP also show reduced values (not as low as at NEEM), which can be interpreted in a similar way as at NEEM – GRIP, GISP2, and NGRIP might have been influenced by Eemian melt as well. The simulated 125 ka TAC measurements for all three locations are strongly reduced (relative to PI levels) but do not reach levels as low as at NEEM. However, these reduced TAC levels could indicate significant surface melt.
Overall, the lack of a better agreement between observed and simulated Eemian TAC (i.e., few TAC observations as low as the simulations) could be related to the sparse number of Eemian peak warming observations (most ice core records only start after the peak warming; particularly at GRIP, GISP2, NGRIP, and Dye-3). However, another possible explanation could be a shift of the precipitation rates in central Greenland towards much higher values during the Eemian interglacial period. Unfortunately, accumulation rates are unconstrained for the Eemian sections of Greenland ice cores.
Furthermore, another uncertainty for the interpretation of the simulations is the effect of the higher Eemian summer insolation on the TAC. An anti-correlation between local summer insolation and TAC is known in ice core records from East Antarctica during the last 400 000 years (Raynaud et al., 2007), and the insolation signal is also found in Greenlandic TAC (NGRIP; Eicher et al., 2016). NEEM community members (2013) estimates (based on data from the Holocene optimum) that the summer insolation could account for 50 % of the observed Eemian TAC changes at NEEM.
Nevertheless, the possibility of a melt-induced reduction of TAC should be considered for the interpretation of Eemian air content to estimate ice surface elevation changes. An early interpretation of the first Greenland ice cores (Camp Century, Dye-3) suggested an extreme scenario for Eemian Greenland with extensive melt and a much smaller ice sheet, leading to a sea level rise of 6 m (Koerner, 1989). However, this scenario was rejected by later ice core studies showing evidence of Eemian ice (especially NGRIP and NEEM; North Greenland Ice Core Project members et al., 2004; NEEM community members, 2013). Furthermore, GRIP TAC measurements (Raynaud, 1999) have been interpreted as evidence for the elevation of the summit sites having remained above 3000 m of altitude during the Eemian and GRIP deuterium excess measurements remain in the normal range during the Eemian (Landais et al., 2003). However, this last interpretation can be challenged by measurements of a NEEM Holocene melt layer, suggesting that the melt layer mainly influences TAC and CH4 observations, while other variables like deuterium excess may be less influenced by melt (NEEM community members, 2013).
The climate simulations show surface melt at all deep ice core locations and at the Agassiz ice cap under 125 ka climate conditions (Figs. 4 and A2; orange column). Even locations near the summit of Greenland (GRIP, GISP2, and NGRIP) show a few melt days per year on average (defined as >8 mm d−1) during these warm Eemian simulations. NEEM, the ice core location with the longest Eemian record, shows ∼8 melt days per year. While the presence of Eemian surface melt at NEEM was acknowledged previously (NEEM community members, 2013), the lower TAC observations at GRIP, GISP2, and NGRIP could as well be related to Eemian surface melt, rather than stable or higher elevations.
Finally, it should be emphasized that a robust estimate of Eemian Greenland surface melt is challenging to obtain with a single climate model. Ideally, there should be an ensemble of climate models to explore model biases and uncertainties. However, as pointed out earlier in this discussion, there are several reasons why the presented climate simulations could be on the lower end of available climate model in terms of the amount of simulated Eemian melt. It is likely that there are other climate models which show more extensive Eemian surface melt.
In the future, an analysis of individual or ensemble Eemian climate simulations would benefit from a comparison of the observed extreme melt event in 2012 (and similar events in the recent past) with simulated extreme Eemian melt events. Relationships in the Eemian simulations between air temperature and local wind patterns and the resulting simulated melt could be analyzed and used to identify specific weather patterns leading to high surface melt in the simulations (e.g., similar analysis performed by Neff et al., 2014; Keegan et al., 2014; Fettweis et al., 2013b; Tedesco and Fettweis, 2020).
Using regional climate simulations (including a full surface energy balance), this study shows surface melt at all Greenland ice core locations during the Eemian interglacial period (e.g., GRIP, GISP2: ∼60 ; NGRIP: ∼150 ). The amount of refreezing exceeds 25 % of the annual accumulation at the summit of Greenland (GRIP, GISP2) and reaches values as high as 90 % at less central locations like Dye-3 and EGRIP. The simulated air pressure, temperature, and refreezing are used to estimate Eemian TAC and high melt rates could explain the low corresponding ice core TAC observations. This is true even though the discussed simulations could show conservative melt estimates (several potentially melt-increasing processes are neglected). Therefore, the possibility of widespread surface melt should be considered for the interpretation of Greenlandic total air content records (as an elevation proxy) from warm periods such as the Eemian interglacial period. Finally, a robust map of Eemian melt estimates in Greenland in combination with accumulation patterns could be used to identify potential future ice cores sites on Greenland. Such a procedure would increase the chances of finding Eemian ice influenced by a minimum number of melt layers. These sites will have relatively high accumulation combined with low surface melt.
The source code of MAR, version 3.7, is available on the MAR website: https://mar.cnrs.fr (last access: 20 January 2021, CNRS, 2021). An description on how to retrieve the source code is given in the download section of the MAR website: https://mar.cnrs.fr/index.php?option_smdi=presentation&idm=10, last access: 20 January 2021.
The Eemian MAR simulations are available from the corresponding author upon request. MEaSUREs Greenland Surface Melt Daily 25 km EASE-Grid 2.0, version 1 is freely available at https://doi.org/10.5067/MEASURES/CRYOSPHERE/nsidc-0533.001 (Mote, 2014). For more information and to request the T19Hmelt data (Fettweis et al., 2011, https://doi.org/10.5194/tc-5-359-2011), please contact Xavier Fettweis (email@example.com). For more information and to request the collection of Greenland shallow ice core and weather station data (Faber, 2016), please contact Anne-Katrine Faber (firstname.lastname@example.org). The TAC observations at NEEM are freely available at 10.1038/nature11789 (NEEM community members, 2013). The GRIP TAC is freely available at https://doi.org/10.1594/PANGAEA.55086 (Raynaud, 1999). The GISP2 TAC is freely available as a Supplement to Yau et al. (2016) (https://doi.org/10.1073/pnas.1524766113). The NGRIP TAC is freely available at https://www.ncdc.noaa.gov/paleo/study/20569 (last access: 27 November 2020, Eicher et al., 2016, https://doi.org/10.5194/cp-12-1979-2016). The Dye-3 TAC data was extracted from Fig. 4 in Herron and Langway (1987). For more information and to request the extracted data please contact Sindhu Vudayagiri (email@example.com) or Thomas Blunier (firstname.lastname@example.org).
AP and BMV designed the study with contributions from KHN. Dye-3 total air content data were extracted from Fig. 4 in Herron and Langway (1987) by SV. AP made the figures and wrote the text with input from BMV, KHN, SV, and TB.
The authors declare that they have no conflict of interest.
We thank Chuncheng Guo for performing the Eemian NorESM simulations and Sébastien Le clec'h for downscaling the NorESM simulations with the regional model MAR. Furthermore, we would like to thank Anne-Katrine Faber for valuable discussions and providing the shallow ice core data she compiled during her PhD. We also thank Xavier Fettweis for providing the T19Hmelt data. Furthermore, we very much thank the editor Eric Wolff and two anonymous referees for their comments and suggestions, which significantly improved the manuscript.
This research has been supported by the European Research Council under the European Community's Seventh Framework Programme (FP7/2007-2013) (ICE2ICE (grant no. 610055)).
This paper was edited by Eric Wolff and reviewed by two anonymous referees.
Alley, R. B. and Anandakrishnan, S.: Variations in melt-Iayer frequency in the GISP2 ice core: implications for Holocene summer temperatures in central Greenland, available at: https://www.igsoc.org/annals/21/igs_annals_vol21_year1995_pg64-70.pdf (last access: 13 January 2021), 1995. a
Alley, R. B. and Koci, B. R.: Ice-Core Analysis at Site A, Greenland: Preliminary Results, Ann. Glaciol., 10, 1–4, https://doi.org/10.3189/S0260305500004067, 1988. a
Brun, E., David, P., Sudul, M., and Brunot, G.: A numerical model to simulate snow-cover stratigraphy for operational avalanche forecasting, J. Glaciol., 38, 13–22, https://doi.org/10.3189/S0022143000009552, 1992. a
CAPE Last Interglacial Project Members: Last Interglacial Arctic warmth confirms polar amplification of climate change, Quaternary Sci. Rev., 25, 1383–1400, https://doi.org/10.1016/j.quascirev.2006.01.033, 2006. a, b
Capron, E., Govin, A., Stone, E. J., Masson-Delmotte, V., Mulitza, S., Otto-Bliesner, B., Rasmussen, T. L., Sime, L. C., Waelbroeck, C., and Wolff, E. W.: Temporal and spatial structure of multi-millennial temperature changes at high latitudes during the Last Interglacial, Quaternary Sci. Rev., 103, 116–133, https://doi.org/10.1016/j.quascirev.2014.08.018, 2014. a, b
CNRS: MAR code, version 3.7, available at: http://mar.cnrs.fr/, last access: 20 January 2021. a
Eicher, O., Baumgartner, M., Schilt, A., Schmitt, J., Schwander, J., Stocker, T. F., and Fischer, H.: Climatic and insolation control on the high-resolution total air content in the NGRIP ice core, Clim. Past, 12, 1979–1993, https://doi.org/10.5194/cp-12-1979-2016, 2016 (available at: https://www.ncdc.noaa.gov/paleo-search/study/20569, last access: 27 November 2020). a, b, c, d, e
Faber, A.: Isotopes in Greenland precipitation: Isotope-enabled AGCM modelling and a new Greenland database of observations and ice core measurements, PhD thesis, Centre for Ice and Climate, University of Copenhagen, 2016. a, b
Fettweis, X.: Reconstruction of the 19792006 Greenland ice sheet surface mass balance using the regional climate model MAR, The Cryosphere, 1, 21–40, https://doi.org/10.5194/tc-1-21-2007, 2007. a
Fettweis, X., Tedesco, M., van den Broeke, M., and Ettema, J.: Melting trends over the Greenland ice sheet (19582009) from spaceborne microwave data and regional climate models, The Cryosphere, 5, 359–375, https://doi.org/10.5194/tc-5-359-2011, 2011. a, b, c
Fettweis, X., Franco, B., Tedesco, M., van Angelen, J. H., Lenaerts, J. T. M., van den Broeke, M. R., and Gallée, H.: Estimating the Greenland ice sheet surface mass balance contribution to future sea level rise using the regional atmospheric climate model MAR, The Cryosphere, 7, 469–489, https://doi.org/10.5194/tc-7-469-2013, 2013a. a
Fettweis, X., Hanna, E., Lang, C., Belleflamme, A., Erpicum, M., and Gallée, H.: Brief communication “Important role of the mid-tropospheric atmospheric circulation in the recent surface melt increase over the Greenland ice sheet”, The Cryosphere, 7, 241–248, https://doi.org/10.5194/tc-7-241-2013, 2013b. a
Fettweis, X., Box, J. E., Agosta, C., Amory, C., Kittel, C., Lang, C., van As, D., Machguth, H., and Gallée, H.: Reconstructions of the 19002015 Greenland ice sheet surface mass balance using the regional climate MAR model, The Cryosphere, 11, 1015–1033, https://doi.org/10.5194/tc-11-1015-2017, 2017. a
Fisher, D., Zheng, J., Burgess, D., Zdanowicz, C., Kinnard, C., Sharp, M., and Bourgeois, J.: Recent melt rates of Canadian arctic ice caps are the highest in four millennia, Global Planet. Change, 84–85, 3–7, https://doi.org/10.1016/j.gloplacha.2011.06.005, 2012. a
Gallée, H. and Duynkerke, P. G.: Air-snow interactions and the surface energy and mass balance over the melting zone of west Greenland during the Greenland Ice Margin Experiment, J. Geophys. Res.-Atmos., 102, 13813–13824, https://doi.org/10.1029/96JD03358, 1997. a
Gallée, H. and Pettré, P.: Dynamical Constraints on Katabatic Wind Cessation in Adélie Land, Antarctica, J. Atmos. Sci., 55, 1755–1770, https://doi.org/10.1175/1520-0469(1998)055<1755:DCOKWC>2.0.CO;2, 1998. a
Guo, C., Bentsen, M., Bethke, I., Ilicak, M., Tjiputra, J., Toniazzo, T., Schwinger, J., and Otterå, O. H.: Description and evaluation of NorESM1-F: a fast version of the Norwegian Earth System Model (NorESM), Geosci. Model Dev., 12, 343–362, https://doi.org/10.5194/gmd-12-343-2019, 2019. a
Herron, S. L. and Langway, C. C.: Derivation of paleoelevations from total air content of two deep Greenland ice cores, in: The Physical Basis of Ice Sheet Modelling, Proceedings of the Vancouver Symposium, August 1987, IAHS Publ. no. 170, 14 pp., 1987. a, b
Keegan, K. M., Albert, M. R., McConnell, J. R., and Baker, I.: Climate change and forest fires synergistically drive widespread melt events of the Greenland Ice Sheet, P. Natl. Acad. Sci. USA, 111, 7964–7967, https://doi.org/10.1073/pnas.1405397111, 2014. a, b, c
Kessler, E.: On the Distribution and Continuity of Water Substance in Atmospheric Circulations, in: On the Distribution and Continuity of Water Substance in Atmospheric Circulations, American Meteorological Society, Boston, MA, Meteorological Monographs, 1–84, https://doi.org/10.1007/978-1-935704-36-2_1, 1969. a
Koerner, R. M.: Ice Core Evidence for Extensive Melting of the Greenland Ice Sheet in the Last Interglacial, Science, 244, 964–968, https://doi.org/10.1126/science.244.4907.964, 1989. a
Landais, A., Chappellaz, J., Delmotte, M., Jouzel, J., Blunier, T., Bourg, C., Caillon, N., Cherrier, S., Malaizé, B., Masson Delmotte, V., Raynaud, D., Schwander, J., and Steffensen, J. P.: A tentative reconstruction of the last interglacial and glacial inception in Greenland based on new gas measurements in the Greenland Ice Core Project (GRIP) ice core, J. Geophys. Res., 108, 4563, https://doi.org/10.1029/2002JD003147, 2003. a, b, c
Landais, A., Masson-Delmotte, V., Capron, E., Langebroek, P. M., Bakker, P., Stone, E. J., Merz, N., Raible, C. C., Fischer, H., Orsi, A., Prié, F., Vinther, B., and Dahl-Jensen, D.: How warm was Greenland during the last interglacial period?, Clim. Past, 12, 1933–1948, https://doi.org/10.5194/cp-12-1933-2016, 2016. a, b
Lecavalier, B. S., Fisher, D. A., Milne, G. A., Vinther, B. M., Tarasov, L., Huybrechts, P., Lacelle, D., Main, B., Zheng, J., Bourgeois, J., and Dyke, A. S.: High Arctic Holocene temperature record from the Agassiz ice cap and Greenland ice sheet evolution, P. Natl. Acad. Sci. USA, 114, 5952–5957, https://doi.org/10.1073/pnas.1616287114, 2017. a
Le clec'h, S., Charbit, S., Quiquet, A., Fettweis, X., Dumas, C., Kageyama, M., Wyard, C., and Ritz, C.: Assessment of the Greenland ice sheet–atmosphere feedbacks for the next century with a regional atmospheric model coupled to an ice sheet model, The Cryosphere, 13, 373–395, https://doi.org/10.5194/tc-13-373-2019, 2019. a, b
Lin, Y.-L., Farley, R. D., and Orville, H. D.: Bulk Parameterization of the Snow Field in a Cloud Model, J. Clim. Appl. Meteorol., 22, 1065–1092, https://doi.org/10.1175/1520-0450(1983)022<1065:BPOTSF>2.0.CO;2, 1983. a
Martinerie, P., Raynaud, D., Etheridge, D. M., Barnola, J.-M., and Mazaudier, D.: Physical and climatic parameters which influence the air content in polar ice, Earth Planet. Sc. Lett., 112, 1–13, https://doi.org/10.1016/0012-821X(92)90002-D, 1992. a
Martinerie, P., Lipenkov, V., Raynaud, D., Chappellaz, J., Barkov, N. I., and Lorius, C.: Air content paleo record in the Vostok ice core (Antarctica): A mixed record of climatic and glaciological parameters, J. Geophys. Res.-Atmos., 99, 10565–10576, https://doi.org/10.1029/93JD03223, 1994. a
Merz, N., Born, A., Raible, C. C., Fischer, H., and Stocker, T. F.: Dependence of Eemian Greenland temperature reconstructions on the ice sheet topography, Clim. Past, 10, 1221–1238, https://doi.org/10.5194/cp-10-1221-2014, 2014a. a
Merz, N., Gfeller, G., Born, A., Raible, C. C., Stocker, T. F., and Fischer, H.: Influence of ice sheet topography on Greenland precipitation during the Eemian interglacial, J. Geophys. Res.-Atmos., 119, 10749–10768, https://doi.org/10.1002/2014JD021940, 2014b. a
Morcrette, J.-J., Barker, H. W., Cole, J. N. S., Iacono, M. J., and Pincus, R.: Impact of a New Radiation Package, McRad, in the ECMWF Integrated Forecasting System, Mon. Weather Rev., 136, 4773–4798, https://doi.org/10.1175/2008MWR2363.1, 2008. a
Mote, T. L.: MEaSUREs Greenland Surface Melt Daily 25 km EASE-Grid 2.0, Version 1, Greenland subset, NASA National Snow and Ice Data Center Distributed Active Archive Center, Boulder, Colorado USA, https://doi.org/10.5067/MEASURES/CRYOSPHERE/nsidc-0533.001, 2014. a, b
NEEM community members: Eemian interglacial reconstructed from a Greenland folded ice core, Nature, 493, 489–494, https://doi.org/10.1038/nature11789, 2013. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o
Neff, W., Compo, G. P., Ralph, F. M., and Shupe, M. D.: Continental heat anomalies and the extreme melting of the Greenland ice surface in 2012 and 1889, J. Geophys. Res.-Atmos., 119, 6520–6536, https://doi.org/10.1002/2014JD021470, 2014. a, b
Nghiem, S. V., Hall, D. K., Mote, T. L., Tedesco, M., Albert, M. R., Keegan, K., Shuman, C. A., DiGirolamo, N. E., and Neumann, G.: The extreme melt across the Greenland ice sheet in 2012: EXTREME MELT ACROSS GREENLAND ICE SHEET, Geophys. Res. Lett., 39, L20502, https://doi.org/10.1029/2012GL053611, 2012. a
North Greenland Ice Core Project members: Andersen, K. K., Azuma, N., Barnola, J.-M., Bigler, M., Biscaye, P., Caillon, N., Chappellaz, J., Clausen, H. B., Dahl-Jensen, D., Fischer, H., Flückiger, J., Fritzsche, D., Fujii, Y., Goto-Azuma, K., Grønvold, K., Gundestrup, N. S., Hansson, M., Huber, C., Hvidberg, C. S., Johnsen, S. J., Jonsell, U., Jouzel, J., Kipfstuhl, S., Landais, A., Leuenberger, M., Lorrain, R., Masson-Delmotte, V., Miller, H., Motoyama, H., Narita, H., Popp, T., Rasmussen, S. O., Raynaud, D., Rothlisberger, R., Ruth, U., Samyn, D., Schwander, J., Shoji, H., Siggard-Andersen, M.-L., Steffensen, J. P., Stocker, T., Sveinbjörnsdóttir, A. E., Svensson, A., Takata, M., Tison, J.-L., Thorsteinsson, T., Watanabe, O., Wilhelms, F., and White, J. W. C.: High-resolution record of Northern Hemisphere climate extending into the last interglacial period, Nature, 431, 147–151, https://doi.org/10.1038/nature02805, 2004. a
Otto-Bliesner, B. L., Rosenbloom, N., Stone, E. J., McKay, N. P., Lunt, D. J., Brady, E. C., and Overpeck, J. T.: How warm was the last interglacial? New model–data comparisons, Philos. T. R. Soc. A, 371, 20130097, https://doi.org/10.1098/rsta.2013.0097, 2013. a, b
Pfeffer, W. T. and Humphrey, N. F.: Formation of ice layers by infiltration and refreezing of meltwater, Ann. Glaciol., 26, 83–91, https://doi.org/10.3189/1998AoG26-1-83-91, 1998. a
Plach, A., Nisancioglu, K. H., Le clec'h, S., Born, A., Langebroek, P. M., Guo, C., Imhof, M., and Stocker, T. F.: Eemian Greenland SMB strongly sensitive to model choice, Clim. Past, 14, 1463–1485, https://doi.org/10.5194/cp-14-1463-2018, 2018. a, b
Plach, A., Nisancioglu, K. H., Langebroek, P. M., Born, A., and Le clec'h, S.: Eemian Greenland ice sheet simulated with a higher-order model shows strong sensitivity to surface mass balance forcing, The Cryosphere, 13, 2133–2148, https://doi.org/10.5194/tc-13-2133-2019, 2019. a, b
Raynaud, D.: GRIP total air content, PANGAEA, https://doi.org/10.1594/PANGAEA.55086, 1999. a, b, c
Raynaud, D., Chappellaz, J., Ritz, C., and Martinerie, P.: Air content along the Greenland Ice Core Project core: A record of surface climatic parameters and elevation in central Greenland, J. Geophys. Res.-Oceans, 102, 26607–26613, https://doi.org/10.1029/97JC01908, 1997. a, b, c
Raynaud, D., Lipenkov, V., Lemieux-Dudon, B., Duval, P., Loutre, M.-F., and Lhomme, N.: The local insolation signature of air content in Antarctic ice. A new step toward an absolute dating of ice records, Earth Planet. Sc. Lett., 261, 337–349, https://doi.org/10.1016/j.epsl.2007.06.025, 2007. a, b
Sander, R.: Compilation of Henry's law constants (version 4.0) for water as solvent, Atmos. Chem. Phys., 15, 4399–4981, https://doi.org/10.5194/acp-15-4399-2015, 2015. a, b
Souchez, R., Bouzette, A., Clausen, H. B., Johnsen, S. J., and Jouzel, J.: A stacked mixing sequence at the base of the Dye 3 Core, Greenland, Geophys. Res. Lett., 25, 1943–1946, https://doi.org/10.1029/98GL01411, 1998. a
Tedesco, M. and Fettweis, X.: Unprecedented atmospheric conditions (19482019) drive the 2019 exceptional melting season over the Greenland ice sheet, The Cryosphere, 14, 1209–1223, https://doi.org/10.5194/tc-14-1209-2020, 2020. a
Van de Berg, W. J., van den Broeke, M., Ettema, J., van Meijgaard, E., and Kaspar, F.: Significant contribution of insolation to Eemian melting of the Greenland ice sheet, Nat. Geosci., 4, 679–683, https://doi.org/10.1038/ngeo1245, 2011. a
Vinther, B. M., Clausen, H. B., Fisher, D. A., Koerner, R. M., Johnsen, S. J., Andersen, K. K., Dahl-Jensen, D., Rasmussen, S. O., Steffensen, J. P., and Svensson, A. M.: Synchronizing ice cores from the Renland and Agassiz ice caps to the Greenland Ice Core Chronology, J. Geophys. Res., 113, D08115, https://doi.org/10.1029/2007JD009143, 2008. a
Vinther, B. M., Buchardt, S. L., Clausen, H. B., Dahl-Jensen, D., Johnsen, S. J., Fisher, D. A., Koerner, R. M., Raynaud, D., Lipenkov, V., Andersen, K. K., Blunier, T., Rasmussen, S. O., Steffensen, J. P., and Svensson, A. M.: Holocene thinning of the Greenland ice sheet, Nature, 461, 385–388, https://doi.org/10.1038/nature08355, 2009. a, b
Willerslev, E., Cappellini, E., Boomsma, W., Nielsen, R., Hebsgaard, M. B., Brand, T. B., Hofreiter, M., Bunce, M., Poinar, H. N., Dahl-Jensen, D., Johnsen, S., Steffensen, J. P., Bennike, O., Schwenninger, J.-L., Nathan, R., Armitage, S., Hoog, C.-J. d., Alfimov, V., Christl, M., Beer, J., Muscheler, R., Barker, J., Sharp, M., Penkman, K. E. H., Haile, J., Taberlet, P., Gilbert, M. T. P., Casoli, A., Campani, E., and Collins, M. J.: Ancient Biomolecules from Deep Ice Cores Reveal a Forested Southern Greenland, Science, 317, 111–114, https://doi.org/10.1126/science.1141758, 2007. a
Yau, A. M., Bender, M. L., Robinson, A., and Brook, E. J.: Reconstructing the last interglacial at Summit, Greenland: Insights from GISP2, P. Natl. Acad. Sci. USA, 113, 9710–9715, https://doi.org/10.1073/pnas.1524766113, 2016. a, b
Yin, Q. Z. and Berger, A.: Insolation and CO2 contribution to the interglacial climate before and after the Mid-Brunhes Event, Nat. Geosci., 3, 243–246, https://doi.org/10.1038/ngeo771, 2010. a, b, c