Climate Noise Influences Ice Sheet Mean State

Evidence from proxy records indicates that millennial‐scale abrupt climate shifts, called Dansgaard‐Oeschger events, happened during past glacial cycles. Various studies have been conducted to uncover the physical mechanism behind them, based on the assumption that climate mean state determines the variability. However, our study shows that the Dansgaard‐Oeschger events can regulate the mean state of the Northern Hemisphere ice sheets. Sensitivity experiments show that the simulated mean state is influenced by the amplitude of the climatic noise. The most likely cause of this phenomenon is the nonlinear response of the surface mass balance to temperature. It could also cause the retreat processes to be faster than the buildup processes within a glacial cycle. We propose that the climate variability hindered ice sheet development and prevented the Earth system from entering a full glacial state from Marine Isotope Stage 4 to Marine Isotope Stage 3 about 60,000 years ago.


Introduction
Millennial-scale abrupt shifts between cold and warm states, called Dansgaard-Oeschger (DO) events (Dansgaard et al., 1993), are observed in Greenland ice core records through the last glacial cycle. The DO events began with a rapid warming of 8 to 15 • C within a few decades, followed by gradual cooling before a jump back to a cold state (Huber et al., 2006). Spectral analysis shows a typical recurrence time of about 1,470 years (Schulz, 2002). However, whether DO events are periodic or noise induced is still debated (Alley et al., 2001;Benzi et al., 1982;Ditlevsen et al., 2007;Ganopolski & Rahmstorf, 2002). In Antarctic ice cores, millennial-scale oscillations with nearly opposite phase compared to Greenland are found, showing that this phenomenon is global in scale and might be related to the Atlantic ocean heat transport (EPICA Community Members, 2006;Dima et al., 2018). Other proxy records from different areas (e.g., North Atlantic Ocean sediments, European pollen, and Chinese loess) also contain these signals (Bond et al., 1993;Porter, 2001;Woillard & Mook, 1982). Despite large research efforts, no consistent theory of the physical mechanism of the cause of DO events is widely accepted (Ganopolski & Rahmstorf, 2001;Petersen et al., 2013;Zhang, Lohmann, et al., 2014). Global and regional change in surface temperature variability on decadal to millennial timescales (here called "climate noise") declined when the Earth system transitioned from a glacial to an interglacial condition (Rehfeld et al., 2018;Shao & Ditlevsen, 2016). The Northern Hemisphere ice sheets were restricted to Greenland during interglacials and advanced further south to cover large portions of North America and northern Europe during glacials. Larger climatic variability (DO events) happened when ice sheets were intermediate in size, indicating that DO cycles are likely more closely linked to ice sheet evolution than to other potential drivers within the Earth system (Rahmstorf, 1995). Various studies applied numerical climate models to uncover the mechanism behind DO events, in terms of freshwater perturbations or other potential drivers (Bond & Lotti, 1995;Broecker, 1994;Clark et al., 2001;Ganopolski & Rahmstorf, 2001;Rahmstorf, 1995;Zhang, Lohmann, et al., 2014;Zhang, Prange, et al., 2014). What is common to all these theories is that climate mean state determines climate variability. In the majority of these studies, the simulated climate states (here called "climate mean state") are prescribed with predefined insolation, greenhouse gases, and ice sheet configuration. Under this premise, the thermohaline circulation exhibits several modes, which could be linked to millennial-scale climate variability. With massive discharge of freshwater from ice sheets, the North Atlantic meridional overturning circulation (AMOC) is weakened and leads to a reduction of heat transport to the Northern Hemisphere, constituting a cold state (stadial). With a reduction of freshwater discharge, the AMOC is strengthened, constituting a warm state (interstadial).

