The causes of the change in amplitude of El Niño–Southern Oscillation (ENSO) during the mid-Holocene were investigated by diagnosing the model simulations that participated in the Paleoclimate Modelling Intercomparison Project phases 2 and 3. Consistent with paleoclimate records, 20 out of the 28 models reproduced weaker-than-preindustrial ENSO amplitude during the mid-Holocene. Two representative models were then selected to explore the underlying mechanisms of air–sea feedback processes. A mixed layer heat budget diagnosis indicated that the weakened ENSO amplitude was primarily attributed to the decrease in the Bjerknes thermocline feedback, while the meridional advective feedback also played a role. During the mid-Holocene, the thermocline response to a unit anomalous zonal wind stress forcing in the equatorial Pacific weakened in both models because of the increased ENSO meridional scale. A further investigation revealed that the greater ENSO meridional width was caused by the strengthening of the Pacific subtropical cell, which was attributed to the enhanced mean trade wind that resulted from the intensified Asian and African monsoon rainfall and associated large-scale east–west circulation in response to the mid-Holocene orbital forcing.
El Niño–Southern Oscillation (ENSO), which originates in the tropical Pacific region, is the most prominent mode of climate variability at interannual time scales in the Earth climate system (Bjerknes 1969; Rasmusson and Carpenter 1982; Philander 1990). ENSO is a coupled ocean–atmosphere mode that involves both oceanic dynamical processes (Zebiak and Cane 1987; Suarez and Schopf 1988; Battisti and Hirst 1989; Jin 1997) and atmospheric responses to sea surface temperature (SST) (Gill 1982; Lindzen and Nigam 1987; Wang and Li 1993). ENSO has a great impact on weather and climate not only in the tropics but also in the extratropics. Understanding the changes in ENSO behavior is thus of broad scientific and socioeconomic importance.
How ENSO activity would change in response to global warming has been extensively studied based on observational and model data. However, it remains unclear whether future ENSO amplitude is going to weaken, strengthen, or remain unchanged (van Oldenborgh et al. 2005; Zelle et al. 2005; Guilyardi 2006; Merryfield 2006; Philip and van Oldenborgh 2006; Yeh and Kirtman 2007; Latif and Keenlyside 2009; Park et al. 2009; Collins et al. 2010; DiNezio et al. 2012; Christensen et al. 2013; Cai et al. 2014; Chen et al. 2015, 2017). Different models also project a wide range of responses from El Niño–like to La Niña–like changes in the mean state of the tropical Pacific zonal SST gradients (Collins 2005; Merryfield 2006). These uncertainties thus spur us to further understand ENSO and the link between its variability and the background climate state. Studying ENSO variability under past climate in which a mean state differs significantly from the present climate offers a unique opportunity to achieve this goal.
Paleoclimate reconstructions provide the possibility of examining the theories of ENSO response to global warming and evaluating climate models used to project this response. The climate of the mid-Holocene (MH), approximately 6000 years ago, has received much attention because a large body of evidence suggests a substantially weaker ENSO at that time (see section 7 for details). To explore the mechanism of ENSO changes during the MH, several modeling studies have been performed using coupled atmosphere–ocean simulations (Bush 1999; Clement et al. 2000; Liu et al. 2000, 2014; Kitoh and Murakami 2002; Otto-Bliesner et al. 2003; Brown et al. 2008a,b; Zheng et al. 2008; Chiang et al. 2009; Braconnot et al. 2012b; Luan et al. 2012; An and Choi 2014; Roberts et al. 2014; Emile-Geay et al. 2016). In accord with reconstructed records, these simulations have consistently indicated that ENSO activity was reduced during the MH, although the metric used to measure ENSO amplitude and the magnitude of the reduction varied greatly across the models. Excluding the study of Emile-Geay et al. (2016), who proposed that there was no clear link between orbital forcing and ENSO amplitude during the Holocene, it has been widely shown that orbitally induced change in the seasonal cycle of insolation was the root driver of the reduction in the MH ENSO variance; however, the responsible dynamic processes have remained inconclusive (Clement et al. 2000; Liu et al. 2000; Otto-Bliesner et al. 2003; Brown et al. 2008a; Zheng et al. 2008; Chiang et al. 2009).
Previous theoretical and modeling studies have proposed that the change in the background mean state is an important factor that affects the ENSO characteristics (Li 1997a; Li and Hogan 1999; Fedorov and Philander 2000; An and Jin 2001; Wang and An 2001, 2002; Chen et al. 2013, 2015; Chung and Li 2013; Xiang et al. 2013). However, there is little consensus on which of the mean state changes was essential for the changes in ENSO during the MH. Some studies attributed the reduced MH ENSO activity to the enhanced climatological mean trade winds and their changed seasonality (Clement et al. 2000; Brown et al. 2008a; Zheng et al. 2008). In contrast, some argued that changes in the mean state of the thermocline along the equator associated with the tropical Pacific zonal SST gradient were the source of the weakened ENSO variability at that time (Liu et al. 2000; Otto-Bliesner et al. 2003). More recently, Chiang et al. (2009) proposed that the reduction in the amplitude of stochastic forcing from the extratropics by teleconnection was the cause of the reduced ENSO variance during the MH. Additionally, Roberts et al. (2014) put forward the impact of the stability of ocean–atmosphere interactions, which demonstrated the importance of considering the combined effects of a mean state change on the coupled ocean–atmosphere system. Therefore, it is a challenge to uncover what part of the mean state changes was crucial and responsible for the ENSO amplitude change during the MH in state-of-the-art climate models, which calls for a thorough and quantitative diagnosis of air–sea coupling processes and feedbacks that dominate the change in ENSO amplitude.
Based on the above and motivated by the study of Chen et al. (2015), this paper presents a diagnostic analysis of the cause of the ENSO amplitude change during the MH using coupled climate models under the protocol of the Paleoclimate Modelling Intercomparison Project phase 2 (PMIP2) and phase 3 (PMIP3). The remainder of this paper is organized as follows. In section 2, we describe the models, data, and analysis methods, and evaluate the model performance for tropical climate and ENSO characteristics. Then, in section 3, two representative models are selected for the analysis of the mixed layer heat budget to diagnose the air–sea coupling processes and feedbacks during the ENSO developing phase. The causes of the changes in the two leading budget terms, the Bjerknes thermocline feedback and meridional advective feedback, are investigated in sections 4 and 5, respectively. The mechanism responsible for the thermocline response including the role of the mean state is highlighted in section 6. A comparison between model results and proxy data is performed in section 7, and conclusions and a short discussion are finally given in section 8.
2. Data and methods
a. Models and data
Based on data availability, the present study used 28 models under the framework of PMIP2 and PMIP3 for the MH climate simulations, including 17 PMIP2 atmosphere–ocean or atmosphere–ocean–vegetation coupled models and 11 PMIP3 state-of-the-art coupled models with overall higher spatial resolutions than the previous generation. Basic information about the 28 models is provided in Table 1. The boundary conditions for the MH experiment include changes in Earth’s orbital parameters (Berger 1978) and atmospheric concentrations of greenhouse gases with respect to the preindustrial (PI) period. In PMIP2 and PMIP3, the orbital parameters of eccentricity, obliquity, and angular precession during the MH (PI) were 0.018682 (0.016724), 24.105° (23.446°), and 0.87° (102.04°), respectively. The atmospheric methane concentration was reduced from the PI level of 760–650 ppb for the MH. The atmospheric concentration of carbon dioxide and nitrous oxide was held constant at 280 ppm and 270 ppb, respectively. SSTs were computed by oceanic general circulation models. The vegetation was simulated with equilibrium or dynamic vegetation models in atmosphere–ocean–vegetation models but was fixed at the present state in atmosphere–ocean models. More details about the models and experiments are provided by Braconnot et al. (2007) and Taylor et al. (2012) (and available online at http://pmip3.lsce.ipsl.fr/). Monthly model output for the last 100 yr (all years) from both PI and MH experiments were analyzed in this study using simulations with the length of run longer (shorter) than 100 yr (last column in Table 1). Note that all monthly anomalies were first obtained by subtracting the annual climatological cycle from the original time series and were then filtered using a 6–96-month Lanczos bandpass filter (Duchon 1979) with 121 weights, which eliminates the influence of any variability occurring on time scales of 8 yr or longer.
To assess model performance for the PI climate, we used monthly SSTs from the HadISST1.1 dataset (Rayner et al. 2003) and the NOAA Extended Reconstructed Sea Surface Temperature (ERSST), version 4, dataset (Huang et al. 2015) for the period 1900–99 (hereafter referred to as “observation”). To detect the ENSO amplitude and its change at a fine horizontal resolution, all model and observation data were aggregated to a grid resolution of 1° latitude × 1° longitude using bilinear or area-weighted interpolation.
b. Analysis methods
The description of the analysis methods parallels that of Chen et al. (2015), as follows. The ENSO amplitude is defined by the standard deviation (STD) of SST anomaly (SSTA) over the Niño-3 (5°S–5°N, 90°–150°W) region. A similar result can be obtained when the Niño-3.4 (5°S–5°N, 120°–170°W) SSTA is applied (figure not shown). The criterion for selecting El Niño (La Niña) events is that the Niño-3 SSTA during ND(0)J(+1) is greater (less) than one (minus one) STD. Here, ND(0)J(+1) represents the average from November–December (ND) in the El Niño or La Niña year (0) to January (J) of the following year (+1).
A mixed layer heat budget analysis was conducted to diagnose the relative role of the dynamic and thermodynamic processes in the SSTA growth during the ENSO developing phase, which is defined from April(0) to November(0). Following Li et al. (2002) and Su et al. (2010), and the same as derived from Chen et al. (2015), the mixed layer temperature tendency equation is written as
where T represents the mixed layer temperature, V = (u, υ, w) is the three-dimensional ocean current, ∇ = (∂/∂x, ∂/∂y, ∂/∂z) denotes a three-dimensional gradient operator, the prime represents the anomaly variables, and the overbar denotes the climatological annual cycle variables; is the anomalous temperature advection vertically integrated from the ocean surface to the bottom of the mixed layer, which can be decomposed into two linear advection terms and a nonlinear advection term ; Qnet represents the sum of the net downward shortwave radiation, net downward longwave radiation, and surface latent and sensible heat fluxes at the ocean surface; ρ0 is the density of seawater; Cp is the specific heat of water; H represents a constant mixed layer depth (approximately 50 m but model dependent); and r denotes the residual term. A positive heat flux represents a heating to the ocean. The variation of the mixed layer depth in space and time was not considered here because the simplicity of a mean constant mixed layer depth shows a reasonable framework for the heat budget, as is displayed in section 3.
To understand the cause of Bjerknes thermocline feedback change, we explored the air–sea coupling processes associated with the thermocline feedback. The analysis approach here is the same as derived from Chen et al. (2015), and closely resembles that applied in the Bjerknes coupled stability index proposed by Jin et al. (2006). Following Liu et al. (2011), the growth rate σ related to the thermocline feedback term is written as
where represents the mean vertical upwelling velocity, H is the constant oceanic mixed layer depth (approximately 50 m but model dependent), denotes the zonal wind stress anomaly, T′ represents the SSTA, D′ denotes the thermocline depth anomaly represented by the 20°C isothermal depth anomaly, and denotes the subsurface ocean temperature anomaly represented by the temperature anomaly at 70 m. Equation (2) shows that in addition to the mean upwelling velocity, the ENSO growth rate relevant to the thermocline feedback is determined by three feedback coefficients: , representing the zonal wind stress response in the equatorial Pacific (EP; 5°S–5°N, 120°E–80°W) to a unit SSTA forcing in the eastern equatorial Pacific (EEP; Niño-3 region); , denoting the ocean thermocline response in the EEP to a unit forcing in the EP; and , representing the ocean subsurface temperature response to a unit thermocline depth anomaly in the EEP.
Following Merryfield (2006) and Chen et al. (2015), a meridional width index was applied to calculate the meridional width of a certain or coupled field (e.g., , SSTA, or the nonlinear product of two fields) as defined by the following formula:
where W denotes the width index, x represents a certain or coupled field, and y is the latitude. The absolute value of latitude |y| serves as a weight of the field x. Considering the air–sea coupling between the zonal wind stress and SST, we calculated the width index using the regressed × SSTA field averaged over the equator between 120°E and 80°W for both PI and MH simulations.
c. Model selection
First, we examined the ability of the 28 models from PMIP2 and PMIP3 in reproducing the climatological mean SSTs over the tropical Pacific for the PI period (see Fig. S1 in the supplemental information). It was found that most models could reasonably capture the observed large-scale distribution of the SST climatology, such as the warm pool region over the western Pacific and the cold tongue in the equatorial eastern Pacific. However, the simulated equatorial cold tongue extended far too westward into the western Pacific, leading to a westward shift and even an erosion of the warm pool as compared with observations, which are the common biases of current climate models particularly for the non-flux-adjusted (nfa) ones (Latif et al. 2001). Warmer-than-observed SSTs appeared in the eastern Pacific at approximately 10°S, which implies the classical double-ITCZ bias. Note that although the spatial resolutions in PMIP3 models are overall higher than those in PMIP2 models, the above biases generally existed in the models of both two PMIP stages. We further investigated the model performance for modern ENSO characteristics, including its amplitude, frequency, and time evolution. Most models for the PI simulations reproduced a magnitude of ENSO amplitude similar to that from observations, with the values from 16 models being within ±20% of the observations and a weaker- (stronger-) than-observed amplitude occurring in seven (five) models (Fig. 1). The spectral analysis of observations for the period 1900–99 showed two statistically significant peaks at around 3.6 and 5 yr (figure not shown). Most models simulated a broadly correct dominant frequency of approximately 2–5 yr, although they exhibit too regular events as shown by a single peak. The observed time evolution of composite Niño-3 SSTA for warm and cold episodes could be properly reproduced by the PMIP2 and PMIP3 models, which generally showed the ENSO developing phase from April(0) to November(0) and the seasonal phase lock during December(0)–February(+1) (Fig. 2).
Next, we explored the response of ENSO amplitude to the MH forcing in the 28 PMIP2 and PMIP3 models. Figure 3 shows the geographical distribution of the MH minus PI (MH − PI) changes in the STD of SSTA. Although there was a large spread in spatial pattern and magnitude among the models, most [e.g., CCSM3.0, FGOALS-1.0g, MRI-CGM2.3.4 flux adjusted (fa) with global vegetation model (gvm) (MRI-CGCM2.3.4fa-gvm), and MPI-ESM-P] displayed a decrease in ENSO amplitude, whereas others showed an increase [e.g., ECBilt–CLIO–VECODE ESM with gvm (ECBILTCLIOVECODE-gvm), University of Bristol (UBRIS)-HadCM3M2-gvm, and CSIRO Mk3.6.0] or nearly no change. Quantitatively, relative to the PI period, the MH ENSO amplitude was weakened (strengthened) by 2%–21% (2%–29%) in 20 (8) models, with 15 (6) of them being significant at the 95% confidence level (Figs. 1 and 3). Note that the sample of models that simulated a weakening of the MH ENSO amplitude is not independent. There are many models in the PMIP2 and PMIP3 archive that are closely related to each other and share common components, which also show similar response to the MH forcing. For example, the two CSIRO, FOAM, and MRI-CGCM2.3.4nfa models in PMIP2 and the two HadGEM2 models in PMIP3 both simulated a weakened ENSO amplitude; CCSM3.0 and the two CSIRO models in PMIP2 and their successors, CCSM4 and CSIRO Mark 3L, version 1.2 (CSIRO Mk3L-1.2), in PMIP3, also showed a consistent weakening. However, the MH change in ENSO frequency was less consistent between models; some showed a tendency to produce either a lower- or higher-frequency peak of the Niño-3 SSTA power spectra, and others showed almost no change in the broadly dominant frequency (figure not shown). Approximately half of the models displayed a decrease in the spectral energy. The composite results for El Niño and La Niña events (Fig. 2) were generally in agreement with the change in ENSO amplitude, although only a few models showed a prominent reduction of Niño-3 SSTA during the developing phase (e.g., MIROC3.2, MRI-CGCM2.3.4fa-gvm, MRI-CGCM2.3.4nfa, and UBRIS-HadCM3M2 in PMIP2 and MPI-ESM-P in PMIP3).
Because most of the models reproduced weakened ENSO amplitude during the MH, we selected two representative models for detailed investigation, one from PMIP2 and the other from PMIP3 to shed light on the PMIP models from different stages. The criteria for model selection were as follows. First, the magnitude and spatial pattern of the tropical Pacific annual SST for the PI period were reasonably reproduced as compared with observations; second, the simulated ENSO amplitude for the PI period was within ±20% of the observations; third, the weakening of the ENSO amplitude during the MH was significant at the 95% confidence level, with an obvious negative center in the equatorial Pacific; fourth, the power spectra of Niño-3 SSTA for the period of 2–5 yr were reduced during the MH; and last and most importantly, the composite Niño-3 SSTA of El Niño or La Niña events was markedly decreased from PI to MH during the developing phase. Based on these criteria, MRI-CGCM2.3.4fa-gvm was chosen as a representative of PMIP2 models mainly because, first, the monthly ocean data are available only for this model in the PMIP2 output, and second, it showed significant and strongest weakening of ENSO amplitude by 21% during the MH with respect to PI (Fig. 3). At the same time, MPI-ESM-P was chosen from PMIP3 mainly due to its more prominent reduction of composite Niño-3 SSTA from PI to MH during the ENSO developing phase than other PMIP3 models (Fig. 2), which also showed significantly weaker-than-PI ENSO amplitude by 9% during the MH. The subsequent analysis will be based on the results of these two representative models. It is worth mentioning that CCSM3.0, CSIRO Mk3L-1.0, FGOALS-1.0g, and MIROC3.2 in PMIP2 and GISS-E2-R in PMIP3 showed ENSO behaviors similar to those of the two representative models, including a reasonable SST climatology for the PI period, a significant weakening of ENSO amplitude during the MH, and reduced composite Niño-3 SSTA from PI to MH during the ENSO developing phase (see Fig. S1 and Figs. 1–3).
3. Diagnosis of mixed layer heat budget
To investigate the cause of the weakening of the ENSO amplitude during the MH, a mixed layer heat budget analysis was performed to diagnose the mixed layer temperature anomaly (MLTA) tendency during the developing phase of the composite ENSO events for MRI-CGCM2.3.4fa-gvm and MPI-ESM-P. Figure 4 shows the composite changes in the mixed layer temperature budget terms from the two models. Each budget term during the PI and MH was first obtained for the El Niño, La Niña, and El Niño minus La Niña composites. Then, we calculated the corresponding changes from PI to MH (MH minus PI). Note that the changes in estimated (Fig. 4, bar 11) and actual MLTA tendency (Fig. 4, bar 12) were close or similar in magnitude in each composite, especially for MRI-CGCM2.3.4fa-gvm, which implies that the results of the mixed layer heat budget diagnosis are overall credible. However, the estimated MLTA budget has a relatively large error in MPI-ESM-P (Figs. 4d,f), which may arise from the mean state SST cold bias in the equatorial eastern Pacific for the PI period (Fig. S1). The changes in heat budget terms for the El Niño composite generally showed opposite signs to those for the La Niña composite. The decrease (increase) in the MLTA tendency for the El Niño (La Niña) composite means a weaker El Niño (La Niña) event. As a result, the MLTA tendency for the El Niño minus La Niña composite decreased, which was responsible for the weakening of the ENSO amplitude during the MH. More specifically, the most important terms from Eq. (1) that contributed to the weakened ENSO for both models were (term 5; Fig. 4, bar 5) and (term 8; Fig. 4, bar 8), which denote the advection of anomalous temperature by mean upwelling and by mean meridional current, representing the Bjerknes thermocline feedback and the meridional advective feedback, respectively. The contributions of terms 5 and 8 were of similar importance for MRI-CGCM2.3.4fa-gvm, while term 5 was more important than term 8 for MPI-ESM-P. For both models, the zonal advective feedback (, term 1; Fig. 4, bar 1) made a minor positive contribution to the change in the MLTA tendency due to consistent sign of their changes, whereas the advection of mean temperature by anomalous meridional current (, term 7; Fig. 4, bar 7) made a negative contribution due to their opposite changes. The thermodynamic heat flux feedback [, term 10; Fig. 4, bar 10] also made a negative contribution, which was more prominent in MRI-CGCM2.3.4fa-gvm than in MPI-ESM-P.
Note that the above two most important terms that contributed to the weakening of ENSO amplitude during the MH are the product of the mean and anomalous states. To determine whether the mean state change or the perturbation change was the major contributor, we separated their relative contributions based on the following relationship as in Chen et al. (2015):
where d(⋅) denotes the change from PI to MH (MH minus PI), A′ represents the anomalous state or perturbation part ( or ), and represents the mean state part ( or ). As deduced from Eq. (4), the change in the product of A′ and can be decomposed into three parts: the change in perturbation , the mean state change , and the change in the residual part .
Based on Eq. (4), Fig. 5 shows the composite changes in relative contributions of the mean state and the perturbation to terms 5 and 8 for MRI-CGCM2.3.4fa-gvm and MPI-ESM-P. For both leading terms in the two models, the major contributor to the MLTA tendency change for each composite was the perturbation part [ from term 5 and from term 8 in Eq. (1)]. The direct impact of the mean state change [ from term 5 and from term 8 in Eq. (1)] on the mixed layer budget difference was minor; the mean state affected the ENSO amplitude during the MH mainly through its impact on the perturbations. Therefore, we will focus the following analysis on the change in interannual variabilities.
4. Cause of Bjerknes thermocline feedback change
As mentioned above, the first leading budget term that contributed to the decrease in the MLTA tendency during the MH was the Bjerknes thermocline feedback. To understand the dynamics of the weakened ENSO amplitude, we explored the ENSO growth rate associated with the thermocline feedback term according to Eq. (2). Because the change in ENSO perturbation plays a major role, the overall thermocline feedback strength is thus determined by the product of the three air–sea coupling coefficients, represented by .
To quantitatively examine the strength of each of the three feedback coefficients, Fig. 6 shows the monthly scatter diagrams for both the PI and MH simulations for MRI-CGCM2.3.4fa-gvm and MPI-ESM-P, in which the linear slope value calculated by linear regression for each diagram denotes the air–sea coupling coefficient and is followed by its error range. The error range of the linear regression coefficient is calculated following Chen et al. (2013). Note that the wind stress–SST coupling coefficient became stronger in MH than in PI for MRI-CGCM2.3.4fa-gvm but slightly declined for MPI-ESM-P (Figs. 6a,d). Such an inconsistent change cannot explain why the thermocline feedback was weaker during the MH in both models. In contrast, the thermocline–wind stress feedback coefficient became markedly weaker in MH than in PI for both models (Figs. 6b,e), which indicates an important contribution to the consistently weakened thermocline feedback. There was little change in the subsurface temperature–thermocline feedback coefficient from PI to MH (Figs. 6c,f). Note that all the errors are small as compared to the corresponding linear regression coefficients for both models, although the error range for the thermocline–wind stress feedback coefficient is larger in MPI-ESM-P than in MRI-CGCM2.3.4fa-gvm.
The strength of each of the three feedback coefficients and their product in PI and MH are shown in Fig. 7. The magnitude in the combined effect of the three feedback coefficients represented by was greater in both PI and MH for MRI-CGCM2.3.4fa-gvm than for MPI-ESM-P, mainly due to the greater strengths in and . Although with different magnitudes, the combined strength was reduced from PI to MH for both models, confirming that the perturbation part of the thermocline feedback was responsible for the weakened ENSO amplitude during the MH. The response of D′ to (Fig. 7, bar 3) exhibited consistent change with the change in the overall thermocline feedback strength, whereas there was no consistent change in the response of to the SSTA (Fig. 7, bar 2) or in the response of to D′ (Fig. 7, bar 4) during the MH for the two models, with both two responses slightly enhanced for MRI-CGCM2.3.4fa-gvm but weakened for MPI-ESM-P. Taken together, relative to the PI period, the reduction of the thermocline feedback during the MH was primarily due to a weaker thermocline response to the zonal wind stress forcing .
Figure 8 displays the horizontal distributions of D′ regressed onto averaged in the EP, in which the regression coefficient (change) exceeding the 99% (95%) confidence level using the Student’s t test (Z test) is stippled for PI and MH (MH − PI). In response to a positive EP forcing, D′ in both PI and MH showed a consistent eastward tilting with significantly positive anomalies to the east and negative anomalies to the west for the two models (Figs. 8a,b,d,e), in accordance with the Sverdrup balance relationship between anomalous zonal wind stress and thermocline tilting (Neelin 1991; Jin 1997; Li 1997a). The difference between MH and PI (Figs. 8c,f) exhibited significantly negative D′ changes in the EEP for both models, which indicates a weaker thermocline response to the anomalous zonal wind stress and thus weakened ENSO amplitude during the MH.
5. Cause of meridional advective feedback change
The second leading term that contributed to the decrease in the MLTA tendency during the MH was the meridional advective feedback, (Fig. 4). As shown in Fig. 5, the MH change in the advection of anomalous temperature by mean meridional current in the two models was primarily caused by the change in perturbation temperature gradient , which relied on the change in SSTA itself and meridional SSTA structure. Therefore, reduced SSTA amplitude or increased SSTA meridional scale due to the weakened thermocline feedback would further weaken the meridional advective feedback because the mean meridional current in the EEP always points away from the equator. This indicates that the Bjerknes thermocline feedback was the major contributor to the weakening of ENSO amplitude during the MH, whereas the meridional advective feedback acted as an amplifier to enlarge the above change.
6. Mechanism responsible for the thermocline response
Why did the two models both reproduce a weaker-than-PI thermocline response in the EEP during the MH? Based on a theoretical Sverdrup balance relationship (Neelin 1991; Jin 1997; Li 1997a), the relation between D′ and at the equator is as follows:
where D′ and H1 denote anomalous and mean thermocline depth, ρ0 is the density of seawater, and g is the reduced gravity.
According to Eq. (5), two factors may influence the anomalous zonal thermocline gradient at the equator. One is the mean thermocline depth change and the other is the strength of at the equator. It was found that there was little change in the mean thermocline depth at the equator between PI and MH both for MRI-CGCM2.3.4fa-gvm and MPI-ESM-P (figure not shown). This implies that the zonal-mean thermocline depth was not responsible for the weaker thermocline response and thus weakened ENSO amplitude during the MH, which disagrees with the ideas of Liu et al. (2000) and Otto-Bliesner et al. (2003).
Next, we examined whether there was an obvious change in the strength of at the equator or in the meridional structure of during the MH for the two models. Because ENSO is an air–sea coupled system, the changes in the meridional scales of and SSTA are closely linked to each other. To better reflect the nature of air–sea coupling, we examined the meridional structure of a coupled field represented by the product of and SSTA. This new field was then regressed onto the time series of the product of the EP index and the Niño-3 SSTA index in PI and MH, respectively. As shown in Figs. 9a,b,e,f, the significantly maximum change in the coupled field shifted slightly south of the equator, which could drive southward Ekman ocean current in situ in both models. With respect to the PI period, the coupled field was in general weakened at the equator but significantly strengthened off the equator during the MH (Figs. 9c,g), which led to a minimum of the meridional profile of the difference averaged over 120°E–80°W slightly south (north) of the equator for MRI-CGCM2.3.4fa-gvm (MPI-ESM-P) (Figs. 9d,h). This indicates that the meridional structure of the coupled and SSTA field became wider and flatter during the MH (Figs. 9d,h), which implies a weakened and a decrease in the SSTA amplitude at the equator because of the air–sea coupling and thus a consistent change in the meridional structure of the two fields. According to Eq. (5), a weakened at the equator induced a shallower thermocline response in the EEP. Quantitatively, we calculated the meridional width of the coupled field using the formula in Eq. (3), which is averaged over the equator between 120°E and 80°W as shown in Figs. 9d,h. It was found that the width index of the zonal-averaged × SSTA was increased from 1.72° (2.16°) in PI to 1.89° (2.33°) in MH in the meridional domain of 5°S–5°N and increased from 1.00° (1.19°) to 1.82° (2.94°) in the domain of 10°S–10°N for MRI-CGCM2.3.4fa-gvm (MPI-ESM-P). This implies that the meridional distribution of the coupled and SSTA field became less confined to the equator in the two models, which contributed to less efficiency in forcing anomalous thermocline depth response in the EEP during the MH.
Previous studies have indicated that the climatological mean meridional current plays an important role in determining the ENSO meridional scale (Zhang et al. 2009, 2013; Zhang and Jin 2012). Following Li (1997b) and Chen et al. (2015), an antisymmetric component of the mean meridional current relative to the equator, which is associated with the change in the meridional scale of the SSTA, may be derived as follows:
where V(y) represents the mean meridional current at the latitude of y.
Figure 10 displays the meridional–vertical distribution of the antisymmetric mean meridional current averaged over 160°E–90°W, in which the equatorial part of the Pacific subtropical cell (STC) with poleward surface flow was clearly exhibited both in PI and MH. Compared with the PI period, the surface meridional current became stronger during the MH both for MRI-CGCM2.3.4fa-gvm and MPI-ESM-P, and the overall change in magnitude was greater in the latter than in the former model (Figs. 10c,f). When the mean meridional current or the Pacific STC was strengthened, the equatorial SSTA was more effectively spread away from the equator, which led to a flatter but wider ENSO width during the MH.
To further determine what caused the change in the mean Pacific STC and thus the MLTA tendency during the ENSO developing phase, Fig. 11 shows the meridional profile of April–November mean zonal wind stress averaged over 120°E–80°W in PI and MH as well as their difference. On the one hand, westward wind stress was enhanced during the MH for both MRI-CGCM2.3.4fa-gvm and MPI-ESM-P, especially south of the equator (solid lines in Fig. 11), which might strengthen the mean poleward current by Ekman transport. Note that the MH change in the April–November mean zonal wind stress was mainly attributed to the change in the boreal summer season during May–September (dashed lines in Fig. 11). On the other hand, the convergence of the meridional wind stress (with northward wind stress in the south but southward wind stress in the north) was enhanced for MRI-CGCM2.3.4fa-gvm but slightly weakened for MPI-ESM-P (figure not shown), which might weaken the surface poleward current in the former but strengthen it in the latter model. Because the mean Pacific STC was consistently strengthened during the MH for the two models (Figs. 10c,f), it is reasonable to deduce that the zonal wind stress field played a principal role in controlling the mean meridional current, whereas the effect of the meridional wind stress field was minor.
A key question is what caused the enhancement of the mean zonal wind stress during the MH boreal summer. In response to the MH orbital forcing (Berger 1978), the Northern Hemisphere received more insolation at the top of the atmosphere from May to September (Fig. 12e). This led to a large-scale increase (decrease) in land surface temperatures at the middle and high (low) latitudes, as well as relatively colder SSTs over most of the tropical and South Pacific and Indian Ocean in the two models (Figs. 12a,c). Because of the greater land–ocean thermal contrast, Asian and North African monsoon rainfall was enhanced during the MH boreal summer (Figs. 12b,d). The associated convective heating intensified a large-scale thermally direct east–west circulation (Krishnamurti 1971; Chang and Li 2000) that led to anomalous easterlies over the tropical western and central Pacific (Figs. 12a,c). Moreover, the strengthened Indian monsoon could also induce a Kelvin wave response to the east. The easterly anomalies to the east associated with the Kelvin wave response could further suppress the western North Pacific monsoon through induced anomalous anticyclonic shear (Wu et al. 2009), and a negative precipitation anomaly and an anomalous anticyclone thus appeared over the tropical western North Pacific (Fig. 12). The significant negative correlation between the Indian monsoon and western North Pacific rainfall anomalies was found in previous studies (e.g., Gu et al. 2010). The enhanced trade winds in the tropical Pacific discussed above were consistent with the previous modeling studies by Kitoh and Murakami (2002) and Tian et al. (2017, manuscript submitted to J. Climate).
To sum up in the schematic diagram of Fig. 13, the orbitally induced intensification of Asian and African monsoon rainfall and associated large-scale east–west circulation enhanced the trade winds both south and north of the equator, which further strengthened the mean Pacific STC and thus increased the ENSO meridional scale, leading to both a weaker thermocline–zonal wind stress feedback and meridional advective feedback, which eventually contributed to the weakening of ENSO amplitude during the MH. The contribution of the mean trade wind change to the reduced ENSO amplitude is in agreement with previous studies (Clement et al. 2000; Brown et al. 2008a; Zheng et al. 2008), although the responsible air–sea feedback processes are different.
7. Model–data comparison
The mean state of the MH tropical Pacific SST has been reconstructed from a variety of paleoclimate records. In accord with the model results from PMIP2 and PMIP3 (Figs. 12a,c), some reconstructions showed relative cooling of SSTs near the Galápagos Islands in the EEP from magnesium to calcium ratios in planktonic foraminifera (Koutavas et al. 2002). However, other evidence indicated that the MH SSTs were warmer or little changed from the present day in coastal northern Peru in the eastern Pacific based on mollusk assemblages (Sandweiss et al. 1996) and oxygen isotope in otoliths and planktonic foraminifera (Andrus et al. 2002; Lea et al. 2006), as well as in the tropical western Pacific from fossil corals (Gagan et al. 1998) and foraminifers (Stott et al. 2004). These proxy data are generally consistent with the simulated warming in MRI-CGCM2.3.4fa-gvm (Fig. 12a) but inconsistent with the overall colder SSTs as simulated in MPI-ESM-P (Fig. 12c), which shows an extent of uncertainty between climate models.
Records of interannual variability associated with ENSO include high-resolution natural archives such as corals and lake sediments (Rodbell et al. 1999; Tudhope et al. 2001; Moy et al. 2002). In general, the simulated weaker-than-PI ENSO amplitude during the MH in the PMIP2 and PMIP3 models has been supported by multiple reconstructions. Coral records from Papua New Guinea (Tudhope et al. 2001; McGregor and Gagan 2004) indicate that ENSO activity was 15%–60% weaker than at present during the MH. Similarly, inferences of considerably weaker or absent ENSO events in the early-to-mid-Holocene have also been drawn from pollen data from Australia (Shulmeister and Lees 1995), fossil mollusk shells and geoarchaeological evidence from Peru (Sandweiss et al. 1996; Carré et al. 2014), and condensed and low-frequency laminated deposition in lake sediments from Ecuador (Rodbell et al. 1999; Moy et al. 2002) and the Galápagos Islands (Riedinger et al. 2002). Many other studies that analyzed proxy records also supported the suppression of ENSO activity during the MH (McGlone et al. 1992; Cole 2001; Rein et al. 2005; Koutavas et al. 2006; Conroy et al. 2008; Donders et al. 2008; Cobb et al. 2013; McGregor et al. 2013), and the reconstructed reductions in the MH ENSO variance were found to vary between the western, central, and eastern tropical Pacific (Emile-Geay et al. 2016). As a whole, the overall weakening of ENSO amplitude during the MH in the model results from PMIP2 and PMIP3 is consistent with multiproxy data, although the magnitude of the reduction in simulations (2%–21%; Figs. 1 and 3) was smaller than that in reconstructions (15%–60%).
8. Conclusions and discussion
In this study, the ENSO amplitude change during the MH was investigated using 28 coupled models in the framework of the PMIP2 and PMIP3 simulations. Most models could capture the main features of the observed tropical climate and ENSO characteristics, although several biases still systematically existed, such as a westward extension of the equatorial cold tongue, a double ITCZ, and too regular ENSO events exhibited by only one dominant frequency. Compared with the PI period, 20 out of the 28 models reproduced a weakening of the ENSO amplitude by 2%–21% during the MH. Two representative models, MRI-CGCM2.3.4fa-gvm from PMIP2 and MPI-ESM-P from PMIP3, were then selected to explore the cause of the weakened ENSO amplitude during that time. A mixed layer heat budget analysis was performed to diagnose the main contributors to the composite ENSO amplitude change in the two models. It was found that the weakening of the MH ENSO amplitude was primarily attributed to the decrease in the Bjerknes thermocline feedback (i.e., the advection of anomalous temperature by mean upwelling), whereas the meridional advective feedback (i.e., the advection of anomalous temperature by mean meridional current) acted as an amplifier to enlarge the above change.
A separation of the relative contributions from the mean state and perturbation indicated that the major contributor to the MH changes in the two leading air–sea feedback terms was the perturbation part, not the mean state change. This encouraged us to further investigate the air–sea interaction processes between the atmospheric and oceanic anomalies. Quantitative examination of the strength of the air–sea feedback coefficients showed that the MH reduction in the growth rate associated with the thermocline feedback primarily came from the anomalous thermocline–zonal wind stress relationship. Given a unit zonal wind stress anomaly forcing in the EP, the ocean thermocline response in the EEP was consistently weakened in the two models. According to the Sverdrup balance relationship, the weaker thermocline response during the MH was mainly attributed to the change in the ENSO meridional scale, not the mean thermocline depth change. The meridional structure of the air–sea coupled zonal wind stress and SST anomaly field became less confined to the equator during the MH, which led to less efficiency in forcing anomalous thermocline depth response in the EEP. Further inspection of the mean state change indicates that the wider and flatter ENSO width was caused by the strengthening of the basic-state surface poleward current or the Pacific STC near the equator because of the enhancement of the mean zonal wind stress in the tropical Pacific by Ekman transport, which was eventually attributed to the intensified Asian and African monsoon rainfall and associated large-scale east–west circulation in response to the orbital forcing during the MH.
Note that the simulated weaker ENSO amplitude during the MH in the PMIP2 and PMIP3 models has been supported by multiple paleoclimate records from corals, lake sediments, and pollen data. However, the magnitude of the reduction in the simulations was generally smaller than in the reconstructions. The discrepancy may partly arise from the difference in the reference period. The control run for the MH experiment referred to the PI simulation in the PMIP2 and PMIP3 models, whereas the reconstructed suppression of ENSO activity during the MH was compared with the present condition. The simulated biases for the tropical climate and ENSO characteristics during the PI period may have also been responsible for the underestimated MH change in ENSO amplitude in climate models, possibly because of their insufficient sensitivity to external forcings, their underestimation of internal natural variability, and their lack of important feedbacks among the components of the climate system (e.g., Braconnot et al. 2012a; Hargreaves et al. 2013). On the other hand, there was less consistency between the simulated and reconstructed changes in the mean state of the tropical SST during the MH at the site scale, which shows some uncertainties both in model simulations and in proxy data. Additionally, available proxy data remain scarce, which calls for more reconstruction work using multiple proxies and methods to test model results.
It is worth mentioning that only two models, MRI-CGCM2.3.4fa-gvm and MPI-ESM-P, were selected from the 28 model simulations in PMIP2 and PMIP3 to diagnose the cause of the weakening of the ENSO amplitude during the MH. Although MRI-CGCM2.3.4fa-gvm is a flux-adjusted model, it is still a good choice to serve as a representative of PMIP2 models since the flux-adjusted method has no direct impact on air–sea feedbacks on interannual time scale. In an indirect way, instead, the PI cold bias in the equatorial eastern Pacific is reduced when the flux-adjusted method is applied, which can affect the MH change in estimated MLTA tendency and thus the interannual air–sea coupling to certain extent. For MPI-ESM-P, an obvious error existed in the estimated MLTA tendency change for the El Niño and El Niño minus La Niña composites (Figs. 4d,f), which may primarily arise from the mean SST cold bias in the equatorial eastern Pacific for the PI period. This PI cold bias is inherent to current climate models and stronger for the non-flux-adjusted model MPI-ESM-P than for the flux-adjusted model MRI-CGCM2.3.4fa-gvm (Fig. S1), which affects the MH change in the mean mixed layer temperature in Eq. (1), leading to greater decrease in estimated (Figs. 4d,f, bar 11) than in actual MLTA tendency (Figs. 4d,f, bar 12) for the El Niño and El Niño minus La Niña composites. This implies that the ability of the PMIP2 and PMIP3 models in reproducing the present SST climatology over the tropical Pacific, particularly in the equatorial cold tongue region, may have important influence on the result of the ENSO change during the MH. Thus further efforts should be devoted to improve the model performance on the present-day mean state SST for the next stage of PMIP. In addition, the discrepancy between the change in estimated and actual MLTA tendency for both models may also be caused by other processes as represented by the residual term in Eq. (1). As mentioned before, there was a large spread both in spatial pattern and magnitude among the individual models, with the MH change in the ENSO amplitude varying from −21% to 29% relative to the PI period (Figs. 1 and 3). Although the two models are representative of the majority of models in the sense of the amplitude change based on the criteria for model selection, it is desirable to extend the present diagnostic methodology to all PMIP2 and PMIP3 models to shed light on the mechanisms involved in the most common changes among the models. Moreover, the specification of a constant mixed layer depth is an approximation in the current study, and a better way is to use a temporally and spatially varying mixed layer depth in the heat budget. Finally, because the background mean state change plays an important role in affecting the ENSO amplitude, further analysis is required to understand which fundamental physical processes caused the tropical-mean state changes and which mean state changes were physically reliable (Zhang and Li 2014) during the MH.
We sincerely thank the three anonymous reviewers for their insightful and constructive comments. We also thank the climate modeling groups (listed in Table 1) for making their model output available and Seiji Yukimoto for providing the monthly ocean data of MRI-CGCM2.3.4fa-gvm. This research was jointly supported by the National Natural Science Foundation of China (Grants 41625018 and 41421004), the China National 973 project 2015CB453200, NSF Grant AGS-1565653, and NNSFC Grant 41630423.
Supplemental information related to this paper is available at the Journals Online website: http://dx.doi.org/10.1175/JCLI-D-16-0899.s1.