A modified seasonal cycle during MIS31 super-interglacial favors stronger interannual ENSO and monsoon variability

It has long been recognized that the amplitude of the seasonal cycle can substantially modify climate features in distinct timescales. This study evaluates the impact of the enhanced seasonality characteristic of the Marine Isotope Stage 31 (MIS31) on the El Niño–Southern Oscillation (ENSO). Based upon coupled climate simulations driven by present-day (CTR) and MIS31 boundary conditions, we demonstrate that the CTR simulation shows a significant concentration of power in the 3–7-year band and on the multidecadal timescale between 15 and 30 years. However, the MIS31 simulation shows drastically modified temporal variability of the ENSO, with stronger power spectrum at interannual timescales but the absence of decadal periodicity. Increased meridional gradient of sea surface temperature (SST) and wind stress in the Northern Hemisphere subtropics are revealed to be the primary candidates responsible for changes in the equatorial variability. The oceanic response to the MIS31 ENSO extends to the extratropics, and fits nicely with SST anomalies delivered by paleoreconstructions. The implementation of the MIS31 conditions results in a distinct global monsoon system and its link to the ENSO in respect to current conditions. In particular, the Indian monsoon intensified but no correlation with ENSO is found in the MIS31 climate, diverging from conditions delivered by our current climate in which this monsoon is significantly correlated with the NIÑO34 index. This indicates that monsoonal precipitation for this interglacial is more closely connected to hemispherical features than to the tropical–extratropical climate interaction.

1 Introduction

The Marine Isotope Stage 31 (MIS31; early Pleistocene 1085–1055 ka) is a key paleoclimate period to simulate and analyze the global environmental response to a significantly modified climate forcing (Lisiecki and Raymo, 2005; Yin and Berger, 2012). This interval was characterized by boreal summer temperatures that were several degrees greater than modern climate (up to 6 ^∘C), with a substantial recession of the Northern Hemisphere (NH) sea ice (Melles et al., 2012; Justino et al., 2017). On long timescales, Earth's climate is primarily controlled by external and internal processes related to the astronomical forcing and the atmospheric concentration of greenhouse gases (IPCC, 2013; Erb et al., 2015). Internal modes of climate variability, such as the El Niño–Southern Oscillation (ENSO), the Pacific Decadal Oscillation (PDO), and the Northern Annular Mode (NAM), also induce climate anomalies on interannual and decadal timescales (Bjerknes, 1964; Mantua et al., 1997; Thompson and Wallace, 2001).