10.1029/2019GL083717
However, how the climate variability regulates the mean state has not been fully investigated (as is shown in Figure 1 in Timmermann & Lohmann, 2000). Cessi (1994) pointed out that the magnitude of stochastic noise is capable of modifying the equilibrium state of the ocean circulation in a Stommel two-box model (Stommel, 1961). Timmermann and Lohmann (2000) confirmed that noise-induced transitions could be produced in a simplified ocean circulation model. Marshall et al. (2000) and Charbit et al. (2002) found that large century-to millennial-scale climate variability strongly impacts the evolution of simulated ice sheets. In a recent study, Mikkelsen et al. (2018) showed that the incorporation of interannual temperature fluctuations results in a smaller equilibrium Greenland ice sheet volume. These studies suggest that climate variability could in turn drive mean state changes.
In this paper, we first investigate the influence of DO-like climate fluctuations on the Northern Hemisphere ice sheet evolution during glacial-interglacial cycles. By using a glacial index method, we show how the ice sheet mean state is regulated by the fluctuations. Second, we conducted sensitivity experiments to see how the amplitude of climate noise affects the simulated ice sheets. Finally, a physical interpretation of our results is presented.

Methods
Abrupt temperature shifts imprinted in Greenland ice cores between cold (stadials) and warm (interstadials) states range from 8 to 15 • C in amplitude (Huber et al., 2006). The two distinct climate states are likely related to the AMOC, where interstadial (stadial) climate with strong (weak) overturning circulation is close to the preindustrial (PI; the Last Glacial Maximum [LGM]) condition (Ganopolski & Rahmstorf, 2001;Zhang, Prange, et al., 2014). From this perspective, a glacial index method is suitable for investigating the influence of these climate fluctuations on ice sheet evolution, since it synthesizes the necessary boundary conditions (Charbit et al., 2002(Charbit et al., , 2007Greve et al., 1999;Marshall et al., 2000Marshall et al., , 2002Rodgers et al., 2004;Tarasov & Peltier, 2004;Zweck & Huybrechts, 2005). Note that millennial-scale climate fluctuations cannot be fully represented with this method, since orbital-scale signals may be misrepresented. However, matching the simulated ice sheets to the geological constraints under different time or spatial scales is not the main focus of this study. An advantage of the glacial index method is that it is computationally fast and therefore is the best suited method for the purpose of our study. In our study, Northern Hemisphere ice sheet simulations are driven by a glacial index I combined with monthly near-surface air temperature and precipitation fields at two time slices, the PI warm state (I = 0) and the LGM cold state (I = 1; Figure 1a). The PI and the LGM climate states are taken from equilibrium simulations using the General Circulation Model COSMOS (Figures S1 and S2 in the supporting information; Stepanek & Lohmann, 2012;Zhang et al., 2013). The temperature and precipitation fields at other time periods are scaled by the index I and linearly interpolated between the LGM and PI states (see supporting information). The isotope records from the Greenland ice cores are often used as a glacial index, but they only date back to the Last Interglacial (127,000 years BP) due to the relatively high accumulation rate over Greenland. To extend our simulation, we use an 800,000-year synthetic record of Greenland climate variability from Barker et al. (2011), which is based on the Antarctica ice core records and a thermal bipolar seesaw model (Figure 1a, black line). The Greenland synthetic record exhibits a larger variability within the last two glacial cycles than in the earlier cycles. This is likely because the signals are attenuated in the older and thinner ice layers.
For our ice sheet simulation, we use the three-dimensional thermomechanically coupled Parallel Ice Sheet Model (version 0.7.3 ; The PISM authors, 2016). The stress balance in the ice is computed with a combination of the shallow ice and shallow shelf approximations (Bueler & Brown, 2009). The spatial domain covers the Northern Hemisphere with a spatial resolution of 40 × 40 km (same domain as in Figures S1 and S2). Solid Earth deformation from glacial isostatic adjustment is calculated using the Lingle and Clark method (Bueler et al., 2007;Lingle & Clark, 1985). The climatic mass balance (hereafter called "surface mass balance" [SMB]) is computed based on the semiempirical positive degree day (PDD) scheme (Calov & Greve, 2005;Hock, 2005;Cogley et al., 2011;Reeh, 1991), which is calibrated with modern observations (Ritz, 1997).
The assumption in this method is that the melt rate of snow and ice is proportional to the sum of the positive surface air temperature values over the year. Although the PDD scheme may underestimate the contribution of shortwave radiation, it requires only two variables as input (temperature and precipitation) and is still suitable for the purpose of our study (Niu et al., 2019). More details of the model parameters can be found in the supporting information   Ritz, 1997;Rohling et al., 2014;Stepanek & Lohmann, 2012;Winkelmann et al., 2011;Zhang et al., 2013). The model setup is the same as in Niu et al. (2019), where a more detailed discussion of the model performance is presented.

