Responses of Tropical Background State and ENSO Behaviors to Mid-Holocene Forcing Simulated by PMIP3 and PMIP4 Models

The mid-Holocene (MH), a period about 6,000 years ago, provides an opportunity to understand climate change in response to orbital forcing change. Numerical model simulation is an effective tool through which we can study the climate change in the MH, although the climate in the MH can be partly inferred from proxy data. As the Paleoclimate Model Inter-comparison Project phase 4 (PMIP4) recently released the latest simulations for different past climate scenarios, we investigated tropical climate changes, including both the basic state and interannual variability, and tried to find out whether the PMIP phase 3 (PMIP3) and PMIP4 results can be reconciled. Almost all the modelling results show that the inter-hemisphere contrast was enhanced over the tropical Pacific, with warmer and wetter condition on the northern side of the equator and an intensified cross-equatorial flow in the MH than at present, and the annual cycle of the sea-surface temperature (SST) during the MH was reduced. Such background mean state changes arose from the seasonal changes of the solar incident radiation in the MH. In addition to the consistent changes in background mean state changes, some El Niño-Southern Oscillation (ENSO) features, such as the seasonal phase locking feature and periodicity, show consistent changes across the PMIP3 and PMIP4 models, that is, both suites of models exhibit no marked difference in the MH with respect to present-day simulations. In contrast, the modeling results show only agreement on the sign of the ENSO amplitude change (i.e., decrease in the ENSO amplitude in the MH), while the range of reduction varied with model and region. Additionally, the occurrence probability of central Pacific El Niño events increases in the MH, whereas the significance is quite marginal. The modeled changes in the mean state and ENSO serve as a test bed for studying tropical climate system’s response to natural warming, which may provide some insights into understanding climate changes in response to the current anthropogenic warming.


INTRODUCTION
In the mid-Holocene (MH) of~6,000 years before present, the climate system differed greatly with respect to the present-day climate. It is accepted that the root cause leading to the distinctive MH climate that differed from the present-day climate lies in different orbital parameters in the MH epoch (e.g., Kutzbach and Otto-Bliesner 1982;Kutzbach and Liu 1997;Zheng and Braconnot 2013;Zheng and Yu 2013). Understanding the paleoclimate changes, especially the change in the typical warming period-MH period, has profound meaning.
During the MH, the climate in the Earth system experienced obvious changes. For instance, the insolation at high latitudes in the boreal summer was greater in the MH than in present-day climate; hence, a lot of ice sheets in the Northern Hemisphere (NH) melted (Gong et al., 2015;Gierz et al., 2020). The global monsoon areas and precipitation in the MH were enlarged compared with the present-day based on the evidence from pollen data and model simulations (Kohfeld and Harrison 2000;Jiang et al., 2015). Moreover, the tropical Pacific region also witnessed the orbital-induced climate changes in the MH. Previous studies have shown that the mean state, annual cycle and dominant interannual variability in the tropical Pacific, which is called the El Niño-Southern Oscillation (ENSO), exhibited great changes in response to the orbital forcing alteration in the MH. However, these studies also underscored uncertainty regarding the climate changes in the tropical Pacific region in the MH.
Based on the evidence of fossil corals, Gagan (1998) indicated that parts of the tropics were warmer during the MH than at present-day. In contrast, through measuring δ 18 O from proxy data, Koutavas et al. (2006) and Carre et al. (2014) found that the mean sea-surface temperature (SST) during the MH was cooler than that at present-day by about 0.5°o-3°C. In addition to the mean state, ENSO variability in the MH also had some differences from the present-day. A number of proxy records, such as corals, fossil mollusk shells and ice cores, provide evidence for a weakening ENSO variability during the MH, but the magnitude of this weakening varies in a range of 29-80% (e.g., Tudhope et al., 2001;Koutavas et al., 2006;Koutavas and Joanides 2012;Carre et al., 2014;Emile-Geay et al., 2016). However, an interesting study conducted by Cobb (2013) argued that the fossil coral sequences may be too short to reveal the ENSO activity change for a long period of time, and they even argued that one may not easily obtain the conclusion that a significant reduction of ENSO variability is found in the MH. To summarize, there is no consensus on the magnitude of ENSO variability reduction based on the proxy data. This promotes the paleoclimate community to turn to the model simulation's help.
During the past several decades, more and more efforts have been devoted to paleoclimate simulations, e.g., the continuous Paleoclimate Model Inter-comparison Project phase (PMIP) projects. They found that the mean SST in the tropical Pacific was around 1°C cooler than the present-day, especially in the tropical Pacific cold tongue region (e.g., Otto-Bliesner 1999;Liu et al., 2000;Otto-Bliesner et al., 2003;Brown et al., 2006), which supports part of the proxy records. Same as proxy records, plenty of modeling studies showed that the ENSO activity in the MH was reduced, but the magnitude of the reduction varied from 5 to 80% (Otto Bliesner et al., 2003;Zheng et al., 2008;Zheng and Yu 2013;Tian et al., 2017;Chen et al., 2019a;Chen et al., 2019b;Brown et al., 2020). Although most modeling studies show an agreement on the sign of ENSO amplitude change during the MH, there are still large uncertainty regarding the specific magnitude of ENSO amplitude change among these modelling studies.
Recently, the PMIP phase 4 (hereafter PMIP4) simulation results were released (Kageyama et al., 2018), it is necessary to find out whether the PMIP phase 3 (hereafter PMIP3) and PMIP4 results can be reconciled. Here we focused on the mean state, annual cycle and ENSO activity reconstructed by the PMIP4 models over the tropical Pacific during the MH. The simulations of the PMIP3 models were used to be compared with those of the PMIP4 models, and the discrepancies of simulation results between PMIP3 and PMIP4 models were analyzed. The data used in this study is introduced in Section 2. The simulation results of mean state and annual cycle are shown in Section 3, and the ENSO variability is examined in Section 4. We summarize this study in Section 5.

