A modified seasonal cycle during MIS31 superinterglacial 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 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 signicant concentration of power in the 3-7 year band and on the multidecadal time scale between 15-30 years. However, the MIS31 simulation shows drastically 5 modified temporal variability of the ENSO, with stronger power spectrum at interannual time scales but absence of the decadal periodicity. Increased meridional gradient of 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 distinct global monsoon system and its link to the ENSO in respect to current conditions. In particular, the 10 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 significatly correlated with the NINO34 index. This indicates that monsoonal precipitation for this interglacial is more closely connected to hemispherical features than to the tropical-extratropical climate interaction.

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, 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 paleore-5 constructions 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 multi-decadal timescales in distict epochs, such as inverstigated in present study, is vital to understand past interglacial intervals The implementation of MIS31 Antarctic topography differs from the CTR counterpart primary 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 2 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 5 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) for the 1980-2010 interval. Our CTR climate is warmer than the is slightly warmer than ERA-I and 10 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).
Visiting the global distribution of surface temperature ( Fig. 1 Supp material) 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 of 15 temperatures as shown by the standard deviation (STD). It is demonstrated that higher values of STD 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 limitation is noted over steep topographies such as Tibet plateau, Andes and Rocky mountain.
Temperature differences between the MIS31 and the CTR show that most of warming occurs in the boreal summer, reaching 20 1.2 • C in 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 2 (Supp. Material) shows the monthly averaged hemispheric pattern for surface solar radiation (SSR) and surface 25 temperatures delivered by MIS31 and CTR simulations. This figure demonstrates an inter-hemispheric seesaw emphasizing the substantial increase in the boreal SSR during the summer season in the MIS31 experiment, and similar situation occurs in the Southern Hemisphere during DJF in the extra-tropics. It has to be argue 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 2 Supp. Material). Thus, much warmer summer conditions and colder winter/spring in the SH increase the annual amplitude.