Indeed, changes in the physical and dynamical characteristics of the ENSO have been related to seasonal and interannual global-climate disturbances (Cai et al., 2014). Moreover, changes in the amplitude of the seasonal cycle can substantially modify climate features on distinct timescales. The impact of equatorial dynamics and ENSO have been found in different equilibrium climates forced by glacial and interglacial conditions (Karamperidou et al., 2015). For instance, paleoreconstructions have demonstrated a significant reduction in the climate variability associated with ENSO during the mid-Holocene (≈6000 BP (years before present; Karamperidou et al., 2015). The glacial maximum climate was also affected by distinct ENSO variability as well; however, during this time the ENSO demonstrated larger-amplitude self-sustained interannual variations compared to current conditions (Tudhope et al., 2001; An et al., 2004; Toniazzo, 2006; Zhu et al., 2017).

The effect of ocean dynamics also modifies the tropical–extratropical interaction due to different ENSO flavors (Wilson et al., 2014, 2016), which results in anomalous atmospheric and oceanic circulations (Steig et al., 2013; Wilson et al., 2014, 2016). However, it can be argued that there is no preferentially dominant region in the equatorial Pacific, because global climate disturbances have been found in response to NIÑO3, NIÑO4, or NIÑO34 anomalies. For instance, Yin et al. (2014) indicate that warmer conditions during the MIS13, an interglacial that occurred at approximately 0.5 million years ago in the Indian Pacific warming pool, amplify the insolation effect and contribute to a large increase in summer precipitation in southern China, whereas dryer conditions occur in northern China.

The far-reaching effect of equatorial dynamics on climate has also been found by Karami et al. (2015). They argued that lower summer sea surface temperatures (SSTs) in the central tropical Pacific during MIS13 contribute to the strengthening of the northern Pacific subtropical high, increasing the transport of moisture into the East Asia summer monsoon (EASM). Moreover, they highlight the significant influence of the east–west SST differences in the tropical Pacific in maintaining the link between the tropical Pacific and EASM. Larger differences in the east–west SST gradient in the equatorial Pacific that began 1.17 million years ago has also been claimed to support the onset and intensification of the modern Walker circulation (McClymont and Rosell-Melé, 2005).

Based on paleoreconstruction of wind and precipitation on the Chinese Loess Plateau, Sun et al. (2010) demonstrated that monsoonal fluctuations at orbital-to-millennial scales is dynamically linked to changes in solar insolation, and internal boundary conditions. Therefore, it can be assumed that changes in insolation or increased temperature, as occurred during interglacial stages, may trigger a distinct pattern of global monsoon, and likewise can be expected in the future (Hsu et al., 2012).

Significantly modified periodicity and amplitude of past ENSO regimes, and their global influence, shed light on the potential effect of human-induced climate change on the equatorial Pacific, and consequently on future ENSO-like climate. Furthermore, it should be argued that disagreement in the magnitude of cooling or warming among coupled climate models and paleoreconstructions may be related to the local responses of temperature and precipitation elicited by distinct ENSO in both spatial and temporal variability (Peltier and Solheim, 2004; Jost et al., 2005; Yin and Berger, 2012; Dolan et al., 2015; Justino et al., 2017).

Thus, understanding of the air–sea interaction related to the equatorial Pacific and its climate response at interannual and multidecadal timescales in distinct epochs, such as those investigated in the present study, is vital to understand past interglacial intervals that are characterized by depleted ice sheets. Verifying the potential effect of atmospheric CO₂ on the stability of the West Antarctic Ice Sheet (WAIS) will be a key climate factor in decades to come (Nicolas et al., 2017).

2 Coupled climate simulations

2.1 Climate model and experimental design

Climate simulations have been performed with the International Centre for Theoretical Physics – Coupled Global Climate Model (ICTP-CGCM;  Kucharski et al., 2016). ICTP-CGCM consists of the atmospheric global climate model “SPEEDY” version 41 (Kucharski et al., 2006) coupled to the Nucleus for European Modelling of the Ocean v3.3 (NEMO) model (Madec, 2008) with the OASIS3 coupler (Valcke, 2013).

The atmospheric component runs at T30 horizontal resolution, and there are eight levels in the vertical. NEMO is a primitive equation z-level ocean model based on the hydrostatic and Boussinesq approximations. This version applies a horizontal resolution of 2^∘ and a tropical refinement to 0.5^∘. The ocean component has 31 vertical levels with layer thicknesses ranging from 10 m at the surface to 500 m at the ocean bottom (16 levels in the upper 200 m). Additional details on the ICTP-CGCM coupled model are discussed by Justino et al. (2017). Farneti et al. (2014) used the ICTP-CGCM to examine the interaction between the tropical and subtropical northern Pacific at decadal timescales, suggesting that extratropical atmospheric responses to tropical forcing feed back onto the ocean dynamics leading to a time-delayed response of the tropical oceans.

Two simulations are evaluated: a modern climate driven by present-day boundary conditions (CTR) and a second experiment for the MIS31 forcing. The CTR simulation was run to equilibrium for 2000 years, and our modern climate is the time average of the last 500 years of the CTR simulation. The CTR is run under present-day orbital forcing and CO₂ concentration of 325 ppm as it characterizes emission by the year 1950. The MIS31 run starts from equilibrated CTR conditions, including modifications of the WAIS topography based on Pollard and DeConto (2009), and the planetary astronomical configuration of 1.072 Ma according to Coletti et al. (2015). It has been carried out for 1000 years and the analyses take into account the last 500 years of the simulation.

The implementation of MIS31 Antarctic topography differs from the CTR counterpart primarily by the absence of the WAIS, which, according to Pollard and DeConto (2009), was induced by changes in ocean melt via the effect on ice-shelf buttressing that coincides with strong boreal summer insolation anomalies. In all experiments, the CO₂ concentration was set to 325 ppm, which is based on boron isotopes in planktonic foraminifera shells for the MIS31 interval (Honisch et al., 2009). The MIS31 and CTR experiments have been described in further detail elsewhere by Justino et al. (2017), but a brief discussion of the global climate differences between these two runs is provided below.

2.2 Mean climate conditions for MIS31 and CTR

2.2.1 Atmospheric conditions

Table 1 shows the global and hemispheric surface temperature values for the CTR and MIS31 simulations, and ERA-Interim (ERA-I;  Dee et al., 2011) and HadCRUT4 for the 1980–2010 interval. Our CTR climate is slightly warmer than ERA-I and larger differences are noticed for the NH summer when the CTR climate is 1 ^∘C warmer. These differences arise from higher temperatures over land, because the SST and sea-ice distributions in the CTR simulation fit well with the ERA-I dataset, as shown by Justino et al. (2017).

Table 1 Averaged surface temperatures (^∘C) for the CTR, ERA-I, and HadCRUT4 (1980–2010). Differences between MIS31 and the CTR runs are also shown. Values in brackets are for JJA (June, July, and August). Otherwise values are computed for DJF (December, January, and February).

When viewing the global distribution of surface temperature (Fig. S1 in the Supplement), it is demonstrated that the ICTP-CGCM is able to reproduce the main features of the ERA-I data. The ICTP-CGCM performs fairly in reproducing the monthly variability in temperatures as shown by the standard deviation (SD). Higher values of SD are over Asia and North America related to the high seasonality due to the effect of continentality. Larger values are also observed over oceanic regions along the storm track. However, due to the model resolution, a limitation is noted over steep topographies such as the Tibetan Plateau, Andes, and Rocky Mountains.

Temperature differences between the MIS31 and the CTR show that most warming occurs in the boreal summer, reaching 1.2 ^∘C for the global mean, 2.2 ^∘C in the NH, and 0.4 ^∘C in the Southern Hemisphere (SH). Lower temperatures are demonstrated during DJF in the MIS31 run compared to the CTR simulation, clearly showing the hemispheric seesaw effect of the astronomical forcing. Zonally averaged, the MIS31 climate is remarkably warmer than the CTR during JJA except poleward of 45^∘ S. During DJF, the MIS31 is slightly cooler between 45^∘ S and 50^∘ N.

Figure S2 shows the monthly averaged hemispheric pattern for surface solar radiation (SSR) and surface temperatures delivered by MIS31 and CTR simulations. This figure demonstrates an interhemispheric seesaw emphasizing the substantial increase in the boreal SSR during the summer season in the MIS31 experiment, and a similar situation occurs in the Southern Hemisphere during DJF in the extratropics. It has to be argued that the reason for larger seasonality in the SH for the CTR run is related to the excess of SSR in DJF but deficit in JJA as compared to the NH (Fig. S2). Thus, much warmer summer conditions and colder winter/spring conditions in the SH increase the annual amplitude.

Due to astronomically driven reduced sea ice, larger changes are located in the NH extratropics (see Justino et al., 2017). It has to be mentioned that differences between the MIS31 and CTR simulations (the last 500 years of both runs) deliver negative surface temperature anomalies over southern Asia, the western equatorial Pacific, and South Atlantic. Compared to tropical and extratropical paleoreconstructions, the MIS31 simulation performs well with values that depart from paleoproxies by ±1 ^∘C in the tropics between 20^∘ N and 3^∘ S (McClymont and Rosell-Melé, 2005; Medina-Elizalde et al., 2008; Herbert et al., 2010; Li et al., 2011; Russon et al., 2011; Dyez and Ravelo, 2014) by up to −3 ^∘C within the 41 to 67^∘ N belt (Raymo et al., 1996; Herbert et al., 2010; Li et al., 2011; Naafs et al., 2013) and ±1.5 ^∘C in the SH between 23 and 42^∘ S (McClymont et al., 2005; Crundwell et al., 2008; Scherer et al., 2008; Naish et al., 2009; Martínez-Garcia et al., 2010; Russon et al., 2011; Voelker et al., 2015).

The inclusion of distinct astronomical forcing leads to NH peak summer (June/July) insolation, with an opposite effect in the SH due to the interhemispheric seesaw relationship of the precession cycle (Scherer et al., 2008; Erb et al., 2015).

The inclusion of MIS31 boundary conditions also results in changes in sea level pressure (SLP) and the vertical structure of the atmosphere. Figure 1a and b show the eddy SLP (SLP with the zonal mean removed, SLP_e) and the Z200 (geopotential height at 200 hPa with the zonal mean removed). We have shown the SLP_e because differences between high- and low-pressure-dominant systems in the MIS31 and CTR, such as at the subtropical North Pacific and Azores High are enhanced. This facilitates the interpretation of wind anomalies at the subtropics and equatorial region (e.g., the trade wind anomalies).

Figure 1 (a) Annual differences between MIS31 and CTR runs for SLP_e and (b) for Z200 (mbar). Panel (c) shows differences in monthly SST (^∘C) and wind stress (vector, N m⁻²) between MIS31 and CTR runs. (d) Differences in thermocline depth between MIS31 and CTR simulation (meter). Panels (e) and (f) show differences between MIS31 and CTR zonal wind stress (N m⁻²), and precipitation differences between MIS31 and CTR runs (mm d⁻¹). Panel (g) is the Sverdrup transport in the CTR and (h) the difference between the Sverdrup transport in the MIS31 and CTR runs.

Moreover, this strategy is important because changes in circulation are dictated by changes to the gradient of geopotential rather than absolute magnitude anomalies (He et al., 2018). At the surface, SLP_e anomalies exhibit an increase in the western North Pacific subtropical high, opposite to the drop in the eastern North Pacific (Fig. 1a). This partially supports previous results by Mantsis et al. (2013), who found a large strengthening and a northward and westward expansion of the northern Pacific summer anticyclone, driven by changes in the timing of perihelion.

According to Cook and Held (1988) and Timmermann et al. (2004) the meridional circulation v^′ is proportional to the mean westerly circulation u>0, which is also modulated by the seasonal cycle of the SLP. In the upper troposphere (200 hPa), this induces southward anomalies over the eastern Pacific and northward, and low-pressure anomalies on the downstream side of the Tibetan Plateau in MIS31 (Fig. 1b); hence, weakening the jet stream in the MIS31 climate compared to the CTR counterpart. This vertical structure with a baroclinic anomalous pattern in particular over East Asia and western Pacific may be related to the ENSO dynamics in the MIS31 climate, as will be verified later.

2.2.2 Oceanic conditions

These changes in the stationary wave induce substantial modifications in the wind stress and SST near-surface air temperature features delivered by the MIS31 climate. Indeed, Fig. 1c depicts warmer SSTs in the northeast Pacific and subtropical Pacific but cooler temperatures in the west Pacific. These changes along the equatorial belt are primarily induced by weaker northeast trade winds that reduced evaporative cooling and lead to less vigorous equatorial upwelling between 0 and 20^∘ N. Moreover, wind stress changes in the eastern Pacific reduce the cold tongue strength (Fig. 1c, d).

Modification in the near-surface atmospheric circulation can also modify the oceanic vertical characteristics affecting the thermocline depth and ENSO (Wen et al., 2014; Bush, 2013). As discussed by Yang and Wang (2009) for the equatorial Pacific, changes in the depth of the thermocline determine the SST magnitude and the behavior of the air–sea interaction, influencing the phase, amplitude, and timescale of the tropical climate.

The ICTP-CGCM properly reproduces the equatorial thermocline depth (using the depth of maximum vertical temperature gradient) compared to the Levitus dataset (Levitus et al., 2000) and Glorys reanalysis (Parent et al., 2013); however, our CTR climate shows a much shallow thermocline off the equatorial region in the SH. The MIS31 forcing leads to a shallower thermocline and reduction of its zonal gradient (Fig. 1d), which is primarily related to the anomalous wind flow (e.g.,  Zebiak and Cane, 1987; An et al., 1999).

A deeper thermocline, however, is observed in part of the NIÑO3 region (Fig. 1d, contour). In the eastern Pacific, thermocline dynamics have been associated with changes in SST, the air–sea coupling, and ENSO (Leduc et al., 2009; Yang and Wang, 2009). This implies a weaker Walker circulation during the MIS31 interval that is supported by SST reconstructions (from Ocean Drilling Program, ODP, sites 849, 847, 846, and 871) in the western and eastern equatorial Pacific (McClymont and Rosell-Melé, 2005).

Over the western Pacific, stronger equatorward winds (Fig. 1c, e) lead to cooler SSTs and an enhanced subtropical cell, in concert with an intensified subtropical gyre (Fig. 1g, h). The wind-driven circulation may be evaluated by the Sverdrup transport defined as

$$\begin{array}{}\text{(1)}& \mathit{\psi }\left(x\right)={\displaystyle \frac{\mathrm{1}}{\mathit{\beta }\mathit{\rho }}}\underset{x_e}{\overset{x}{\int }}{\displaystyle \frac{\partial {\mathit{\tau }}_x}{\partial y}}\mathrm{d}x,\end{array}$$

where β is the meridional derivative of the Coriolis parameter, ρ is the mean density of sea water, and τ_x is the zonal component of the wind stress. The integral is computed from the eastern to the western boundary in the North Pacific using modeled atmospheric wind stress data. The ICTP-CGCM model simulates the Sverdrup transport quite well (Fig. 1g) compared to the magnitude of the Sverdrup transport estimated from the International Comprehensive Ocean-Atmosphere Data Set (ICOADS;  Woodruff et al., 2011). It has to be mentioned that even though the wind grid-resolution is coarser than the horizontal scale of the western boundary current (i.e., the Kuroshio Current), the τ_x used in the calculation is a representation of the zonal-averaged wind stress so that it is still fine for this analysis.

The intensification of the Sverdrup transport by up to 6 Sv between 20–40^∘ N in the Kuroshio region induces negative SST anomalies due to the intrusion of colder subsurface water related to the speed up of the subtropical cell (Fig. 1h). These processes are in phase with increased precipitation in the central Pacific region, but dryer conditions are noted in the Tropical Warm Pool region (Fig. 1f). The convergence of wind anomalies (Fig. 1c) also indicates westerly flow in concert with a shallow thermocline as delivered by the MIS31 simulation, with potential implications for the ENSO dynamics (Eisenman et al., 2005).

3 Enhanced seasonality in MIS31

The use of harmonic analysis allows for the identification of dominant climate signals in the space–time domain, separating small- and high-frequency processes (e.g., diurnal cycle) from large-scale features (e.g., seasonal). Analyses conducted on the frequency domain can capture and differentiate the contribution of all timescales. Thus, different climate regimes and transition regions can be characterized. The first harmonic shows the dominance of the annual cycle when most of the variance is represented by this harmonic. It has to be stressed that investigations based upon area-averaged time series are embedded with small- and large-scale processes dictated by distinct periodicity; this in turn hampers the identification of periodic climatic signals in the space–time domain (Justino et al., 2010, 2016).

Thus, further evaluation on modifications of the annual and semiannual cycle in the MIS31 and CTR simulations are provided below. Changes in the harmonic variance and amplitude are highly correlated with the amount of incoming shortwave radiation (SSR) in the MIS31 climate, as shown by differences in the first harmonic (Fig. 2a–f) of SSR, SST, SLP, and the oceanic heat flux (HF). It has long been recognized that the equatorial climate exhibits an annual component that strongly dominates SST, wind stress, and precipitation (Li and Philander, 1996, 1997). Nevertheless, the western equatorial Pacific and to a lesser extent the western Atlantic temporal variability present largest variance in the semiannual component (second harmonic). The semiannual component is strongly influenced by ocean–atmosphere interactions in which the surface atmospheric flow and SST feed back on the cloud structure, further modifying the SSR and oceanic heat flux.

Figure 2 Differences in the first harmonic variance between the MIS31 and CTR run for (a) surface solar radiation (SSR; W m⁻²), (b) SST (^∘C), (c) SLP (mbar), and (d) surface net heat flux (W m⁻²). Panels (e), (f), (g) and (h) are the same as (a), (b), (c) and (d) but for the first harmonic amplitude.

Figure 2a–d reveal that the semiannual component, which is dominant in the western equatorial Pacific under current conditions, is weaker in the MIS31 climate allowing larger variance in the annual harmonic. This is highlighted in particular by SST and SLP distributions that potentially impact on ENSO characteristics (Fig. 2b–c). It is also shown that similar patterns are displayed by the SSR and HF, and SST and SLP harmonics. The former (SSR and HF) experiences an interhemispheric distribution, whereas the latter (SST and SLP) is dominated by an equatorial east–west dipole in the Pacific.

Figure 2e–h show changes in the amplitude of the first harmonic between the MIS31 and CTR simulations. The main features are shown as larger (smaller) NH (SH) amplitudes delivered by the MIS31 run, in particular along the continental margins and in the Tropical Warm Pool or western Pacific area (Fig. 2f). Insofar as the western Pacific changes are concerned, it has been found that the local increase in wind stress during JJA driven by the seasonality of SLP over the central and western Pacific are in concert with the higher SSR, SST, and HF amplitudes. These changes in seasonality dramatically alter the MIS31 climate compared to the CTR climate in both spatial patterns and the main mode of variability (further discussed below).

This structure is not seen in the equatorial Atlantic where variance differences between the MIS31 and the CTR are meridional. In fact, under CTR conditions this can be interpreted as the tropical Atlantic variability (TAV) related to the continental monsoon forcing, wind stress, and air–sea interaction (Deser et al., 2010). However, due to orbitally driven changes in SSR (Fig. 2a), the MIS31 climate in the tropical Atlantic shows weakening of the annual component southward to 10^∘ N, and an intensification of the semiannual oscillation between 10 and 20^∘ N compared to the CTR run (Fig. 2b).

The SLP differences are more complex, showing a pattern that differs from zonal or meridional features (Fig. 2c), even though they are correlated with SSTs in the western Pacific (Tropical Warm Pool region). In the Atlantic, the first harmonic weakens, allowing for subseasonal temporal variability (lower-order harmonics) enhanced near the African coast (Fig. 2c).

4 MIS31 – Temporal and spatial characteristics of ENSO

It is expected that those changes in the atmospheric zonal and meridional circulations and the wind-driven oceanic flow can result in modifying ENSO frequency and power. Moreover, a shallow thermocline as delivered by the MIS31 simulation indicates reduced upper-ocean heat content that may intensify the high-frequency cycle in the equatorial region, in particular the interannual ENSO variability (An et al., 2004; An and Jin, 2000), as further discussed.

Figure 3 Multi-taper power spectrum of NIÑO3, NIÑO34, and NIÑO4. Panels (a), (c), and (e) show the power spectrum for the CTR simulation. Panels (b), (d), and (f) are the same as (a), (c), and (e) but for the MIS31 simulation. Red, green, and blue lines show the 99 %, 95 %, and 90 % significance levels, respectively.

The following explores the influence of the MIS31 forcing on ENSO indices. Among several mechanisms related to ENSO dynamics, the magnitude of the seasonal cycle in the equatorial region characterizes its onset, intensity, and frequency (Liu, 2002a; Nonaka et al., 2002; Timmermann et al., 2007). It has been argued that in the case of a strong seasonal cycle, the ENSO signal can be locked in phase and frequency with this external forcing, thus reducing its magnitude. The ENSO signal may also differ in strength and influence if computed over distinct oceanic regions, such as those defined by NIÑO3, NIÑO34, or NIÑO4 (Wilson et al., 2014, 2016).

Figure 3 shows the ENSO power spectrum computed for the NIÑO34, NIÑO3, and the NIÑO4 using the CTR and MIS31 simulations dataset. This is achieved by applying the multi-taper method to detrended time series, three tapers have been used to resolve spectral fluctuations at frequencies greater than the Rayleigh frequency (multi-taper method, MTM;  Thomson, 1982).

All periodicities mentioned below are significant at the 95 % confidence level. Compared with the power spectrum delivered by the HadISST, the ICTP-CGCM shows sharper peaks in the 3–7-year band for all NIÑO indices (Fig. 3a). Under CTR conditions, significant concentration of power is also dominant on the multidecadal timescale between 15 and 30 years. Similar periodicity has been previously found by Nonaka et al. (2002). Nonaka et al. (2002) attributes the equatorial decadal variability to the influence of winds in the trade wind bands, which modifies the strength of the subtropical cell. It is interesting to note that NIÑO34, NIÑO3, and NIÑO4 differ in reproducing the decadal frequency, weakest in the NIÑO4. Moreover, the spectrum of NIÑO4 is shifted to higher frequencies compared to the other two indices.

The reason for this slight shift to higher frequency in the NIÑO4 region is not clear; however, because the NIÑO4 is located much closer to the Tropical Warm Pool region, which is dominated by a weak seasonal cycle with the first harmonic explaining about 30 % of the total variance, it may indicate that higher-order harmonics play a role to induce some power at higher frequency.

The weakening of decadal variability in the NIÑO4 region may be related to wind variability in the off-equatorial tropics as proposed by Nonaka et al. (2002). This assumption has been verified by computing the correlation pattern associated with the NIÑO indices. It turns out that the NIÑO4 relationship with the zonal wind stress within 10–30^∘ N is considerably weaker than that of NIÑO34 or NIÑO3. Moreover, this weaker correlation between the NIÑO4 and wind stress is not confined to the equatorial region but extends to the extratropics.

The incorporation of MIS31 boundary conditions drastically modifies the temporal variability of the interglacial ENSO (Fig. 3c). This simulation shows stronger power spectrum at interannual timescales between 3 and 7 years. Evaluation of the main causes related to the strengthening of the interannual variability in the MIS31 climate compared with the CTR counterpart is not straightforward. It has been found that an increased meridional gradient of SST and wind stress in the NH tropics (Fig. 2a), as simulated by the MIS31 run, may lead to stronger interannual equatorial variability in the MIS31 climate (Liu, 2002a, b; Erb et al., 2015). Likewise, the weaker seasonal cycle of the wind stress in the MIS31 simulation may lead to stronger ENSO power at 3–7 years (Chang et al., 1994). The CTR climate Southern Oscillation Index (SOI) power spectrum also shows enhanced power at similar frequencies found for the NIÑO34 and the NIÑO3 indices. This is in line with the spectrum of equatorial winds (0–20^∘ N, Fig. S3) that shows enhanced power also at interdecadal timescales (Fig. S3).

The opposite is delivered by the MIS31 simulation, a fact that usefully serves to support the assumption of weaker decadal air–sea interaction during this interglacial. Indeed, for the MIS31 simulation, correlation values between the NIÑO34 index and the first principal component (PC1) of wind stress computed at 0–20^∘ N are very low; for the CTR run, these values are 0.6 when all timescales are included and 0.4 for conditions in which frequencies below 10 years have been filtered out. Though previous studies have claimed that the equatorial Pacific interannual variability is primary forced by equatorial wind stress (Nonaka et al., 2002; Timmermann and Jin, 2002), and the decadal variability is strongly connected to the off-equatorial wind stress, our results show that the atmospheric flow between 0 and 20^∘ N can induce decadal variability (Fig. S3).

In fact, the decadal variability found in the CTR NIÑO34 power spectrum fits nicely with the proposed mechanism raised by Farneti et al. (2014). The SST anomalies at the Equator induce changes in the wind stress curl over the western Pacific, which generate SST anomalies fluctuating on decadal timescales through tropical–subtropical interactions. Individual analyses to verify the roles of the North Pacific and South Pacific in inducing the decadal variability demonstrate that most changes in power can be explained by the NH contribution. Interestingly, the NIÑO34 and wind stress anomalies between 20–0^∘ S are highly anticorrelated with values of about 0.6 in both simulations. However, only in the MIS31 climate the wind stress spectrum does exhibit enhanced interannual variability, indicating that the MIS31 ENSO dynamics are also driven at the 3–7-year period by the SH flow. A fact that is not seen for the CTR climate.

Turning to the regression patterns induced by the NIÑO34 indices, Fig. 4 shows that our coupled model reproduces the main tropical SST response to NIÑO34 (Fig. 4a), compared for instance with Cai et al. (2015). The patterns are displayed as amplitudes by regressing hemispheric anomalies on the standardized first principal component time series. The intensification of the NIÑO34 signal does not project substantial change in SST; though in the western Pacific, anomalies between ±0.3 ^∘C are noted (Fig. 4b).

Melles et al. (2012)Lawrence et al. (2009)Raymo et al. (1996)Naafs et al. (2013)Herbert et al. (2010)Herbert et al. (2010)Li et al. (2011)Dyez and Ravelo (2014)Medina-Elizalde et al. (2008)McClymont and Rosell-Melé (2005)

Table 2 Differences in reconstructed SST and Lake E temperatures (ΔT_re; Wet et al., 2016) and NOAA-SST and ERA-I. Differences between MIS31 – CTR SST and Lake E temperatures in JJA (ΔT). √ (X) stands for agreement (disagreement) between the ΔT and induced SST anomalies (MIS31-CTR) by regressing the ENSO index. NE indicates that the index was not evaluated at the grid point or anomalies are too close to zero. Based on Wet et al. (2016) and Justino et al. (2017).

Figure 4 Panel (a) is the leading empirical orthogonal function (EOF) of SST anomalies for the CTR simulation displayed as amplitudes (^∘C) by regressing hemispheric SST anomalies upon the NIÑO34 time series. Panel (b) shows SST differences between the MIS31 and CTR regressed patterns. Panels (c) and (d) are as in (a) and (b) but for SLP (mbar). Panels (e) and (f) as in (a) and (b) but for zonal wind stress (N m⁻²). Please note that figures are shown with distinct labels.

Figure 5 Lag correlations between the Asian and Australia monsoon domains and the NIÑO34 index (black for CTR and red for MIS31 simulations). (a) Asia monsoon, (b) western North Pacific monsoon, (c) East Asia monsoon, (d) Indian monsoon, and (e) Australia monsoon.

The impact of NIÑO34 on SLP (Fig. 4c) extends globally and is fairly reproduced by the ICTP-CGCM compared to the National Centers for Environmental Prediction – National Center for Atmospheric Research (NCEP-NCAR) reanalysis (Ji et al., 2015). The zonal dipole results from the contribution of the baroclinic component over the eastern Pacific and barotropic component over western Pacific, both related to the SST anomalies (Fig. 4c, a). The MIS31 NIÑO34 weakens the barotropic and baroclinic patterns of the SLP as shown by the differences in the MIS31-CTR regression (Fig. 4c, d). In the equatorial region, the anomalies are related to the intensified winds nearby the Tropical Warm Pool region but weakening midlatitude westerlies in both the northern Pacific and Atlantic (Fig. 4e, f). In the following, we compare temperature differences between the MIS31 and CTR to compiled data by Wet et al. (2016), but with a focus on the ENSO responses (Table 2).

This is achieved by comparing the modeled SST anomalies for JJA to SSTs differences between the MIS31 and CTR delivered by the regression pattern related to the NIÑO34 index (ΔT). Differences between the reconstructions and the NOAA extended reconstructed SST V3b/ERA-I are also shown (ΔT_re). Overall, model results and reconstructions agree indicating that the ENSO works in line with the astronomical forcing, inducing warming (√ in Table 2); however, in some cases it acts in the opposite way (X in Table 2).

Global and monsoonal precipitation

As shown in Fig. 1f, it is evident that most wet and dry conditions in the MIS31 compared to the CTR are in agreement with the anomalous temperature pattern and diabatic heating, in line with ENSO-related precipitation (Dai and Arkin, 2017). An exception is found over southern Asia that experiences more precipitation despite the drop in temperatures, which may indicate the contribution of extratropical large-scale atmospheric dynamics.

To further investigate the MIS31 climate features, the correlation between precipitation computed over regional Asia and Australia monsoonal domains, as defined by Yim et al. (2014) and the NIÑO34 index, is evaluated. The domains are the Asia monsoon (AM; 10–45^∘ N, 70–150^∘ E), Australia monsoon (AUS; 5–20^∘ S, 110–150^∘ E), East Asia monsoon (EA; 22.5^∘ S–45^∘ N, 110–135^∘ E), western North Pacific monsoon (WNP; 12.5–22.5^∘ N, 110–135^∘ E), and Indian monsoon (IN; 10–30^∘ N, 70–105^∘ E).

Estimates of changes in precipitation for past interglacials are still scarce, but our MIS31 simulation agrees with other studies showing enhanced Asian summer monsoon during interglacials (An, 2000; Sun et al., 2010). Enhanced seasonality and greater annual values of precipitation have also been documented across the eastern Russian Arctic in line with our modeling experiment (Melles et al., 2012). This is also true for the Asian monsoon insofar as seasonality is concerned (Fig. 5a). Controversy is raised by Oliveira et al. (2017) who argue that reduced seasonality in precipitation along the western Mediterranean region during MIS31 leads to forest decline. Reduced seasonality during the MIS31 in that region is not supported by our MIS31 climate simulation, which in fact shows an increase in both summer and fall or winter precipitation.

Turning to individual monsoonal domains, it is demonstrated that during the MIS31 the link between the NIÑO34 and the Asian and Australian monsoons is weakened with respect to the CTR characteristics (Fig. 5). It should be mentioned that the equatorial Pacific seems to have a direct influence on the AM, WNP, and the AUS monsoon precipitation under the CTR climate (Fig. 5a, b, e). This is not the case for the East Asia and Indian monsoons. Under current conditions, the NIÑO34 is negatively correlated with monsoonal indices with values by up to −0.5 for the AM and AUS, in which the NIÑO34 leads by 2 months (Fig. 5a, e). However, despite stronger interannual ENSO in the MIS31 climate the correlation between the SST and precipitation indices is extremely reduced for AM, AUS, and WNP monsoons during the interglacial period (Fig. 5e).

In order to verify the impact of decadal variability on the link between the NIÑO34 and the monsoon, Fig. 5 also show the correlation between the indices but filtering out the interannual periodicity. This reveals that the decadal variability is responsible for about 40 % of the correlation in the CTR climate. As should be expected, by removing the interannual frequency in the MIS31 climate, correlation values are very small, which indicates that changes in AUS, WNP, and AM precipitation for this interglacial are more closely connected to hemispherical features than to the tropical–extratropical climate interaction.

5 Concluding remarks

This investigation centered on a comparison between present-day conditions (CTR) and those characteristics of a super-interglacial epoch, the Marine Isotope Stage 31 (MIS31). Using coupled global climate model simulations (ICTP-CGCM), we have first demonstrated significant changes in the spatial patterns and seasonality of sea-level pressure, sea-surface temperatures, and heat fluxes during the MIS31 climate compared to present-day conditions, and these changes have a significant impact on the main modes of variability. Anomalous equatorial wind stress associated with a modified seasonal cycle in the MIS31 simulation leads to stronger ENSO variability compared to the present-day climate. Moreover, the decadal variability differs dramatically in the MIS31 simulation from that characteristic of present-day conditions. This decadal variability also differs greatly across the ENSO diversity spectrum with off-equatorial atmospheric circulation playing a significant role in inducing decadal variability.

Evaluation between paleoreconstruction and modeling results is a complex task, because reconstructions depict dominant signals in a particular time interval and locale. Thus, they cannot be assumed to geographically represent large-scale domains, and their ability to reproduce long-term environmental conditions should be considered with care.

Discrepancies between modeling results and paleoreconstructions for the MIS31 climate, which occurred under very particular conditions and high seasonality, may unfortunately be expected. The MIS31 may have been dominated for instance by vegetation patterns drastically different than today. This modifies the global evapotranspiration rates and the hydrological cycle, producing precipitation that can differ greatly from model results. This suggests that uncertainties in the model may be reduced when including more realistic boundary conditions that are currently not available.