DATA
The simulation for the MH climate period is one of the important experiments in the PMIP, which is also part of Coupled Model Intercomparison Project (CMIP). According to the data ability, 13 PMIP3 models and 12 PMIP4 models were used in this study. The basic information of the 25 models are presented in Table 1, and the boundary conditions for the MH simulation in both PMIP3 and PMIP4 are shown in Table 2. The primary difference between the PMIP3 and the PMIP4 experiments boundary conditions lies in greenhouse gases (GHG) concentration, and the other boundary conditions such as ice-sheets, aerosols and solar constant are identical. The variables used in this study include SST, precipitation, surface wind stress, sea water potential temperature and top of the atmosphere incident shortwave radiation. Due to the different model horizontal resolutions, all data were interpolated onto a regular 1°× 1°(latitude x longitude) grid.

Mean States
The tropical climate states of SST and wind stress during the PI and MH of 13 PMIP3 models and 12 PMIP4 models and their differences in terms of multi-model ensemble mean (MME) are shown in Figure 1. In the present-day climate simulations (Figures 1A,D) and MH simulations ( Figures 1B,E), both PMIP3 and PMIP4 models show similar large-scale distribution features, including the cold tongue in the eastern equatorial Pacific, the warm pool in the western Pacific, and the surface winds converging toward the warm surface temperature region. In general, the trade wind and mean SST appear approximately symmetric with respect to the equator over the western tropical Pacific, while the mean SST exhibits asymmetric distribution with respect to the equator over the eastern tropical Pacific (i.e., warm SST in the northeastern part but cold SST in the southeastern part), matching the cross-equatorial southerly over the eastern Pacific, especially the southerly near the coast of Peru. The difference maps ( Figures 1C,F) show the mean SST cooling to various extent over the entire tropical Pacific in the MH compared to that in the PI. In the MMEs of both PMIP3 and PMIP4 models, the magnitude of the cooling in the tropical Pacific reached about 0.1°-0.6°C in the MH. In general, the decrease in the mean SST in the MH is greater in the western equatorial Pacific than in the eastern equatorial Pacific, indicating a slightly westward extension of the cold tongue. The obvious feature from the difference maps (Figures 1C,F) lies in the inter-hemisphere contrast, that is, the mean SST cooling over the tropical Pacific is more obvious in the Southern Hemisphere (SH) than in the NH, consistent with the cross-equatorial southerly. Previous studies (e.g., Luan et al., 2012;Zhao and Harrison 2012) suggested that such inter-hemisphere contrast may result from the strengthened monsoon precipitation in the NH and weakened monsoon precipitation in the SH, which fundamentally stems from orbital forcing change, that is, the enhanced summer insolation in the NH and reduced summer insolation in the SH. It is noted that such inter-hemisphere contrast is slightly stronger in the MME of the PMIP4 models than of the PMP3 models. Figure 2 shows the spatial distribution of mean precipitation in the PI and MH simulations, as well as their differences. The simulated mean precipitation in the PI simulations resembles that in the observation (not shown), except for the excessive eastward extension of the South Pacific convergence zone (SPCZ), which is associated with the notorious doubleintertropical convergence zone (ITCZ) bias-a prevalent bias in coupled models. It is worth mentioning that such double-ITCZ bias is slightly alleviated in the PI simulations of the PMIP4 models compared to the counterpart in the PMIP3. Next, we analyze the difference of tropical mean precipitation in the MH with respect to the PI. As shown in Figures 2C,F, the difference in mean precipitation in the MH also exhibits marked inter-hemisphere contrast, that is, more precipitation in the NH than in the SH. Such asymmetric mean precipitation difference distribution is consistent with the aforementioned interhemisphere contrast of the mean SST and surface wind.
Again, the magnitude of the difference in the mean precipitation is slightly larger in the PMIP4 than in the PMIP3. Figure 3 shows the equatorial profile of the mean thermocline depth for the MME results from the PI and MH simulations. We can see that the difference in the equatorial mean thermocline during the MH and PI is negligible; this feature is also true for the individual models (not shown). This may be due to the fact that surface flow changes are insignificant for the zonal wind at the equator during the MH, although the southerly wind difference is obvious. The insignificant difference in surface zonal wind may be responsible for the similar equatorial mean thermocline profiles in the two periods.