5
The inclusion of distinct astronomical forcing leads to NH peak summer (June/July) insolation, with an opposite effect in the SH, due to the interterhemispheric seesaw relationship of the precession cycle (Scherer et al., 2008;Erb et al., 2015). 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 (Fig 2 Supp. Material).
The inclusion of MIS31 boundary conditions also results in changes in sea level pressure (SLP) and the vertical structure 10 of the atmosphere. Figures 1a,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).
Moreover, this strategy is important, because changes in circulation are dictated by changes to the gradient of geopotential 15 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 in 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 20 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 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. 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). As discussed by Yang and Wang (2009) for the equatorial Pacific, changes in the depth of the thermocline determines the SST magnitude and the behavior of the air-sea interaction, influencing the phase, amplitude, and time scale 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 (glo, 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 5 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 10 Ocean Drilling Program 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 (Figs. 1c,e) lead to cooler SSTs and enhanced subtropical cell, in concert with an intensified subtropical gyre (Figs. 1g,h). The wind-driven circulation may be evaluated by the Sverdrup transport defined as: 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; 20 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 (ie, 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 sub-surface water related to the speed up of the subtropical cell (Fig. 1h). These 25 processes are in phase with increased precipitation in the central Pacific region, but dryer conditions are noted in the 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).

Enhanced seasonality in MIS31
The use of harmonic analysis allows the identification of dominant climate signals in the space-time domain, separating small 30 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 time-scales. Thus, different climate regimes and transition regions can be characterized. The 1st 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(Justino et al., , 2016. Thus, further evaluation on modifications of the annual and semi-annual cycle in the MIS31 and CTR simulations are 5 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 1 st 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 which strongly dominates SST, windstress, andprecipitation (Li andPhilander, 1996, 1997). Nevertheless, the western equatorial Pacific and to a lesser extent the western Atlantic temporal variability present largest variance in the semi-annual component (2 nd harmonic). The  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 25 monsoon forcing, windstress 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 intesification of the semi-annual oscillation between 10-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 (Warm Pool region). In the Atlantic, the 1 st harmonic weakens, 30 allowing for sub-seasonal temporal variability (lower order harmonics) enhanced nearby the African coast (Fig. 2c).

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, shallow thermocline as delivered by the MIS31 simulation indicates reduced upper-ocean heat content that may intensify the high-frequency in the equatorial region, in particular the interannual ENSO variability (An et al., 2004;An and Jin, 2000), as further discussed.

5
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 case of 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 10 NIÑO4 (Wilson et al., 2014.
Figures 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 timeseries, 3 tapers have been used to resolve spectral fluctuations at frequencies greater than the Rayleigh frequency (MTM; Thomson, 1982).
All periodicities mentioned below are significant at the 95% confidence level. Compared with the power spectrum deliv-15 ered by the HadISST, the ICTP-CGCM shows sharper peak 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 time scale between 15-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 sub-tropical 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 20 of NIÑO4 is shifted to higher frequencies compared to the other two indices.
The reason to this slightly shift to higher frequency by the NINO4 is not clear, however, because the NIÑO4 is located much closer to the warming pool region, which is dominated by weak seasonal cycle with the 1st harmonic explaining by about 30% of the total variance, 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 25 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 windstress within 10-30 • N is considerably weaker than that of NIÑO34 or NIÑO3. Moreover, this weaker correlation between the NIÑO4 and windstress 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.   30 3c). This simulation shows stronger power spectrum at interannual time scales between 3-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 straitforward. 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 windstress 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. 3 Supp. material) that shows enhanced power also at interdecadal time scales (Supp. Material Fig. 3).
The opposite is delivered by the MIS31 simulation, a fact that usefully serves to support the assumption of weaker decadal previous studies have claimed that the equatorial Pacific interannual variability is primary forced by equatorial windstress (Nonaka et al., 2002;Timmermann and Jin, 2002), and the decadal variability is strongly connected to the off-equatorial 10 windstress, our results show that the atmospheric flow between 0-20 • N can induce decadal variability (Supp. Material Fig. 3).
In fact, the decadal variability found in the CTR NIÑO34 power spectrum fits nicely with the proposed mechanism raised Turning to the regression patterns induced by the NIÑO34 indices, Figure 4 shows that our coupled model reproduces the 20 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).
The impact of NIÑO34 on SLP (Fig. 4c) (Figs. 4e,f). In the following, we compare temperature differences between the MIS31 and CTR to compiled data by Wet et al. (2016), but with 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 35 that the ENSO works in line with the astronomical forcing, inducing warming ( √ in Table 2) however in some cases it acts in the opposite (X in Table 2).

Global and monsoonal precipitation
As shown in Figure 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).

5
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 it is evaluated the correlation between precipitation computed over regional Asia monsoonal domains, as defined by Yim et al. (2014) (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 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/winter precipitation.
Turning to individual monsoonal domains, it is demonstrated that during the MIS31 the link between the NINO34 and 20 the Asia and Australia monsoon 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 monsoon. Under current conditions, the NINO34 is negatively correlated with monsoonal indices with values by up to -0.5 for the AM and AUS, in which the NINO34 leads by 2 months (Fig. 5a,e). However, despite stronger interannual ENSO in the MIS31 climate the correlation between the SST and 25 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 NINO34 and the monsoon, Figure 5 also show the correlation between the indices but filtering out the interannual periodicity.    Wet et al. (2016)) and NOAA-SST/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 Figure 1. a) Annual differences between MIS31 and CTR runs for SLPe and b) for Z200 (mb). c) shows differences of monthly SST ( • C) and windstress (vector, Nm −2 ) between MIS31 and CTR runs. d) differences of thermocline depth between MIS31 and CTR simulation (meter).
e) and f) show differences between MIS31 and CTR zonal windstress (Nm −2 ), and precipitation differences between MIS31 and CTR runs (mm/day). g) is the Sverdrup transport in the CTR and h) the difference between the Sverdrup transport in the MIS31 and CTR runs.