Results
We conducted two experiments. One simulation was forced by the Greenland 800,000-year synthetic record which contains high-frequency DO signals (GL_hf, Figure 1a, black). In the other simulation, the DO signals and other high-frequency variability were eliminated by applying a 5,000-year running mean filter (GL_lf, Figure 1a, red). The results show that the ice sheets grow larger when the DO signals are absent, while the mean forcing states are the same (Figure 1b). The difference of the ice sheet volume between the two experiments, expressed as the equivalent contribution to eustatic sea level change (SLE), can be up to 60 m. The simulated SLE curves are compared to a reconstructed relative sea level record (Rohling et al., 2014). The root-mean-square deviation between the model and the reconstruction is 40.3 m for GL_hf and 52.0 m for GL_lf. During the last two glacial cycles, the consistency between the GL_hf SLE and the reconstruction increased significantly, with root-mean-square deviation value of 21.7 m (42.0 m for GL_lf). A coherence test shows strong association between the simulated SLE and the climate forcing time series within the millennial timescale ( Figure S3). Note that the DO signals are stronger during the last two glacial cycles than for older times. We propose that the DO-like fluctuations significantly inhibit the ice sheet buildup processes. The relationship between the SLE difference of the two simulations and the amplitude of the climate variability (standard deviation of the index) is shown in Figure S4. The simulated ice sheet volume is reduced in experiment GL_hf, which has large climate variability. The larger the variability, the larger the maximum SLE difference (|GL_hf-GL_lf|) becomes.
The spatial distribution of the ice sheets also differs between the two experiments. Figure 1c shows the simulated ice sheet thickness pattern at 158,000 years BP, when the SLE difference (|GL_hf-GL_lf|) is ∼40 m (Figure 1b). The North American ice sheets and Eurasian ice sheets extend much further south to around 40 • N with lower climate variability (GL_lf). In experiment GL_hf, the North American ice sheets only reach the southern margin of the Hudson Bay. The difference in ice sheet thickness can be more than 3,000 m in southern North American ice sheets (Figure 1c). Our results exhibit two distinct simulated ice sheet mean states during the ice sheet buildup stage. Ice sheets reach maximal full glacial conditions much earlier when the DO signals are absent.
To assess the influence of climate variability on the mean state of the ice sheets, we carried out sensitivity experiments. The isotopic record used for GL_hf and GL_lf is replaced by a cosine-based index. It ranges from 0 to 1 with a period of 120,000 years ( Figure 2). For different sensitivity experiments, white noise with different standard deviation (from 0 to 0.4 with 0.05 interval) was added onto the index. The simulations run for 120,000 years. The other parameters of the model setup are the same as before.
The results show hysteresis in all the simulated ice sheets. The glacial maximum occurs around year 70,000, instead of year 60,000 when temperature forcing is the lowest (Figure 2i). The ice sheet retreat period (from year 70,000 to year 120,000) is shorter than the buildup periods (from year 0 to year 70,000). All simulations end with similar SLE (about 8 m) at 120,000 years. With different values, the simulated maximum SLE differs. The index with larger variability results in a smaller SLE. The simulated ice sheet with 0.4 is half the size of the one with 0 (black line). From the sensitivity experiments, we confirm that the simulated ice sheet mean state is influenced by the amplitude of climate noise. Larger amplitude variability slows down the ice sheet buildup process and results in smaller ice sheets at the glacial maximum.