Annual Cycle
In this subsection, we examine the annual cycle of SST in the equatorial (5 o S-5 o N) Pacific derived from the MMEs of the PMIP3 and PMIP4 models. The annual cycle is characterized by the climatological monthly average after removing the long-term annual mean. As shown in Figures 4A,D, the PI simulations in both PMIP3 and PMIP4 models show an obvious annual cycle of SST in the eastern equatorial Pacific, close to that in the observation (not shown). When comparing the annual cycle result in the PI simulations with the MH simulations, we find that the annual cycle in the eastern equatorial Pacific is weaker in the MH simulations ( Figures 4B,E). Additionally, the peak interval time of the annual cycle in the MH shows a slight difference from that in the PI simulations, that is, the annual cycle of SST in the eastern Pacific reached its positive peak in april and negative peak in August in the PI simulations, whereas the positive peak was slightly shifted toward May and the negative peak, shifted toward July in the MH simulations.
Two factors are responsible for the obvious weakening in the annual cycle of SST in the MH. Following Xie (1994) and An and Choi (2013); An and Choi (2014), the physical reason for the reduced annual cycle amplitude can be explained by the following equation that describes the annual cycle amplitude of SST, where MLT indicates mixed layer temperature, subscript AC indicates the annual cycle after removing the annual mean, c is the phase speed associated with air-sea interaction, h represents mixed layer depth, u 0 indicates mixed layer-averaged mean zonal current, and β denotes the vertical mixing effect on the latent heat flux.
u and v indicate the mean zonal and meridional winds, respectively. ε denotes the Newtonian cooling coefficient. LH and SW are the mean latent heat flux and the perturbed solar radiation, respectively. It is worth noting that the amplitude of SST annual cycle is directly defined as the difference between the maximum and the minimum (e.g., Chen and Jin 2018), which means that the amplitude of SST annual cycle is closely related to the seasonal evolution of the corresponding SST tendency. Considering that the positive SST tendency will be offset by the negative SST tendency in a calendar year, a mathematical treatment is conducted through multiplying both sides of the original mixed layer temperature (MLT) tendency equation by MLT AC . In this way, we obtained Equation 1 and the term z〈MLT 2 AC 〉 zt in the left-hand side of Equation 1 is proportional to the overall SST (or MLT) tendency that determines the annual cycle magnitude. Previous studies (e.g., An andChoi 2013, 2014) pointed out that the annual cycle amplitude of SST is largely determined by the second and third terms on the right-hand side of the The change in the annual cycle of SST during the MH can be traced back to the orbital parameter change (Clement et al., 2000;Braconnot et al., 2012a), that is, the insolation at the top of the atmosphere, which varies between spring and fall, as a result of the orbital-induced alteration during the MH (Berger 1978). As shown in Figure 5, the incoming solar radiation is less in early calendar months (especially February) but more in the boreal summer (especially August) in the equatorial region. According to the second term on the right-hand side of Eq. 1, such change in the incoming solar radiation can cause a weakened annual cycle of SST in the MH. Consequently, the annual cycle of SST shows a decrease in the MH simulations compared to that in the pI. Moreover, the slightly intensified cross-equatorial wind may lead to a slight deepening in the ocean mixed layer (not shown) due to the strengthened stirring effect. The second and third terms on the right-hand side of Eq. 1 indicate that the annual cycle amplitude is inversely proportional to the mixed layer depth h; thus, the slightly deepened mixed layer can cause the increase in the thermal inertia and the reduction of the annual cycle amplitude of SST.
Note that the shading in Figures 1C,F, Figures 2C,F, Figures  4C,F indicates the difference between the MH and PI simulations is significant based on Student's t-test. To summarize, the modelled tropical Pacific mean states during the MH show an agreement across the PMIP3 and PMIP4 models, including the mean SST cooling in the tropical Pacific with severely cooling in the SH, mean precipitation decrease in the SH but slight increase in the NH, insignificant difference in the mean thermocline depth in the equatorial Pacific, and remarkable weakening in the annual cycle of SST in the eastern equatorial Pacific.

EL NINO-SOUTHERN OSCILLATION
In this section, we focus on whether there were some obvious differences in the dominant interannual mode of the ENSO during the MH. The ENSO activity is commonly described by the standard deviation (STD) of SST anomalies. Figure 6 shows the STDs of monthly SST anomalies averaged in the centraleastern equatorial Pacific in each of the PMIP3 and PMIP4 models. First, the majority of the models show that the mature phase of ENSO is phase-locked to the boreal winter (November-January) in the PI simulations, although a few models fail to show such phase-locking feature. When comparing the seasonal phase-locking feature between the PI and MH simulations, we can clearly see that there are nearly no differences in ENSO phase-locking feature.  Figure 7 shows the spatial distribution of the STD of SST anomalies averaged over the boreal winter (November-January) for the MMEs of PMIP3 and PMIP4 models. The difference plots between MH and PI simulations ( Figures 7C,F) show an obvious decrease in the interannual variation of SST anomalies over the central-eastern equatorial Pacific, indicating a reduced ENSO activities during the MH. Additionally, the extent of the reduction is more obvious in the PMIP4 MME than in the PMIP3 MME. Note that in both PMIP3 and PMIP4 models, the stippling denotes that the majority of the models (i.e., exceeding two thirds) show an agreement with the change sign, indicating the majority show suppressed ENSO variability.
To show the spread of ENSO variability changes across the models, we plot the box and whisker chart. As shown in Figure 8, there is large uncertainty regarding the changes in ENSO-related SST anomalies during the MH. For one thing, the spread among the models is somewhat large. Moreover, the magnitude of reduction in MME is small, although the overall variation of SST anomalies shows decreases in different regions, including the Niño3, Niño4, and Niño3.4 regions. For instance, the STDs of the Niño3.4 index are 0.86 and 1.15 for the MMEs of PMIP3 and PMIP4 PI simulations, respectively; and the counterpart in the MH simulations are only 0.81 and 1.05 for the MMEs of the PMIP3 and PMIP4, showing a reduction of about 5.8% and 8.7%, respectively. Recall the relatively large spread among the models, such minor reduction may cause a loss of certainty regarding the modeled changes in ENSO variability. For the Niño1+2 region, a core region for the ENSO, the difference in the variation of SST anomalies between the MH and PI is negligible ( Figure 8D). Since the reductions in the variation of SST anomalies are somewhat small in the MMEs of the PMIP3 and PMIP4 models, we use Figure 9 to clearly show these small changes. Nearly all the scatters that represent the STDs of the SST anomalies over the Niño1+2 region for the PI and MH simulations are located along the diagonal line ( Figure 9A), indicating little change regarding the variation of SST anomalies in the Niño1+2 region. For the regions to the west (Niño3, Niño3.4, and Niño4), the number of the models that show weaker STDs of SST anomalies in the MH than in the PI increases. Such increase may be partly due to the fact that most of the coupled To summarize the changes of ENSO-related SST anomalies, we present the violin plots for all the PMIP3 and PMIP4 models. Figure 10 shows the probability density of each relevant index. The results are consistent with the aforementioned results building on the PMIP3 or PMIP4 models only, that is, the ENSO-related SST anomalies in the Niño4, Niño3.4, and Niño3 regions exhibit a weakening to some extent in the MH compared to the PI, where the uncertainty remains large.
So far, various mechanisms have been put forward to account for the suppressed ENSO amplitude in the MH. For instance, Clement et al. (2000) indicated the reduced ENSO variance is regarded as a response to the orbital-induced change of the tropical annual cycle. Liu et al. (2000) suggested the deepened mean equatorial thermocline and the strengthened Asian summer monsoon may lead the suppression of ENSO intensity in MH. Moreover, Zheng et al. (2008) pointed out that the reduction of ENSO intensity is interpreted as the influence of strengthened easterly trade wind in tropic Pacific. Chen et al. (2019a); Chen et al. (2019b) suggested that the changes in the mean Pacific subtropical cell (STC) is the key factor responsible for the reduced ENSO intensity in the MH. In a word, the real physical reason for the inhibition of ENSO activity in the MH period is complex, which still needs some more indepth investigation in the future.
Next, we examine whether there is any change in terms of ENSO periodicity. Figure 11 suggests almost no change regarding the spectrum, for either PMIP3 or PMIP4 MME. Such results indicate that the ENSO period is not sensitive to the orbital forcing change in the MH, at least in the current climate models.
In addition to ENSO amplitude and periodicity, we also investigate the change in ENSO diversity in response to the orbital forcing. Several studies uncovered a new flavor of El Niño events in recent years, which is characterized by warm SST anomalies in the central Pacific (CP) and usually called the CP El Niño, El Niño Modoki, warm pool El Niño, or dateline El Niño (Larkin and Harrison 2005;Ashok et al., 2007;Kug et al., 2009;Yu and Kao 2009;. Based on the criteria proposed by Kug et al. (2009), we first compare the SST anomalies in the boreal winter averaged in the Niño3 region with that in the Niño4 region for certain El Niño event, and then label this event as the eastern Pacific (EP) El Niño (CP El Niño), if the Niño3 index is greater (less) than the Niño4 index. Figure 12 shows the ratio of the CP El Niño number to the EP El Niño number for each model and its corresponding MME. Clearly, a divergent performance in simulating the change of the ratio exist across the models. Half of the PMIP3 models (six out of the 13 models) show the ratio of the incidence of CP El Niño events to the EP El Niño events in the MH outstrips the that in the PI, and five among them exhibit a significant increase in the occurrence of CP El Niño during the MH at the 95% confidence level. For the PMIP4 models, eight out of the 12 models record more incidence of CP El Niño events during the MH, and five among them have enhanced ratios exceeding the 95% confidence level. Then, we obtain the MMEs for the PMIP3 and PMIP4 models (the rightmost columns of Figures   12A,B). The increase of CP El Niño events during the MH is minor and insignificant for the PMIP4 MME. Although the increase of CP El Niño events during the MH is relatively large for the PMIP3 MME, it is also not statistically significant at the 95% significant level. Overall, the model results show a larger occurrence of CP El Niño in the MH than in the PI, though the significance is quite marginal.

SUMMARY
Increasing attention has been given to understand tropical climate change, especially in terms of ENSO behaviors, in response to global warming; however, some puzzles remain Frontiers in Earth Science | www.frontiersin.org (e.g., Collins et al., 2010;Chen et al., 2015;Chen et al., 2017). One way to gain confidence in future projections is to evaluate current models' performances in simulating past climate changes (Braconnot et al., 2012a;Braconnot et al., 2012b), because some proxy datasets can be used as reference. Here, we investigated the response of tropical Pacific climate to the orbital forcing in the MH, an epoch about 6,000 years ago, by analyzing model results. Although the climate in MH can be partly inferred from proxy data, model simulation is an effective tool through which we can study the climate change in the MH. Using recently released PMIP4 outputs along with PMIP3 archives, we examined tropical climate changes, including both basic state and interannual variability; we then attempted to find out whether there is any aspect that can be reconciled for the PMIP3 and PMIP4 model results. The main conclusions are summarized here. 1) All the PMIP3 and PMIP4 models show the enhancement of the inter-hemisphere contrast over the tropical Pacific, including the warmer and wetter condition north of the equator and an intensified cross-equatorial flow in the MH than at present. This difference can be attributed to the change in the orbital forcing, or more straightforwardly, the change in the solar incident radiation.
2) Most models show a weakening in the annual cycle of SST during the MH. The relevant air-sea processes may be more complicated, although such background mean state changes can also be linked to the orbital forcing change in the MH. Specifically, the weakened annual cycle amplitude of SST primarily arises from the seasonal change in the incoming solar radiation during the MH, and the slight shoaling in the mixed layer depth plays a secondary role.
3) The majority of the PMIP3 and PMIP4 models show that ENSO amplitude was suppressed during the MH, while the range of reduction varies with model and region. When synthesizing the overall probability distribution of the changes in ENSO indices, we found that the results remain large uncertainty, which varies with space and pattern. Additionally, the occurrence probability of CP El Niño events increased in the MH, whereas the significance is quite marginal. Interestingly, the other ENSO behaviors, such as the seasonal phase-locking feature and periodicity, exhibit no obvious differences between the MH and PI simulations.  Frontiers in Earth Science | www.frontiersin.org March 2022 | Volume 10 | Article 853577