Discussion
Our results show that the simulated ice sheets tend to be smaller when influenced by climate noise. One possible reason could be the nonlinear response of SMB to temperature. Such nonlinearity could cause an asymmetry of the change in SMB in response to positive and negative temperature anomalies. We demonstrate this in a theoretical way as follows.
We assume that the annual temperature cycle in one location follows a cosine function (Reeh, 1991) where T mean and T July are the annual mean (−10 • C) and July mean (0 • C) surface air temperatures respectively. The time t is in days and A is 1 year (365 days). White noise (T noise ) with a standard deviation of 3 • C is added to account for synoptic variations (Figure 3b, red ; Krebs-Kanzow et al., 2018). The potential surface melt rate over the year (M A ) can be computed with the semiempirical PDD method (Reeh, 1991). The related degree day factor value f PDD is a freely tuned parameter to match the observed surface energy balance under different situations. Assuming a constant precipitation (P) of 0.002 m water equivalent (w.e.)/day and for snow the SMB over 1 year (SMB A ) can be estimated with a refreezing rate (f refreeze = 0.6) as follows: Temperature changes are represented by shifting the annual cycle temperature curve up and down. The corresponding annual SMB (SMB A , notated with SMB) can also be calculated (Figure 3a, black line, P = 0.002 m/day, f PDD = 0.003 m· • C −1 ·day −1 ). We find a nonlinear relationship between the annual SMB and the annual mean temperature. To be more specific, SMB is a concave downward function of T, that is, 2 (SMB) T 2 < 0. The larger the increase of temperature, the larger the loss of SMB. We also varied the P and f PDD values to test the sensitivity of the calculated SMB (Figure 3a). When we vary f PDD from 0.003 m· • C −1 ·day −1 for snow to 0.009 m· • C −1 ·day −1 for ice (Ritz, 1997), the resultant SMB ranges from −4.38 to −13.14 m w.e./year with T mean = 10 • C. However, the SMB has a linear response from 0 to ∼1 m w.e./year when precipitation increases from 0 to 1 m/year. Around the ice sheet margins of the nonmarine terminating Northern Hemisphere ice sheets, the annual SMB is negative at lower altitudes and positive at higher altitudes. Now we focus on the black line in Figure 3a and assume that a point is located around the ice sheet margin with annual mean temperature of −6 • C and SMB of 0 m w.e./year (intersection of dotted lines). The change of SMB (| SMB|) responding to positive or negative temperature anomalies (| T|) is shown in Figure 3c. Focusing on the red line (blue line), when the positive (negative) temperature anomaly increases, the SMB decreases (increases) and the rate of change increases (decreases). As a result, the change of SMB due to a warm anomaly cannot be balanced by an identical cold anomaly (d| SMB| ≠ 0, black line). In our case, a temperature variance of ±5 • C  (1). White noise with a standard deviation of 3 • C is added on to the curve to account for the synoptic variations. The red line is with an annual mean temperature of −10 • C. The pink lines are shifted by ±10 • C from the red line. (c) The relation between the annual mean temperature anomalies | T| and the change of the resultant surface mass balance | SMB|, based on the intersection of dashed lines in (a). The red line is with the warm anomaly, the blue line is with the cold anomaly, and the black is the difference between the red line and the blue line.
can cause a potential reduction of 0.4 m w.e./year of SMB (Figure 3c). In summary, for a melt-dominated process, the increase of surface mass budget responding to negative temperature anomaly cannot compensate the decrease of surface mass budget responding to positive temperature anomaly. Thus, climate noise leads to a decrease of mean SMB.
In reality, the surface mass budget is influenced not only by temperature but also by precipitation change. The exponential relationship between water vapor saturation pressure and temperature can result in precipitation increase by 5-7% per degree increase in temperature (Huybrechts, 2002). As is shown in Figure S2e and S2f, more precipitation is simulated in most of the Northern Hemisphere at PI than at the LGM. Such an adjustment could increase snow accumulation in winter which competes with the higher melt in the summer. This effect could be significant especially for accumulation-dominated areas, where SMB becomes a concave upward function of T, that is, 2 (SMB) T 2 > 0. In this case, climate noise would lead to an increase of mean SMB. For example, in some regions in Antarctica, SMB is increasing due to increasing precipitation with present warming (Frieler et al., 2015). However, we propose that once the warming passes a tipping point, the melt process will be dominant and results in considerable surface mass loss.
The total ice sheet mass balance through time is much more complex due to calving and basal melting under grounded ice and ice shelves. However, we argue that SMB is the main process accounting for ice gain, while the other processes account for ice loss and are mainly dependent on the ice sheet geometry (Cuffey & Paterson, 2010). In the sensitivity experiments (Figure 2), the maximum ice sheet volume decreases when large temperature fluctuations are added. The larger the amplitude of the variability, the larger the decrease is. Another notable result from that is the hysteresis of the simulated ice sheet evolution, where the buildup processes are much slower than the retreat processes. One possible explanation is delayed isostatic rebound in which the fast retreat is governed mainly by rapid ablation due to the lowered surface elevation (Abe-Ouchi et al., 2013). We propose that the nonlinear relationship between SMB and temperature could be another explanation. Due to that, the response time of the ice sheet buildup and retreat processes differs substantially (Ruddiman, 2001). Finally, although the maximum states are a function of the magnitude of the variability, all simulations end with similar SLE, showing the importance of climate mean state to interglacial states. We also compared the simulated mean SMB over the ice sheet domain for experiments GL_hf and GL_lf (Figure 4b). We focus on the last two glacial cycles when the DO signals are more pronounced (Figure 4a). For GL_lf, the SMB is relatively constant (around 0.025 m/year) during the ice sheet buildup stages (e.g., 120,000-20,000 and 240,000-140,000 years BP), while negative SMB occurs mainly during glacial terminations. For GL_hf, the SMB shows a strong asymmetric response to DO fluctuations. The positive SMB is no more than 0.064 m/year in response to DO fluctuations, while negative SMB can be as low as −0.205 m/year. The relationship between the 18 O index and SMB is shown in Figures 4c and 4d. Similar to Figure 3a, the simulated SMB has a nonlinear response when climate forcing shifts between cold and warm states. The relation between the difference of 18 O (| 18 O|) and the resultant SMB (| SMB|) for the two experiments (GL_hf-GL_lf, Figure 4e) shows similar features as in Figure 3c, where the positive SMB response to negative temperature anomalies has an upper limit (blue), while the negative SMB is more scattered (red). This can lead to a significant reduction of the mean SMB than without the variability. Mikkelsen et al. (2018) demonstrated that interannual temperature fluctuations cause a smaller steady-state ice sheet volume of 1-m SLE using a minimal ice sheet model and from a theoretical derivation. In our study, we show that it is also valid for transient ice sheets under millennial-scale abrupt climate fluctuations. The maximum deviation of simulated ice sheet volume response to DO fluctuations can be up to 60-m SLE. More generally, the mean state of ice sheets is influenced by the amplitude of climate noise. Forced with high-amplitude climate fluctuations, the ice sheet volume can be reduced significantly compared to the case without climate noise (Figure 2). This phenomenon is attributed to the nonlinearity of SMB to temperature. The surface mass budgets are calculated with the PDD method. Since it is a semiempirical method and can be tuned to match the observations, we hypothesize that using a surface energy balance model would lead to similar results. Siddall et al. (2008) synthesized various reconstructed sea level records and found millennial-scale sea level fluctuations of 20-to 30-m magnitude during Marine Isotope Stage 3 (MIS 3, 60-27 kyr BP). Sea level fluctuated between 40 and 80 m during MIS 4 (70-60 kyr BP) and MIS 3 (Chappell, 2002;Lambeck & Chappell, 2001), which is in line with the GL_hf simulation (Figure 1b, black line). Comparing our results to the reconstructed sea level record from (Rohling et al., 2014; blue line), we found a better match when the DO fluctuations are considered (Figure 1b). The terrestrial ice sheet reconstructions before the LGM are not well constrained, especially for North American ice sheets, due to erosion during the most recent glaciation son Bay Lowlands suggests that the Laurentide Ice Sheet was significantly reduced during MIS 3 and grew rapidly toward the LGM (Dalton et al., 2019). Terrestrial chronologies from southern Greenland during early to late Holocene show that the Greenland ice sheet margins are sensitive to not only long-term (larger than 1,000 years) but also short-term (less than 100 years) climate fluctuations (Reusche et al., 2018). Based on this evidence, we propose that the Northern Hemisphere ice sheets possibly underwent larger variability in response to the millennial-scale abrupt climate shifts than previously assumed (Stokes et al., 2012).
This study presents a novel perspective on the interplay between ice sheet evolution and climate noise. During glacial inception, smaller climate variability created conditions favorable for ice sheet buildup. As ice sheets grew, larger climate variability arose and prevented the ice sheets from further growth. This forms a negative feedback between the climate variability and ice sheet evolution, leading the Earth system toward an intermediate state. The simulated ice sheets without climate variability reached full glacial conditions much earlier than when there was significant climate variability (Figure 1b), which can be more than 40,000 years. During MIS 4, the boreal summer insolation decreased to as low as the LGM, while CO 2 value was relatively high (∼205 ppm) and the ice sheets were at an intermediate size (Ahn & Brook, 2008;Lisiecki & Raymo, 2005). Van Meerbeeck et al. (2009) found a warmer equilibrium climate at MIS 3 than at the LGM and showed that the climate sensitivity to insolation forcing is higher than greenhouse gas forcing. Schaefer et al. (2015) showed evidence that the Southern Hemisphere went from an interglacial condition to a full glacial condition at ∼65,100 years BP within a 15,000-year interval. They proposed that Northern Hemisphere ice sheets require half a glacial cycle to reach full glacial conditions. In our simulation (GL_hf), Northern Hemisphere ice sheets further build up with lower insolation during MIS 4 and reach an intermediate state with 80 m SLE. However, without the DO oscillations (GL_lf), ice sheets reach a full glacial state at 64 kyr BP and maintain their shape through MIS 3. In other words, the climate variations favored maintaining an intermediate ice sheet configuration and corresponding climate. We propose that the large variability slowed down the ice sheet buildup process and prevented the Earth system from entering a full glacial state from MIS 4 to MIS 3.
In our study, we emphasize the important role that DO fluctuations played on Northern Hemisphere ice sheet evolution. Similar to Cessi (1994) and Timmermann and Lohmann (2000), we show that the ice sheet mean state is regulated by DO-like fluctuations and is affected by the climatic noise level. Although the glacial index that we used may not accurately represent the millennial-scale climate variability, the nonlinear dependence of the SMB on temperature is considered to be a robust phenomenon (Niu et al., 2019). The nonlinear interactions between the ice sheet and atmosphere or ocean are not considered within this method (Abe-Ouchi et al., 2007;Liakka et al., 2012;Löfverström et al., 2015) and may affect the influence of the climate noise level on ice sheet development. A more complex coupled model is needed to account for the interplay between the ice sheet, the ocean, and the atmosphere or other climate components and to quantify the effect of climate noise on the climate (Rehfeld et al., 2018).

Conclusions
Our study shows that climate noise can significantly slow ice sheet development and influence the ice sheet mean state. The larger the amplitude of the variability, the smaller the resultant ice sheets. We interpret that the nonlinear dependence of the SMB on the temperature is the main reason for this phenomenon. In addition, this nonlinearity could also result in a hysteresis of the ice sheet evolution, whereby the buildup process is much slower than the retreat process. We propose that the DO events may significantly inhibit the development of Northern Hemisphere ice sheets and prevented the Earth system from reaching a full glacial state from MIS 4 to MIS 3. This gives us a new perspective on understanding the role of climate variations played within the Earth system and shows that the Northern Hemisphere ice sheets are sensitive to millennial-scale abrupt climate shifts.