The hydrological cycle and ocean circulation of the Maritime Continent in the Pliocene: results from PlioMIP2

The Maritime Continent (MC) forms the western boundary of the tropical Pacific Ocean, and relatively small changes in this region can impact the climate locally and remotely. In the mid-Piacenzian warm period of the Pliocene (mPWP; 3.264 to 3.025 Ma) atmospheric CO₂ concentrations were ∼ 400 ppm, and the subaerial Sunda and Sahul shelves made the land–sea distribution of the MC different to today. Topographic changes and elevated levels of CO₂, combined with other forcings, are therefore expected to have driven a substantial climate signal in the MC region at this time. By using the results from the Pliocene Model Intercomparison Project Phase 2 (PlioMIP2), we study the mean climatic features of the MC in the mPWP and changes in Indonesian Throughflow (ITF) with respect to the preindustrial. Results show a warmer and wetter mPWP climate of the MC and lower sea surface salinity in the surrounding ocean compared with the preindustrial. Furthermore, we quantify the volume transfer through the ITF; although the ITF may be expected to be hindered by the subaerial shelves, 10 out of 15 models show an increased volume transport compared with the preindustrial.

In order to avoid undue influence from closely related models that are present in the PlioMIP2 ensemble, we introduce a new metric, the multi-cluster mean (MCM), which is based on cluster analysis of the individual models. We study the effect that the choice of MCM versus the more traditional analysis of multi-model mean (MMM) and individual models has on the discrepancy between model results and data. We find that models, which reproduce modern MC climate well, are not always good at simulating the mPWP climate anomaly of the MC. By comparing with individual models, the MMM and MCM reproduce the preindustrial sea surface temperature (SST) of the reanalysis better than most individual models and produce less discrepancy with reconstructed sea surface temperature anomalies (SSTA) than most individual models in the MC. In addition, the clusters reveal spatial signals that are not captured by the MMM, so that the MCM provides us with a new way to explore the results from model ensembles that include similar models.

1 Introduction

The Maritime Continent (MC; 20^∘ N–10^∘ S and 90–150^∘ E) consists of more than 22 000 islands and lies in the warmest ocean region of the west Pacific Warm Pool (WPWP; Ramage, 1968; Yoneyama and Zhang, 2020). The near-surface air temperature in most areas of the MC is higher than 27 ^∘C throughout the year (Li et al., 2018a). High sea surface temperature (SST) leads to substantial moisture flux into the atmosphere. High levels of moisture and energy characterise the unstable climate in this region and make it “the world's strongest atmospheric convection center” (Yoneyama and Zhang, 2020). Based on satellite data, it has been established that the MC is the region that produces the largest amount of precipitation in the world (Adler et al., 2017; Yamanaka et al., 2018). A large amount of rainfall releases vast quantities of latent heat into the atmosphere (1 mm yr⁻¹ equals 0.08 W m⁻²), which is an important driver of global atmospheric circulation (Trenberth et al., 2009; Yamanaka et al., 2018). According to these characteristics, the MC has been recognised as the “boiler box” of the global climate system (Ramage, 1968; Neale and Slingo, 2003).

The MC also impacts the oceanic hydrological system in that it links the Pacific Ocean and the Indian Ocean via the Indonesian Throughflow (ITF), which is a key component of global ocean circulation. In general, the ITF transports warm and comparably fresh water from the Pacific Ocean into the Indian Ocean. Dramatic tectonic processes and sea level changes have shaped this gateway that plays an important role in influencing climate both locally and in other regions of the world (Gordon, 2005b; Tillinger, 2011; Yoneyama and Zhang, 2020). According to the International Nusantara Stratification and Transport (INSTANT) program, from 2004 to 2006 (Gordon et al., 2010; Sprintall et al., 2009) the ITF transported about 15 Sv of water from the Pacific Ocean into the Indian Ocean; the heat flux (i.e. the effective temperature of the total throughflow) of the total ITF export was 17.6 ^∘C during these 3 years (Sprintall et al., 2009).

The MC acts as a source of sensible heat for the atmosphere over land, and the geographical configuration of the MC redistributes heat in the western Pacific and the Indian Ocean; both the effect of sensible heat and the geography of the MC are crucial for the onset of the Walker Circulation (Dayem et al., 2007). The variation in the Walker Circulation leads to periodic phenomena that modulate the climate of the MC, including the El Niño–Southern Oscillation (ENSO), the Indian Ocean Dipole (IOD) and the Pacific Decadal Oscillation (PDO). The MC forms the western boundary of the ENSO phenomenon. On the interannual scale, ENSO and IOD play important roles in modulating rainfall and the ITF transport by affecting zonal winds (Wang, 2019). During El Niño years, the easterly trade winds in the Pacific are weakened. As a result, rainfall is decreased, and the ITF is weakened (Zhang et al., 2016; Yamanaka et al., 2018), with a peak-to-trough amplitude of about 5 Sv (Feng et al., 2018). During La Niña years, the situation is reversed. The positive phase of the IOD normally occurs during the autumn of an El Niño developing year, which offsets the effects of El Niño on the ITF (Wang, 2019; Sprintall and Révelard, 2014). On decadal timescales, the PDO plays a dominant role in regulating the ITF (Wang, 2019; Li et al., 2018b). Similar to the effect of ENSO, in the warm phase of the PDO, easterly trade winds become weaker, causing the ITF transport to weaken, and vice versa (Wang, 2019). A modelling study from Tan et al. (2022) highlighted the importance of the geometry of the MC and found that the shallow opening of the MC can trigger an active Pacific meridional overturning circulation and enhance the upwelling strength along the Equator from the central to eastern Pacific.

It is possible to obtain an atlas of the future warming world from climate model simulations (e.g. Iturbide et al., 2021; Gutiérrez et al., 2021). However, there is no direct evidence for us to validate how well climate models simulate the future. The mid-Piacenzian warm period (mPWP) is an interval in the Pliocene between 3.264 and 3.025 Ma, which has also been described as mid-Pliocene warm period in previous publications (Dowsett et al., 2016). In this study, we use mPWP as the abbreviation for this warm period. According to the Intergovernmental Panel on Climate Change (IPCC) Sixth Assessment Report (AR6) Working Group I (WGI) Technical Summary (Arias et al., 2021), the mPWP is very likely 2.5 to 4.0 ^∘C warmer than the preindustrial, and the concentration of CO₂ during this period is very likely in the range from 360 to 420 ppmv, which is substantially higher than for the preindustrial (280 ppmv). In terms of temperature and geography, the Earth system of the mPWP is similar to projections for the end of the 21st century (Dowsett et al., 2016). Consequently, the mPWP is considered to be one of the most recent analogues to future warm conditions on a broad-scale perspective (Salzmann et al., 2009; Burke et al., 2018).

During the Pliocene epoch, the Indonesian gateway reorganised due to tectonic processes, and the ITF was restricted progressively (Karas et al., 2009; Auer et al., 2019). This process contributed to variations in the climate throughout the Pliocene. For example, during the mid to late Piacenzian (3.3–2.6 Ma), the reorganisation of the ITF led to a strengthening of the South Asian summer monsoon and acted as a precursor to east African aridification (Cane and Molnar, 2001; Sarathchandraprasad et al., 2021). By using model results, we can explore the mPWP climate in this region. This effort is supported by geological records that provide a variety of evidence towards reconstructing the climate over the MC during the Pliocene. For example, deep-sea deposit samples of the Pliocene can be utilised to reconstruct the SST (Wara et al., 2005; Zhang et al., 2014; O'Brien et al., 2014; Fedorov et al., 2015; McClymont et al., 2020). By comparing results of simulations and geological records, we can validate climate models and gain confidence in the projections of future climate if the models perform well for past warm climates.

The Pliocene has been broadly studied; previous work by the Pliocene Research, Interpretation and Synoptic Mapping (PRISM; e.g. Dowsett et al., 2013, 2016) and the Pliocene Model Intercomparison Project (PlioMIP; e.g. Haywood et al., 2010, 2016a, 2020) provide us with a large number of palaeoclimate reconstructions and simulations of environmental and climatic information. The second phase of PlioMIP (PlioMIP2) that this study relates to is briefly described in Sect. 2 of this paper. Large-scale climate features of the mPWP, as derived from PlioMIP1 and PlioMIP2, are described in Haywood et al. (2013) and Haywood et al. (2020). On a global scale, the simulated annually averaged surface air temperature is elevated, with models suggesting an anomaly above the respective conditions of the preindustrial that ranges between 1.8 and 3.6 ^∘C, based on PlioMIP1 (coupled atmosphere–ocean climate models only), and between 1.7 and 5.2 ^∘C, based on PlioMIP2. The annual mean total precipitation rates are higher by 0.09–0.18 mm d⁻¹ (PlioMIP1; coupled atmosphere–ocean climate models only) and 0.07–0.37 mm d⁻¹ (PlioMIP2) in the mPWP. All the models from PlioMIP1 and PlioMIP2 simulated a clear polar amplification of the warming in the mPWP and show that the warming magnitude in the MC is less than the global average. Regarding precipitation changes in the region, the PlioMIP1 ensemble shows an increase over most parts of the MC, especially in the tropics (Fig. 3b in Haywood et al., 2013). In contrast to PlioMIP1, PlioMIP2 shows spatially heterogeneous changes in precipitation over the MC (Fig. 5b in Haywood et al., 2020). Moreover, there is an increase in precipitation minus evaporation over exposed continents across this region (Fig. 1 in Feng et al., 2022). PlioMIP1 studies have shown that the tropical overturning circulation slows down in the mPWP due to the changes in SST and that the ascending branch of the Walker Circulation over the MC weakens (Corvec and Fletcher, 2017). In terms of climate variability, ENSO exists in the modelled mPWP, but the amplitude of ENSO is weaker than in the preindustrial. This inference is consistent from PlioMIP1 (Brierley, 2015) to PlioMIP2 (Oldeman et al., 2021). However, the mean zonal SST gradient does not decrease consistently across models in PlioMIP2 (Oldeman et al., 2021). Although the MC is important for the regional- and large-scale climate, there are relatively few studies of the mPWP focusing on the MC (e.g. Smith et al., 2020), especially those studies that use state-of-the-art models.

In this paper, we investigate the climate of the MC in the mPWP by employing PlioMIP2 simulations. We study large-scale climate patterns (SST, the net fresh water flux at the surface (precipitation minus evaporation; P−E), salinity at the sea surface (SOS) and wind stress) of the MC and investigate ocean circulation and strength of the ITF in the mPWP. Moreover, we assess the ability of the PlioMIP2 models to simulate the mPWP climate of the MC by cross-validating simulations and models with proxy reconstructions and reanalysis data. We investigate differences in predictions by individual models with a hierarchical clustering method.

The rest of the paper is organised as follows. Section 2 describes the PlioMIP2 models, experiment design and simulation data that are relevant to this study. Section 3 addresses the following main questions:

• Q1

  How do PlioMIP2 models perform in simulating the mPWP and the preindustrial climate of the MC (Sect. 3.1)?

• Q2

  What are the characteristics of the climate (SST, P−E, SOS and wind stress) of the MC in the mPWP in the PlioMIP2 simulations (Sect. 3.2)?

• Q3

  Did the volume transport via the ITF intensify in the mPWP, and what were the characteristics of this throughflow (Sect. 3.3)?

• Q4

  How can we best emphasise results from individual models that may otherwise disappear in the multi-model ensemble mean (MMM)? How can we account for the effect of any bias duplication caused by models from the same “family” within our analysis (Sect. 3.4)?

Finally, we discuss the drivers for changes in the ITF and of the climate of the MC (Sect. 4.1), we summarise performance of individual models and of the MMM in simulating mPWP and preindustrial climate (Sect. 4.2), and we address similarities and differences between results by individual models (Sect. 4.3).

2 PlioMIP2

In this section, we describe the PlioMIP2 models, experiments and the PlioMIP2 project. In Sect. 2.1, we introduce the PlioMIP2 models and the development of models which participated in the PlioMIP phase 1 and phase 2. Then, the two core experiments of PlioMIP2 adopted in this study are described in terms of boundary conditions (Sect. 2.2).

2.1 PlioMIP2 models

PlioMIP (Haywood et al., 2011) is a project aimed at studying climate and environments of the late Pliocene. PlioMIP is part of the Palaeoclimate Modelling Intercomparison Project (PMIP). In consideration of the development of models and proxy data in the past years, this project has progressed to the second phase (PlioMIP2; Haywood et al., 2016b, 2020). In this ongoing phase, there are more models (17 models in comparison to 11 models that participated in PlioMIP phase 1; Table 1). Most of the newly participating models are descendants of existing models. For example, CCSM4 participated in PlioMIP1. The models derived from it, namely the Community Earth System Models of CCSM4-Utr, CCSM4-UoT, CESM1.2 and CESM2, participated in the second phase; IPSL-CM5A from the Institute Pierre-Simon Laplace (IPSL) participated in PlioMIP and then joined PlioMIP2 with IPSL-CM5A2 and IPSL-CM6A-LR; the Norwegian Earth System Model (NorESM-L), which has been developed from CCSM4, participated in the first phase and then participated in the second phase with NorESM1-F. The atmosphere-only model HadAM3 (third Hadley Centre Atmospheric Model) participated in PlioMIP1 as a component of the UK Hadley Centre atmosphere–ocean general circulation model HadCM3; in PlioMIP2, it made contributions via the CMIP6-class model, which is the UK Met Office HadGEM3-GC31-LL (Williams et al., 2021), and via HadCM3 (Hunter et al., 2019). In terms of the atmospheric and oceanic components, new model versions were developed in that the parameterisation changed (CCSM4-UoT; Peltier and Vettoretti, 2014, and CCSM4-Utr use a simplified version of ocean mixing scheme from CCSM4 with parameterisation changes. Moreover, NorESM1-F updated the parameterisation from NorESM-L; model components were updated (CESM1.2 and CESM2 adopted the atmospheric components of CAM5 and CAM6); or the resolution has been enhanced (IPSL models; NorESM1-F applied the tripolar grid in the ocean component compared with NorESM-L). All of these models have run the core mid-Pliocene eoi400 experiment (Eoi⁴⁰⁰) and provide the preindustrial experiment (E²⁸⁰) as a control simulation. For more details on the 16 PlioMIP2 models, see Haywood et al. (2020).

Neale et al. (2010)Danabasoglu et al. (2012)Feng et al. (2020a)Baatsen et al. (2022)Chandan and Peltier (2017)Neale et al. (2010)Feng et al. (2020a)Danabasoglu et al. (2020)Danabasoglu et al. (2020)Feng et al. (2020a)Roeckner et al. (2003)Marsland et al. (2003)Stepanek et al. (2020)Madec et al. (2017)Zhang et al. (2021a)Kelley et al. (2020)Kelley et al. (2020)Hunter et al. (2019)Williams et al. (2021)Hourdin et al. (2013)Madec et al. (2017)Tan et al. (2020)Hourdin et al. (2013); Sepulchre et al. (2020)Sepulchre et al. (2020)Tan et al. (2020)Hourdin et al. (2020); Boucher et al. (2020)Hourdin et al. (2020)Chan and Abe-Ouchi (2020)Kamae et al. (2016)Li et al. (2020)Li et al. (2020)

Table 1 Details of the 16 PlioMIP2 models (Haywood et al., 2020) and of the HadGEM3 model used in this study. NA stands for not available.

2.2 PlioMIP2 experimental design

The global-scale and gridded palaeo reconstruction data from PRISM provide data for mPWP palaeogeography and land ice to enable the preparation of reliable boundary conditions and to carry out climate simulations of the mPWP. For PlioMIP phase 1, experimental boundary conditions for the mPWP have been built based on PRISM3D (Dowsett et al., 2012). In PlioMIP2, new boundary conditions have been adopted, which are derived from the state-of-the-art reconstruction PRISM4 (Dowsett et al., 2016).

In this study, we adopt the E²⁸⁰ simulation as a control experiment, which is forced with preindustrial conditions, and the Eoi⁴⁰⁰ simulation is used as a sensitivity experiment for mPWP conditions, which is forced with all mPWP boundary conditions. These conditions are topography, ice sheets, soil, vegetation, lakes, land–sea mask and a concentration of 400 ppm of CO₂ in the atmosphere. Simulation naming follows the rules established by Haywood et al. (2016a). The numbers 280 and 400 indicate the CO₂ concentration; the character “o” indicates changes to orography, bathymetry, land–sea mask, vegetation, lakes and soils from the preindustrial model set-up towards mPWP conditions in regions free of ice sheets today; and the character “i” indicates such changes in regions of current ice sheets. Figure 1 shows the land–sea masks used in the E²⁸⁰ and the Eoi⁴⁰⁰ experiments indicated by black lines and colours, respectively, and illustrates the location of the MC. From this map, it is clear that in the mPWP, several modern ocean gateways were absent, such as the Bering Strait, the waterbody connecting Canada and Greenland (Canadian Arctic Archipelago), and some straits in the MC. Details of the experimental design are described by Haywood et al. (2016a). The initial conditions of ocean salinity are either derived from Levitus and Boyer (1994), from an equilibrium state of the modern (control) simulation or from the end of the PlioMIP1 experiment (Haywood et al., 2011). Note that HadGEM3 did not change the land–sea mask in the Eoi⁴⁰⁰ experiment, and therefore, for this model, the various gateways are the same as in E²⁸⁰. Data from the experiments E²⁸⁰ and Eoi⁴⁰⁰ used in this study are from the PlioMIP2 data set, which is available at https://esgf-node.llnl.gov/search/cmip6/ (last access: 18 October 2023) and https://geology.er.usgs.gov/egpsc/prism/7.2_pliomip2_data.html (last access: 18 October 2023). In this study, the results are calculated across the final 100 model years of each simulation. Data have been regridded to a 1^∘ × 1^∘ grid in order to calculate MMM and multi-cluster mean (MCM); results for individual models are illustrated based on the original grids. The MC regional average value in this paper is calculated from the region of 25^∘ N–25^∘ S and 80–170^∘ E.

Figure 1 PRISM4 Pliocene land–sea mask. The modern continental outline is indicated with black lines; blue and green shadings indicate ocean and land in the Eoi⁴⁰⁰ experiment. The red rectangle indicates the location of the Timor passage of the Indonesian Throughflow (ITF) that we refer to in Sect. 3.3. The red arrow denotes the position and direction of the ITF illustratively. The Sunda and Sahul shelves are illustrated in this figure via white labels. The continental outline is constructed from Dowsett et al. (2016) and retrodicted with a set of procedures based on the palaeogeographic maps of Markwick (2007).

3 Results

In Sect. 3.1, model results are compared to proxy reconstructions and the Hadley Centre Sea Ice And Sea Surface Temperature reanalysis version 1 (HadISST1) data set to evaluate model performance for reproducing mPWP and modern climate. In addition, the multi-model mean method (MMM) is evaluated in comparison to individual models. We illustrate the most relevant patterns of mean climate for the MC in mPWP relative to the preindustrial in terms of SST, precipitation minus evaporation, wind stress and salinity at the ocean surface (Sect. 3.2). Then we analyse oceanic circulation in the MC (Sect. 3.3). In Sect. 3.4, we introduce a new metric, the multi-cluster mean (MCM), which takes the weight of models from the different clustering branches into consideration to better illustrate the results from multiple models. Finally, we present a comparison between the performance of the MCM, the MMM and individual models.

3.1 Model skill at simulating the preindustrial and mPWP MC climate

In this section, we investigate how well a model performs in reproducing preindustrial climate (Sect. 3.1.1) and mPWP climate (Sect. 3.1.2). Our aim is to address the following questions: (1) are the models that perform well in simulating preindustrial climate of the MC also good at simulating mPWP climate of the MC? (2) Does the multi-model ensemble have improved performance for the MC compared to individual models, as is often the case in weather and climate prediction (e.g. Tebaldi and Knutti, 2007)?

3.1.1 Comparison with HadISST1

Observational data and reconstructed climatic information from proxy data can offer a way to evaluate simulation results. Here, we adopt HadISST1 data set (Rayner, 2003) from January 1870 to January 1900 to evaluate the model performance in reproducing preindustrial climate. The HadISST1 data set is available at https://www.metoffice.gov.uk/hadobs/hadisst/ (last access: 19 October 2023).

Figure 2 shows the spatial distribution pattern of the preindustrial SSTs over the MC from the HadISST1 reanalysed SSTs and from the PlioMIP2 MMM. The MMM reproduces the SST zonal gradient from the Equator to higher latitudes, but simulated SSTs are too warm around the tropical islands and too cold in the northwestern Pacific Ocean when compared to the reanalysis data. By calculating the root mean square error (RMSE) between reanalysed SSTs and model SSTs over the MC, we quantify the model and observation discrepancy for the modern period, which is shown in Fig. S1 and Table S1 in the Supplement. The mean discrepancy in the MC between the MMM SSTs and HadISST1 SSTs is 0.51 ^∘C, which is lower (better) than 13 out of 17 models. By that metric, the MMM better reproduces preindustrial climate in this region than most of the individual PlioMIP2 models. The MMM metric also shows better agreement with observations than individual simulations in many other studies. For example, Williams et al. (2022) show that for the Deep-Time Model Intercomparison Project (DeepMIP) model ensemble, the MMM is the best estimate of the present-day precipitation in terms of spatial patterns relative to individual DeepMIP models.

Figure 2 The annual mean sea surface temperature (SST) over the MC from the (a) PlioMIP2 multi-model ensemble mean (MMM) for the preindustrial (E²⁸⁰), (b) HadISST1 reanalysed SSTs between January 1870 to January 1900 and (c) their difference. White shading indicates the grid boxes that are land for all the models. In panel (a), grid boxes outside the black lines are ocean grid boxes for all the models. Stippled areas show significant differences, according to a one-sample t test, at the 90 % confidence level.

In terms of individual models, CCSM4 shows the lowest discrepancy of any individual model (0.39 ^∘C); IPSL-CM5A and HadCM3 show the highest discrepancies (1.29 ^∘C). In order to evaluate the spatial patterns of the simulated SSTs, the discrepancy between the MMM SSTs and the HadISST1 SSTs is also calculated after bias adjustment (i.e. removing the mean SST over the MC). The values are as shown in Fig. S1.

3.1.2 Comparison with proxies

For the mPWP, the data of multi-proxy reconstructed SST anomalies (SSTAs) compiled by McClymont et al. (2020) are available for a narrow time slice centred on 3.205 Ma. These data (McClymont et al., 2020) are compared with the simulated mPWP SSTAs in this study. Data and the full details for reconstructed SSTAs are available at https://doi.org/10.1594/PANGAEA.911847 and from McClymont et al. (2020). There are four sites with available data over the MC, as shown in Fig. 4a. Sites ODP214 and ODP763 are located at the eastern Indian Ocean; site ODP806 is located at the western tropical Pacific Ocean; and site ODP1143 is located near the islands of the MC. The SSTA for ODP1143 is reconstructed from U${}_{\mathrm{37}}^{K^{\prime }}$ with the alkenone-derived (BAYSPLINE) calibration (Tierney and Tingley, 2018); the SSTAs for ODP763, ODP214 and ODP806 are reconstructed from planktonic foraminifera Trilobatus sacculifer Mg $/$ Ca, with the Mg $/$ Ca-derived (BAYMAG) calibration (Tierney et al., 2019a). Figure 3a shows reconstructed and modelled SSTAs for every site individually. Model results for each site are provided as the average of the nine grid boxes centred on the site in order to reduce local bias.

Figure 3 Discrepancy between modelled and reconstructed SST anomalies (SSTAs) for the mid-Piacenzian warm period (mPWP). (a) Comparison between simulated and reconstructed SSTAs for four sites in the MC. The location of these four sites can be seen in Fig. 4a. The vertical orange bars and shading indicate the double-standard errors (±2σ) of the reconstructed SSTs. The horizontal orange bars and the asterisk indicate the reconstructed SSTAs of the site, individually. The short horizontal black bars among the results of the models indicate the MMM of the simulated SSTAs. Triangles refer to the results by different clusters (coloured) and to the MCM (black). (b) Box plot of model to proxy discrepancies over the MC for every model. Vertical lines inside the boxes are the median values (E) of the probability density function of the discrepancies derived via Monte Carlo simulations, which is explained in Sect. 3.1.2. Lower and upper box boundaries indicate the interval at which the model and proxy discrepancies are likely to be (probability range from 66 % to 100 %). Lower and upper error bars indicate the interval where the model and proxy discrepancies are most likely to be (probability range from 90 % to 100 %). Models are ordered by their value of E. Also shown are the E values for the five clusters and the MCM.

Figure 4 The MMM anomalies over the MC in the Eoi⁴⁰⁰ simulation relative to the E²⁸⁰ simulation for (a) sea surface temperature (SST), (b) precipitation minus evaporation (P−E), (c) zonal surface wind stress (TAUU) and (d) salinity at the ocean surface (SOS). In panel (a), circles mark the location of the SSTA proxies. Their colour represents the reconstructed SSTAs, the value of which is also shown in Fig. 3a, and the name of each site is indicated on the top of each circle. White shading and black lines in panels (a) and (d) are as in Fig. 2. In panels (b) and (c), hatches indicate land grid boxes in the Eoi⁴⁰⁰ experiment, and the black lines indicate the continental outlines of the preindustrial experiment. In panel (b), filled circles denote proxy data in terms of wetter or no-change conditions. Details of each numbered site are shown in Table S2. In panel (c), the green shading indicates westward anomaly, and the purple shading indicates eastward anomaly. Stippled areas show significant multi-model differences, according to a t test, at the 90 % confidence level.

There are large standard errors (σ) for the proxy reconstructed SSTAs, with ±2σ shown with the orange shading in Fig. 3a. Models align well with the reconstructed SSTAs at site ODP806 in the WPWP and at site ODP1143 in the South China Sea. For the sites ODP763, ODP806 and ODP1143, all of the model results are within the uncertainties in the reconstructions. For the site ODP214, a few (two) of models suggest comparably warm mPWP SSTAs that are outside the uncertainty range of the reconstruction. We note that there is a persistent mismatch between Mg $/$ Ca-derived (BAYMAG) and alkenone-derived (BAYSPLINE) SSTAs. McClymont et al. (2020) show that generally, and with few exceptions, Mg $/$ Ca-derived SSTAs are lower than their alkenone-based counterparts. In particular at lower latitudes, there is a substantial mismatch. For example, at site ODP1143, which has both UK37 and Mg $/$ Ca data available, Mg $/$ Ca suggests a negative SSTA, but a result from this proxy is not supported by alkenones. For the Mg $/$ Ca data at sites ODP763, ODP214 and ODP806, the Mg $/$ Ca-derived SSTAs are all cooler than the MMM SSTA. In particular, for site ODP763 and for site ODP214, all models simulate warmer SSTAs than the reconstructions, with a larger distance to the reconstructions for site ODP214. For this site in the Indian Ocean, the proxy data show the mPWP SST that is cooler than in the preindustrial, but all the models show warmer SSTAs. The cold biases of the Mg $/$ Ca-derived SSTAs may be attributed to additional environmental factors, such as the seasonality and the calcification depth of foraminifera (Hertzberg and Schmidt, 2013; Hönisch et al., 2013; Tierney et al., 2019b; McClymont et al., 2020), which can explain larger biases in the tropics where the thermocline is relatively deeper than at higher latitudes. Moreover, the foraminifera Trilobatus sacculifer adopted in ODP214 and ODP763 broadly derives larger negative SST anomalies compared to another surface-dwelling warm water foraminifera, Globigerinoides ruber, in the study of McClymont et al. (2020), which aligns with results from some studies (Fairbanks et al., 1980; Ravelo and Fairbanks, 1992; Curry et al., 1983) which show that, in the tropics, Trilobatus sacculifer is often found in a slightly deeper habitat than Globigerinoides ruber (Tierney et al., 2019b). Furthermore, studies such as Zhang et al. (2014) and Smith et al. (2020) show reconstructed SSTAs from organic geochemical proxies and also suggest that the SSTAs in the site ODP806 and region near the ITF outlet could potentially be warmer.

In order to quantify discrepancies between proxy reconstructed and simulated SSTAs from individual models, we use Monte Carlo simulations, following a similar method to that used in Kageyama et al. (2021). The results are shown in Fig. 3b and sorted in ascending order of the model and proxy discrepancy median value E. Since there are uncertainties for proxy data, we generate an array containing 1000 random SSTAs, using a uniform range distribution derived from the reconstructed SSTAs and ±2σ for every site. Uncertainties in the simulated SSTAs are derived over the last 100 years of annual mean SSTAs. By subtracting the reconstructed SSTA arrays from the 100 simulated SSTAs for every site separately and combining the arrays for the four sites, we derive arrays of model and proxy discrepancy for individual models which contain 400 000 ($\mathrm{1000}\times \mathrm{100}\times \mathrm{4}$) SSTAs. A box plot illustrating the likely and very likely range of the probability density distribution function of the model and proxy absolute discrepancy is shown for every model in Fig. 3b. The vertical lines inside the boxes are the medians of the probability density distribution function, which means that there is a probability of 50 % that the model and proxy discrepancy are less than this value. The box in Fig. 3b shows the range of the likely model and proxy discrepancy; the dashed error bar shows the range of the very likely model and proxy discrepancy. The likelihood definition adopted here is that from the IPCC Guidance Note (Mastrandrea et al., 2010); the term likely means the probability of the outcome ranges from 66 % to 100 %, and the term very likely means the probability of the outcome ranges from 90 % to 100 %. By using this method, we take uncertainties in the proxy data and model simulation into account. From Fig. 3b, we find that the E values between proxy and models range from 1.52 to 2.28 ^∘C. As the majority of the model SSTAs are warmer than the reconstructed SSTAs (Fig. 3a), the models showing relatively less regional warming also show a lower discrepancy value E with proxy SSTAs over the MC than the models that produce higher warming (Fig. S2). The PlioMIP2 models' MMM shows a discrepancy value E of 1.58 ^∘C, which is lower than that of most of the PlioMIP2 individual models.

In order to explore whether those models that perform well in simulating preindustrial SSTs over the MC are also good at simulating the mPWP SSTAs over the MC, the discrepancies in the models and observed SSTs on the one hand and reconstructed SSTAs on the other hand are shown in Fig. S3 for both experiments. In terms of simulating the climate of the MC, there is no clear correlation between model performance in preindustrial climate and in the mPWP climate, implying that models, which reproduce the preindustrial climate of the MC well, are not always good at simulating mPWP climate of the MC in accordance with the proxy reconstructions. However, we note that the uncertainty in the discrepancies between model and reconstructed SSTAs in the MC is significant due to limited sites for comparison and high uncertainties in the proxy reconstruction. From Fig. S3, models seem to fall into the following two groups: models where the discrepancy in the Eoi⁴⁰⁰ experiment appears decoupled from that in the E²⁸⁰ experiment, and models where regional biases are linked between the E²⁸⁰ and Eoi⁴⁰⁰ simulations (including CESM2, EC-EARTH3.3, HadGEM3, NorESM1-F and CESM1.2). Models in the second group also appear to produce relatively warm temperatures in the Eoi⁴⁰⁰ simulation (Fig. S4). In terms of simulating the global SSTAs, the comparison of the model and data discrepancies on a global scale between the E²⁸⁰ and Eoi⁴⁰⁰ simulations is shown in Fig. S5, which shows a better correlation than the regional SST simulation.

3.2 Mean climate features

In this section, we describe the modelled mean climatic characteristics of the mPWP relative to the preindustrial over the MC in terms of SST, wind stress, P−E and SOS.

The MMM shows a regional average warm SSTA of 1.88 ^∘C over the MC, the spatial pattern of which is as shown in Fig. 4a. The WPWP and the southern sector of the Indian Ocean in the MC the anomaly is less warm compared with other regions. SSTAs are comparably warm over the northern Indian Ocean in the MC region. Oceans surrounding the land masses of the MC also show apparently warm SSTAs, but this is not a robust signal because the land–sea masks of the models differ in the region of the MC. Models showing the greatest regional average warm SSTAs in the MC also show the largest global average warmings (Fig. S4). The linear regression coefficient between regional (the MC) and global annual SSTAs is ∼1.06, implying that the SSTAs over the MC change at a similar magnitude to that of global SSTAs. The SSTA spatial patterns over the MC from individual models are shown in Fig. S6.

Figure 5 shows the correlation of equilibrium climate sensitivity (ECS) and MC SSTAs from models. Regarding the SSTAs, models consistently show warmer SSTs over the MC in the mPWP relative to the preindustrial; regional average SSTAs range from 0.43 ^∘C produced by NorESM1-F to 3.33 ^∘C produced by HadGEM3. HadGEM3, EC-Earth3.3 and CESM2 all produce SSTAs exceeding 3 ^∘C, representing the three largest warm anomalies in the region among the PlioMIP2 models; three of the lowest SSTAs, which are all colder than 1 ^∘C, are from NorESM1-F, NorESM-L and MRI-CGCM2.3 and are also the models with the lowest-equilibrium climate sensitivity (ECS) among PlioMIP2 models (Haywood et al., 2020), as shown in Fig. 5. On long timescales, ECS is constrained by pattern effects i.e. radiative feedbacks change as SST patterns evolve (Ceppi and Gregory, 2017; Zhou et al., 2016; Armour, 2017; Gregory and Andrews, 2016; Dong et al., 2020; Andrews et al., 2015; Stevens et al., 2016; Zhou et al., 2021). For CMIP5 models, the WPWP is a key region, a warming in this region dominates the difference in the radiative feedback by enhancing Earth’s albedo (Zhou et al., 2021; Dong et al., 2020). In the study of Haywood et al. (2020), which adopted CMIP6 and non-CMIP6 models, there is a relationship between the ECS of the models and their global temperature anomalies (Eoi⁴⁰⁰ minus E²⁸⁰). With regard to warming of the WPWP in CMIP6 models, the MC contributes to the climate feedback differences, and besides the effect of the MC, there is a stronger sensitivity of extratropical clouds to surface warming that may also result in differences (Dong et al., 2020). As such, the relationship of ECS and SSTAs on the MC is not exactly linear. For example, GISS-E2-1-G shows an SSTA that is colder than 1 ^∘C, but its ECS is relatively high. CCSM4-UoT and CCSM4-Utr have a modest ECS but show a high global near-surface air temperature anomaly (Haywood et al., 2020) and a relatively high SSTA over the MC in mPWP. The best-fit line of the linear relationship between ECS and SSTAs on the MC is ECS = 2.43 ^∘C + 0.74 × SSTA at the 95 % confidence level with R²=0.56.

Figure 5 SST anomalies over the MC versus equilibrium climate sensitivity (ECS) from all the PlioMIP2 models.

The warmer MC SSTs not only influence the atmospheric circulation in the mPWP but also bring higher evaporation due to higher water-vapour-holding capacity of warm air (the Clausius–Clapeyron equation) and more atmospheric moisture, which contributes to the change in the simulated precipitation. As an essential part of the hydrological cycle, the freshwater budget over land and ocean (P−E) influences ocean salinity and shapes the climate of the MC. Figure 4 illustrates the MMM anomalies of P−E, zonal wind stress (TAUU) and SOS over the MC in the mPWP relative to the preindustrial.

For the MMM P−E anomalies, if we take the MC as a whole, then the P−E is more positive in the mPWP over the MC relative to the preindustrial (Fig. 4b); 13 out of 17 models and the MMM consistently show an increase in P−E (Fig. S7). Unlike SST, P−E does not change consistently over the MC but varies throughout the region (Fig. 4b), increasing over the western Pacific Ocean and the northeastern Indian Ocean, while decreasing over the western oceanic sector of the MC gateway. The difference in the distribution of land and sea in the Eoi⁴⁰⁰ and the E²⁸⁰ simulations is evident from hatching and black lines in Fig. 4b, which denote land grid boxes in the mPWP experiment and continental outlines of the preindustrial experiment, respectively. P−E increases over most land areas of the MC and also increases generally over grid cells which are ocean in the preindustrial but land in the mPWP, leading to a land bias of the precipitation increase. The spatial patterns of P−E anomalies from individual models are quite different (Fig. S7). In the case of some models, the simulated P−E anomaly is nearly opposite to each other, such as for COSMOS and EC-Earth3.3. Even the same family of models can show different results (e.g. CCSM-Utr and CCSM-UoT).

A large amount of precipitation provides energy to the atmosphere by releasing a large amount of latent heat, which fuels the atmospheric circulation. In order to explore the changes in the atmospheric circulation over the MC, we analyse the TAUU anomalies, as shown in Fig. 4c. Green shading indicates a westward anomaly. There are westward anomalies over the northwestern Pacific Ocean and eastward anomalies over the ocean to the west of Australia. This means that the tropical trade wind has been enhanced on the side of the northern WPWP during the mPWP. This signal can be seen from half of the available PlioMIP2 models (Fig. S8). There are extreme values over some regions which belonged to a land mass during the mPWP but have turned into ocean in modern times. The extreme values found there may therefore result from the change in the land–sea distribution that also impacts on the prevalence of surface roughness. The latter impacts atmospheric flow more over land than over ocean.

As a result of the changes in the hydrological system, SOS changes surround the MC. Figure 4d illustrates the spatial distribution of SOS anomalies. We find a decrease in the SOS in most of the ocean regions belonging to the MC. This pattern is especially pronounced in the northern part of the Indian Ocean, where SOS decreases by more than 1.2 PSU. For the Pacific Ocean sector and the southern Indian Ocean sector of the MC, SOS decreases by less than 1 PSU. Overall, we find that seawater becomes fresher at the surface around the MC in the mPWP in comparison to the preindustrial. This finding is consistent in 14 out of 16 PlioMIP2 models and in the MMM. Spatial patterns of SOS anomalies for individual models are shown for reference in Fig. S9.

3.3 Ocean currents

Temperature and density differences are factors that can contribute to changes in the ITF, which is an important process for regulating the climate of the MC. As the seawater characteristics of Pacific Ocean and Indian Ocean are so different, oceanic flow through the MC not only transports volume but also creates a net transport of salinity and heat between these two oceans. This process can directly influence both local and remote climate. Therefore, we quantify and investigate the changes in the oceanic flow through the MC in the mPWP. In this section, we address the following questions: is there more fresh and warm water injected into the Indian Ocean from the Pacific Ocean in the mPWP relative to the preindustrial? Where do the anomalies originate from? Do the anomalies result from the surface ocean or rather from processes and characteristics that are shaped in deeper layers of the sea?

We study the oceanic flow with the available data. During the calculation of ocean currents, we find that there are some disagreements between results derived from regridded data and analyses performed on the native ocean model grids of individual models. Therefore, to reduce errors and improve the accuracy, we use native grid data to calculate the oceanic flow intensity. In order to include all the seawater transport through this gate and retain consistency for all the models, we calculate the seawater volume transport through a broad transect, as shown with the red line in Fig. 6a. In the modelled mPWP, some regions, which are islands in the preindustrial, are still a peninsula that is connected to the Asian mainland in some models (CESM2, EC-Earth3.3, CESM1.2, MIROC4m, HadCM3, GISS-E2-1-G, CCSM4, CCSM4-Utr and CCSM-UoT). In these models, except for HadGEM3, which adopted the same land–sea mask of the MC in both experiments, the Malacca Strait and some other small straits in the northern part of the MC are closed, as shown in Fig. 1. Except for this gateway, the strait between New Guinea and Australia is also closed. Therefore, the Sunda and Sahul shelves caused some straits of the modern MC that are present in the E²⁸⁰ experiment to be absent in the Eoi⁴⁰⁰ experiment.

Figure 6 (a) The MMM zonal seawater volume transport anomaly in the Eoi⁴⁰⁰ experiment relative to the E²⁸⁰ experiment. Blue shading indicates a westward anomaly, and brown shading indicates an eastward anomaly. The red line is the transect to calculate the ocean volume transport intensity in panel (b). The hatches indicate land grid boxes in the mid-Pliocene-Eoi⁴⁰⁰ experiment; the black lines indicate the continental outlines of the preindustrial experiment. (b) The integrated ocean volume transport intensity through the gateway. Negative values indicate westward transport. Stippled areas show significant multi-model differences, according to a t test, at the 90 % confidence level. The maps of zonal seawater volume transport for individual models are as shown in Fig. S13.

The spatial pattern of the vertical integrated zonal seawater volume transport anomaly is shown in Fig. 6a. There is a westward transport anomaly along the northeastern coast of the New Guinea, which favours the ocean mass transport of the ITF. For the oceanic transport through the gateway, the westward water transport via the Timor passage has been strengthened over the north but weakened over the south in the mPWP compared to the preindustrial. Models and modern observation both show a strong flow through the gateway between MC and Australia. According to the observation between 2004 and 2006 (Sprintall et al., 2009), the ocean volume transport via the Indonesian Throughflow to the Indian Ocean is 15 Sv. By integrating the water volume transport through this gateway vertically, we derive the total seawater transport volume through the MC, which is shown in Fig. 6b. Based on the model results, the flow through the MC ranges from 10.9 Sv (CCSM4-Utr) to 23.7 Sv (MIROC4m) westward in the preindustrial and from 9.5 Sv (HadGEM3) to 29.5 Sv (NorESM-L) westward in the mPWP. The result of the E²⁸⁰ experiment from every model shows the same transport direction and a magnitude comparable to the observations. Compared to the preindustrial, the strength of the flow shows an increase in 10 out of 15 models, ranging from 5.6 % (CCSM4) to 43.8 % (IPSL-CM5A). Five reductions are observed in the simulation by CESM2, IPSL-CM6A-LR, COSMOS, MIROC4m and HadGEM3. As mentioned before, the gateway geometry in HadGEM3 is the same in both E²⁸⁰ and Eoi⁴⁰⁰ experiments, which may lead to a result that is different from those of the other models. Therefore, 10 out of 15 of the PlioMIP2 models show agreement on an increase in ocean volume transport through the MC in the mPWP in comparison to the preindustrial (NorESM1-F and MRI-CGCM2.3 are absent). Moreover, 5 of 15 models (CESM2, IPSL-CM6A-LR, COSMOS, MIROC4m and HadGEM3) disagree with the rest of the ensemble on the sign of transport change.

In order to explore which depth of the ocean determines the changes in the throughflow, we investigate the vertical structure of water volume transport through a broad passage of the ITF–Timor passage, as shown by the red line in Fig. 1, the MMM of which is shown in Fig. 7. In the MMM result, the strongest westward water transport is via the current located around 10^∘ S. In general, the direction of water transport is the same in mPWP and in the preindustrial. This encompasses the westward transport through the main ITF passage and the relatively weak eastward transport through the southern branch of the gateway simulated by some models (pale yellow shading in Fig. 7). The strong westward flow counterbalances the weak eastward flow, and there is the net transport of fresh and warm water from the Pacific Ocean to the Indian Ocean. In both experiments, the water mass is mainly transported in this gateway in the upper ocean above 1000 m depth, which is consistent based on the results of the MMM and individual models (Fig. S10).

Figure 7 The MMM profiles of the ocean volume transport through the MC in the mid-Piacenzian warm period (mPWP) (a), in the preindustrial (b) and the anomalies of the mPWP relative to the preindustrial (c). The graphs attached to the panels on the right-hand side show the horizontally averaged ocean volume transport. Values in parentheses indicate the total water volume transport through the gateway; positive values indicate westward transport. Warm shading indicates direction from west to east. Cold shading indicates direction from east to west. The upper colour bar is for absolute values in panels (a) and (b). The colour bar on the bottom is for anomaly values in panel (c). Stippled areas show significant multi-model differences, according to a t test, at the 90 % confidence level.

The profiles of the throughflow anomalies are shown in the bottom panel of Fig. 7. From the difference between Eoi⁴⁰⁰ experiment and E²⁸⁰ experiment, it can be seen that the changes in meridional integrated ocean volume transport through Timor passage occur largely above 1000 m depth, especially strongly at the surface and at a depth of around 800 m but without revolving the direction of water transport. In terms of the structure of the throughflow anomaly, there is a negative difference (stronger flow) in the northern part of the throughflow and a positive change (weaker flow) between 12 and 8^∘ S. This finding suggests that the flow through Timor passage developed closer to the northern boundary of the gateway in the mPWP than it is located during the preindustrial at the layers above 1000 m.

3.4 Clustering

So far we have studied the climatic features of the MC by using the multi-model ensemble method, which has often been employed when analysing results from multiple models. However, we have shown that, although models largely show agreement on a warmer and wetter climate over the MC in the mPWP, the spatial patterns of climate anomalies vary between models (e.g. Fig. S6). These differences might counterbalance each other to some extent, and therefore, potentially important climate signals from individual models will be eliminated when using the MMM. Moreover, some of the employed models belong to a group of general circulation models that share model components or may be similar in other respects. Consequently, similar signals from the same model family may be amplified in the ensemble of multi-model results. In order to investigate differences between models and to show the signals which might disappear in the MMM, we use clustering analyses in this section. Thereafter, we introduce a new metric – the multi-cluster mean (MCM).

In the beginning of this section, we demonstrate the method of classifying models into groups based on the spatial features of SSTAs and precipitation anomalies over the MC from individual models by applying pattern correlation and the hierarchical clustering method, which are derived, for example, from Sierra et al. (2017) and Knutti et al. (2013a). After averaging the results from multiple clusters, we obtain the MCM results. Then we illustrate the spatial anomaly pattern of every group and of the MCM and compare them with the MMM. After that, we investigate the performance of cluster groups and of the MCM at the example of the mPWP and the preindustrial simulation and investigate if these show less discrepancy with observation or reconstruction than the MMM.

Figures 8a and 9a illustrate the procedures of clustering. The first step is to find the two models with the smallest distance and join them into one group. For each subsequent step, the two closest groups or models are merged into a larger group. The distance between two models is the pattern correlation between the multi-annual mean SSTAs or precipitation anomalies of these two models. The x axes of Figs. 8a and 9a indicate the distance of every two joined groups at every step. The pattern correlation measures the similarity of the spatial patterns between two models (Taylor, 2001), which is the Pearson correlation coefficient of the same variable of two different maps from two models. Here, for SSTAs, we only use grid boxes which are ocean in all the models. We choose the average-linkage method (Wilks, 2011) to merge groups, which means that the distance between two groups is the averaged distance between individual members of the groups.

Figure 8 (a) Dendrogram of the SSTA clustering analysis based on the regional SSTAs over the MC from 17 PlioMIP2 models. Below the dendrogram, the clustering step is shown as a function of distance between groups. (b) The SSTA patterns of the composite for different clusters, multi-model mean (MMM) and the multi-cluster mean (MCM) after removing the regional mean SSTAs. Cluster 1 is CESM2; cluster 2 is MIROC4m and HadCM3; cluster 3 is GISS-E2-1-G, CCSM4-UoT, CCSM4-Utr, NorESM1-F, CCSM4, CESM1.2, EC-Earth3.3 and NorESM-L; cluster 4 is COSMOS, HadGEM3, IPSL-CM6A-LR, IPSL-CM5A and IPSL-CM5A2; and cluster 5 is MRI-CGCM2.3. White shading indicates the grid boxes that contain land for all the models of the cluster. The grid boxes outside the black lines are ocean grid boxes for all of the models of the cluster. Stippled areas show significant differences, according to a t test, at the 90 % confidence level.

Figure 9 The same as Fig. 8 but for precipitation. Cluster 1 is CESM2, CESM1.2, CCSM4 and CCSM4-UoT; cluster 2 is IPSL-CM6A-LR, IPSL-CM5A2, IPSL-CM5A, COSMOS, HadGEM3, MIROC4m and MRI-CGCM2.3; cluster 3 is NorESM-L and NorESM1-F; cluster 4 is HadCM3, EC-Earth3.3 and CCSM4-Utr; and cluster 5 is GISS-E2-1-G. The hatching and the black lines are the same as in Fig. 4. Stippled areas show significant differences, according to a t test, at the 90 % confidence level.

In the bottom panel of Fig. 8a, we find a large jump after clustering step 12, suggesting an appropriate stopping point with five groups, the group members of which are shown in Fig. 8a. CESM2 shows an SSTA spatial pattern that is different from that of the other PlioMIP2 models. Consequently, it forms cluster 1 on its own. The same is true for MRI-CGCM2.3, but this model could be merged into another cluster if we moved one step forward. Cluster 2 contains MIROC4m and HadCM3. Cluster 3 contains the largest number of models, with eight members, including almost all of the CCSM series model (except for CESM2), NorESM series models, GISS-E2-1-G and EC-Earth3.3. IPSL-CM5A2 is an updated version of IPSL-CM5A in terms of updated components and technical characteristics (Sepulchre et al., 2020). This is clearly shown in the clustering results, with IPSL-CM5A and IPSL-CM5A2 being the two closest models. All the three of the IPSL-CM models have been classified into cluster 4, together with COSMOS and HadGEM3. Figure 8b illustrates the ensemble SSTA pattern of every cluster after removing the regional mean SSTA. Clusters 4 and 5 show a relatively warmer pattern for the WPWP, which is opposite to the results present in cluster 3. For the Indian Ocean, cluster 4 shows a relatively weak warm pattern, which means the Indian Ocean warms to a lesser extent than the region, especially the WPWP. Spatial distributions of SSTAs in the Indian Ocean are similar for clusters 3 and 2, with temperatures that are warmer in the north but less warm in the south, which is the opposite of the results from cluster 5. Since models are equally weighted in the ensemble, and cluster 3 contains the largest number of PlioMIP2 models, it dominates the MMM SSTA pattern (Fig. 8b). The MCM is the mean of these five clusters; in other words, it is the multi-model mean with weights. The MCM can avoid the signals being overweighted from the same model family; it conserves details of different types of climate realisation that may come from similar or very different models. For example, the centre of the WPWP shows a warmer anomaly, which is a result consistent for all clusters, except for cluster 3. A comparison between the warming of the WPWP and the eastern Indian Ocean on different clusters indicates three possible SSTA patterns, namely that the eastern Indian Ocean warms more than the WPWP (cluster 2, cluster 3 and the MMM); that the eastern Indian Ocean warms less than the WPWP (cluster 4); and that there is no apparent difference in warming between these two regions (cluster 1, cluster 5 and the MCM).

In addition to the cluster based on SSTAs, we also sort models based on the precipitation anomaly spatial patterns, as shown in Fig. 9. Similar to the SSTA analyses, we find a jump after step 12. Therefore, we classify the models into five groups (Fig. 9a). Four of the CCSM series models (CESM2, CESM1.2, CCSM4 and CCSM-UoT) form cluster 1. The IPSL-CM-type models are members of the cluster 2, together with COSMOS, HadGEM3, MIROC4m and MRI-CGCM2.3. This cluster is the biggest one, with seven models. Similar to clustering based on SSTAs, IPSL-CM5A and IPSL-CM5A2 are still the two most similar models; they are also the closest group for IPSL-CM6A-LR. NorESM-L and NorESM1-F form cluster 3. The remaining model in the CCSM series, CCSM4-Utr, joins EC-Earth3.3 and HadCM3 as cluster 4. For precipitation, cluster 5 (containing the single-model GISS-E2-1-G) shows a decrease in precipitation over the centre of the MC. For the WPWP, cluster 1 and cluster 2 show relatively wetter anomalies; cluster 3 and cluster 4 show basically drier anomalies with spatial heterogeneity, whereas cluster 5 (GISS-E2-1-G) alone shows a relatively drier climate in the northwest but wetter conditions in the rest of the WPWP. All clusters except cluster 5 show a relatively dry anomaly over the north of the WPWP. Similar to SSTA-based clustering, the cluster with the largest number of models dominates the spatial pattern of the MMM and leads to some signals in other clusters disappearing. Overall, compared with the eastern Indian Ocean, the WPWP is relatively wetter (cluster 1, cluster 2, the MMM and the MCM), drier (cluster 3 and cluster 4) or receives similar amounts of precipitation (cluster 5). If we compare Fig. 8b with Fig. 9b, the pattern of the MCM is more similar to the pattern of the MMM in precipitation than in the SSTAs.

Figure 10 compares the five clusters with the reconstructed SSTA data. In order to quantify the performance of the model clusters when simulating the mPWP SSTAs, the SSTA-based clusters are validated in the same way as individual models and as the MMM, as shown in Figs. 3 and S3. The performance of the MCM is similar to the MMM for every site. However, in contrast to the MMM, the MCM can take the model weight into account to avoid the duplication of model signals from the same family and also retain the climate signals we might lose in the multi-model ensemble. An example is the warm anomaly signal in the centre of the western Pacific Ocean.

Figure 10 Comparison of reconstructed and modelled SSTAs over the MC in the mid-Piacenzian warm period (Eoi⁴⁰⁰) relative to the preindustrial (E²⁸⁰) from the SSTA-based clusters and the multi-cluster mean (MCM). Cluster 1 is CESM2; cluster 2 is MIROC4m and HadCM3; cluster 3 is GISS-E2-1-G, CCSM4-UoT, CCSM4-Utr, NorESM1-F, CCSM4, CESM1.2, EC-Earth3.3 and NorESM-L; cluster 4 is COSMOS, HadGEM3, IPSL-CM6A-LR, IPSL-CM5A and IPSL-CM5A2; and cluster 5 is MRI-CGCM2.3. Units are in degrees Celsius. The MMM SSTAs shown in Fig. 4a. Note that in this figure, the regional average SSTA has not been removed as in Fig. 8b. White shading and black lines are the same as in Fig. 8b.

4 Discussion

4.1 Physical mechanisms

The MC is a region of intense atmospheric–oceanic coupling and exchange, where slight changes can affect climate afar (Gordon, 2005b; Yoneyama and Zhang, 2020). Due to differences in land distribution, topography and CO₂ concentration, we find that the MC experienced warmer SSTs, a higher freshwater budget, lower SOS, strengthened surface zonal wind in WPWP and increased oceanic volume transport through the mPWP compared to the preindustrial.

In the mPWP, wind and large-scale atmospheric circulation change in response to an overall warmer climate. According to the PlioMIP2 simulations, the Northern Hemisphere warmed more than the Southern Hemisphere (Haywood et al., 2020). As a result of the variation in the interhemispheric meridional heat fluxes and energetic constraints for the Intertropical Convergence Zone (ITCZ) position, the Southern Hemisphere subtropical high pressure system is intensified, and the western Pacific easterly trade wind system is enhanced correspondingly (Pontes et al., 2022), which is consistent with the strengthened zonal wind stress that we find in Fig. 4c. In terms of atmospheric circulation, the Walker Circulation consistently shows a westward shift within PlioMIP2 models, but the change in the strength proved to be highly model dependent (Han et al., 2021). In contrast, PlioMIP1 models suggested a weakened Walker Circulation in the mPWP in relative to the preindustrial (Corvec and Fletcher, 2017).

On multi-decadal timescales, the Pacific trade wind system plays a role in the circulation strength and patterns within the MC and Indian oceans (Sprintall et al., 2014; Vecchi and Soden, 2007). This has been observed in the 1970s, when the weakening of the trade winds resulted in the western Pacific shoaling thermocline anomalies (Vecchi et al., 2006; Wainwright et al., 2008), and this anomaly can be transmitted to the Indian Ocean by the ITF (Alory et al., 2007). Results from coupled models also reveal that the strength of the ITF reduced correspondingly with the weakened Pacific trade winds (Sprintall et al., 2014). On the projected anthropogenic warming world, most models suggest a reduction in the equatorial trade winds and a decreased ITF transport (Sen Gupta et al., 2012). In our study, most models show an increased ITF transport and also suggest an intensified western equatorial Pacific zonal wind stress (Fig. S11). The wind can influence the water properties, thermocline and sea level through wave processes and therefore impact the formation of the WPWP (Wijffels and Meyers, 2004; Drushka et al., 2010; Pujiana et al., 2013; Sprintall et al., 2014). Changes in the WPWP have a direct effect on the ITF, as the ITF is driven by the sea level gradient between the western Pacific (high sea surface height) and the eastern Indian Ocean (low sea surface height; Wyrtki, 1987; Sun and Thompson, 2020; Gordon, 2005b). On decadal timescales, the Pacific Ocean sea surface height (SSH) variability controls more than 85 % of the ITF variation, according to ocean reanalysis data (1959–2015; Shilimkar et al., 2022). Here, we also correlated the annual ITF volume transport with annual SSH by using data available from four models (Fig. S12). We find that the western Pacific SSH is generally positively correlated with the ITF. Furthermore, the SSH gradient between the western Pacific and the eastern Indian Ocean correlates well with the ITF in both the E²⁸⁰ and Eoi⁴⁰⁰ experiments of GISS-E2-1-G and the Eoi⁴⁰⁰ experiments of CESM2 and IPSL-CM6A-LR, but it has insignificant correlation with the ITF in Eoi⁴⁰⁰ experiment of EC-Earth3.3, as shown by values at the top right of Fig. S12. Moreover, the regions in which SSH correlates with the ITF differ between experiments and models. For CESM2, for example, experiments Eoi⁴⁰⁰ and E²⁸⁰ show different spatial patterns of the correlation coefficient map. However, models that show increased ITF in the Eoi⁴⁰⁰ do not always show an increase in the SSH gradient, and vice versa, implying that, at least for these four models, changes in the SSH gradient are not robust drivers for the changes in the ITF between the mPWP and the preindustrial.

The ITF can transfer the signal of ENSO to the Indian Ocean and further afield. As shown in a study of ENSO from PlioMIP2, models GISS-E2-1-G, COSMOS, IPSL-CM6A and CESM2 show an El-Niño-like mean state, while CCSM4-Utr, EC-Earth3.3 and CCSM4 show a “La-Niña-like” mean state in the mPWP compared to the modern Pacific Ocean (Oldeman et al., 2021). All the models that show a clear La-Niña-like mean state also show a stronger Pacific Walker Circulation (Fig. 10a from Han et al., 2021) and corresponding intensified ITF; IPSL-CM6A, CESM2, COSMOS and HadGEM3 show El-Niño-like mean states at weakened ITF. Although GISS-E2-1-G shows an El-Niño-like mean state and unspecific zonal wind stress changes over the western tropical Pacific Ocean, it produces a stronger ITF. The rest of the models (IPSL-CM5A, CESM1.2, HadCM3 and NorESM-L), with intensified ITF, and MIROC4m, with slightly weakened ITF, show a slight La-Niña-like mean state. IPSL-CM5A2 and CCSM-UoT, with intensified ITF, show a slight El-Niño-like mean state (Oldeman et al., 2021; Fig. 4c).

The ITF as a component of the global ocean circulation system also interacts with other circulations. In PlioMIP2, most models simulate a stronger Atlantic Meridional Overturning Circulation (AMOC) in the mPWP than the preindustrial (Zhang et al., 2021b; Weiffenbach et al., 2023). On the centennial scale, the ITF decreases if the AMOC slows down, as a northward surface flux anomaly implies a strengthened convergence of ocean water masses and increases SSH of the Indian Ocean (Sun and Thompson, 2020). Yet, as discussed for SSH above, the connection between changes in AMOC and ITF remains to be studied. Nevertheless, in the study of a simulated future CO₂-induced warming world, changes in the overturning circulation in the Pacific basin explain the change in the ITF (Sen Gupta et al., 2016). The northward flow of deep waters entering the Pacific basin at 40^∘ S is weakened, and the associated net basin-wide upwelling to the north of the southern tip of Australia decreases, ultimately reducing the ITF (Fig. S10 in Sen Gupta et al., 2016). In the mPWP, the Bering Strait is closed, and water cannot be transported through this gateway and thus contribute to the increased ITF volume transport. In contrast, the volume transport through the South Pacific basin and the contribution to the ITF can be checked (Fig. S16). We find that all models that show decreased ITF volume transport show decreased volume transport through the South Pacific basin (shallow or deep) in the mPWP relative to the preindustrial. For models showing increased ITF volume transport, the closed Bering Strait plays a role as the water leaves the Pacific Ocean decreased, except for the HadCM3, which shows insignificant water transport via Bering Strait in the preindustrial. All models show an increased ITF volume transport, but the GISS-E2-1-G model shows a positive anomaly of the volume transport through the South Pacific basin above 2000 m in the mPWP relative to the preindustrial, and it can contribute to the increased ITF volume transport.

The ITF is both impacted by and has an impact on the state of the surrounding oceans (Auer et al., 2019; Brierley and Fedorov, 2016; Gordon, 2005b). By comparing the simulations with the MC gateways that are closed and opened, the removal of the ITF causes warmer SSTs in the tropical Pacific and cooler SSTs in the southern Indian Ocean (Lee et al., 2002). Moreover, simulations in which the ITF is blocked also show an eastward shift in the atmospheric deep convection and precipitation in response to the warmer Pacific SSTs (Sprintall et al., 2014). In the mPWP, the substantially increased western Pacific SST results in higher atmospheric moisture content (Allen and Ingram, 2002), which provides substantial water vapour that contributes to increased precipitation over the MC as part of the “wet-gets-wetter, dry-gets-drier” response. As a consequence of multiple factors, the spatial precipitation anomaly pattern forms, as shown in Fig. 4b (Feng et al., 2022).

Precipitation minus evaporation is the largest freshwater component entering the ocean (Byrne et al., 2016). Increased fresh water flux provides a way to explain the general decrease in the SOS over the MC, especially the decreased SOS in the northern Indian Ocean. The changes in the topography also play an important role. The closure of the gateway between New Guinea and Australia can change the water source from the saline South Pacific to the less saline North Pacific (Cane and Molnar, 2001), which can cause the northern Indian Ocean SOS to decrease significantly. Furthermore, higher precipitation and wetter climate over the Sunda shelves, India, the south slope of the Tibetan Plateau and the surrounding regions (as shown in Fig. 5b of Haywood et al., 2020, and Fig. 1b of Feng et al., 2022) may increase the local run-off, which can also inject more fresh water into the ocean and then dilute salinity.

Large-scale flow patterns seem to be simulated unambiguously across the PlioMIP2 model ensemble. In particular, 10 out of 15 models agree on increased westward ocean volume transport through the MC during the mPWP, and 5 out of 15 models show decreased volume transport.

4.2 Summary of the PlioMIP2 model performance

While MC SST are significantly warmer in the mPWP, it seems that within the PlioMIP2 model ensemble, the MC climate is to some extent decoupled from model ECS (Fig. 5). Nevertheless, we find that models with a large ECS also produce rather warm mPWP conditions around the MC, and models with a small ECS show relatively cooler conditions (Fig. 5). When using models to predict the climate of the MC for different levels of CO₂ or changes in Earth surface conditions, it is important to note that the models that best reproduce a recent regional climate of the MC, as quantified by a comparison to modern observations, are not necessarily the same models as those that produce the best agreement between simulated and reconstructed mPWP SSTAs (Fig. S3). Yet, at global scale, there is a correlation between the model performance on simulating mPWP and modern climate (Fig. S5). From that point of view, the mPWP proves to be a real test bed for the skill of models for climate states that strongly differ from the recent one.

In general, we find agreement between proxy reconstructions of SSTA and mPWP climate simulations. At the ODP214 site, a small number of models fall out of the uncertainty range of the reconstruction at the warm end. Yet, we note that ODP214 indeed appears somewhat cool in comparison to nearby SSTA reconstructions in the region. This leaves us with the possibility that the models are not able to capture climate mechanisms that shaped the mPWP climate recorded at the ODP214 site. An alternative explanation for the discord at this site may be sought within the proxy data or in the interpretation of the data. We note that the model–data comparison is hampered by large uncertainty ranges of reconstructed SSTAs and by the limited availability of suitable SSTA reconstructions in the region. This highlights the value of more regional MC SSTA reconstructions for the mPWP. We also find that the simulated climate in PlioMIP2 has a much smaller uncertainty range than the reconstructed climate, as seen by the spread of the simulated SSTAs across the PlioMIP2 model ensemble, which is much smaller than the uncertainty bars attributed to reconstructions (Fig. 3a). This suggests that the MC proxy reconstructions are shaped by additional aspects that the models are not responsive to.

4.3 Model hierarchy and clustering

The multi-model ensemble mean (MMM) has frequently been adopted to illustrate modelling results. We found the MMM produces smaller discrepancy than more than half of the PlioMIP2 individual models in both mPWP and preindustrial climate simulations (Fig. S3). In the case of model independence, the model biases will be partially cancelled, which results in the MMM outperforming individual models (Solomon et al., 2007; Knutti et al., 2010). Yet, some models, such as models adopted in Coupled Model Intercomparison Project phase 6 (CMIP6), are not fully independent from each other and might have similar biases (Abramowitz et al., 2019). Boé (2018) quantified the relationship between the CMIP5 model component replication and the proximity of their results and found that there is a clear relationship at both global and regional scales. Model clustering is a way to study proximity of the spatial patterns of model results and to weight models to reduce the reduplication of similar biases from the same model family.

In PlioMIP2, as in CMIP6, not all models are independent from each other. There are models that have been developed from the same source, with shared codes or using different generations of the same components, as shown in Table 1. A study from Ho et al. (2016) adopted outputs from 41 CMIP5 models, clustered the models according to the concept of “model genealogy” (Knutti et al., 2013b) and found the difference between the weighted and unweighted multi-model ensemble mean. Yet, general conclusions from both methods did not change. In this study, we cluster models into groups according to the proximity of their results when simulating SSTAs and precipitation (Sect. 3.4). Although some models used the same ocean component, they turned out to be in different SSTA-based cluster groups, such as CESM2. Descendants of CCSM4, CCSM4-Utr, CCSM4-UoT, CESM1.2, CESM2, NorESM-L and NorESM1-F are all in cluster 3, except for CESM2 (Fig. 8a). Compared to the other CCSM4 models, CESM2 has updated the POP ocean component not only with respect to the parameterisation but also with regard to the use of new schemes (Danabasoglu et al., 2020). Moreover, CCSM4, CESM1.2 and CESM2 have very different atmospheric components. In the study of Pacific Walker Circulation, CCSM4 showed a strengthened Walker Circulation, while both CESM1.2 and CESM2 showed a weakened Walker Circulation (Feng et al., 2020b). Nearly all of the models that include the NEMO ocean component (EC-Earth3.3, IPSL-CM5A, IPSL-CM5A2, IPSL-CM6A-LR and HadGEM3), except for EC-Earth3.3, appear in the SSTA-based cluster 4 (Fig. 8a). In the analysis of El Niño variability from PlioMIP2 (Oldeman et al., 2021), the La-Niña-like mean state model group and the El-Niño-like mean state group are not the same as the SSTA-based clusters found in our study. This implies that clustering depends on the region and on the studied physical quantity. In terms of the atmospheric component, some models used the same atmosphere component but have not been clustered into the same group based on precipitation. For example, CCSM4, CCSM4-UoT, CESM1.2 and CESM2 are all in cluster 1, but CCSM4-Utr is not in this group. The two closest models in both SSTA-based and precipitation-based clusters are IPSL-CM5A and IPSL-CM5A2. They adopted the same atmospheric component, oceanic component and vertical and horizontal resolutions. Compared with the IPSL-CM5A, IPSL-CM5A2 updated model components and retuned the cloud radiative forcing (Tan et al., 2020), but both models still produce a similar climate signal over the MC. From the discussion above, we conclude that the model configuration can provide a reference to cluster models, but even models of the same model family may still produce different climatic signals, depending on the analysis region or the studied climate characteristic. These results suggest that there is not always a clear trace of model family apparent via the clustering.

Moreover, the multi-model mean will lead to the loss of climatic signal (Knutti et al., 2010). In our study, there are noticeable variations between the spatial anomaly patterns from individual models. By illustrating the model results in cluster groups, climatic signal loss can be reduced (e.g. the warm anomaly signal in the centre of the western Pacific Ocean; Fig. 8b).

5 Conclusions

The subaerial Sunda and Sahul shelves and relatively high atmospheric CO₂ concentration, combined with other forcings, led to a different climate of the MC in the mPWP than in the preindustrial. In line with the global climate, all the 17 PlioMIP2 models (including HadGEM3) show higher SSTs in the mPWP MC, ranging from +0.43 to +3.33 ^∘C, with a MMM warming by 1.88 ^∘C; 13 out of 17 models show fresh water flux P−E increased in the MC, with the anomaly ranging from +0.04 to +0.50 mm d⁻¹. In terms of atmospheric circulation, the easterlies over the western tropical Pacific Ocean show enhanced intensity in the MMM; there is an anomalous eastward wind stress over the ocean near western Australia. Seawater salinity affects the thermohaline circulation; 14 out of 16 models (MRI-CGCM2.3 is absent here) show a lower SOS in the MC.

Even though the topography of the MC in the mPWP acts as a barrier for volume transport between the Pacific Ocean and the Indian Ocean compared to the preindustrial, 10 out of 15 models show the ITF intensity increasing by 5.6 % to 43.8 %. CESM2, IPSL-CM6A-LR, COSMOS, MIROCm4 and HadGEM3 show a decrease in intensity of the ITF, and the amplitude of westward flow is here lowered by 3.4 % to 39.4 %. In terms of spatial differences, the MMM shows that the ocean volume transport of the Timor passage is stronger in the north but weaker in the south in the mPWP relative to the preindustrial. Moreover, the water mass is mainly transported through the ITF in the ocean above 1000 m depth in both the preindustrial and the mPWP. The changes in the volume transport of ITF also mainly occur in the ocean above 1000 m depth.

The model performance has also been quantified in this work by comparing the discrepancies between model results and both reanalysis data and reconstructed proxy data. The results imply that models which reproduce modern climate well are not always good at simulating the mPWP climate of the MC, according to the proxy reconstructions. Although, at global scale, there is a correlation between model performance in terms of simulating the mPWP climate and modern climate (Fig. S6).

The comparison between individual models and the MMM suggests that the MMM could reproduce the preindustrial SST of the reanalysis better than 13 out of 17 models and show less of a discrepancy with reconstructed SSTAs than 12 out of 17 models.

There is a big difference in the spatial distribution of anomalies between models. Because some models will amplify similar signals and because some climatic signals from individual models will get lost in the MMM, a new metric, MCM, has been introduced. It takes the weight of the model families into account, provides a new perspective on the ensemble of multiple model results and demonstrates different climate states by means of clusters. The evaluation of the MCM shows that it outperforms 12 out of 17 models in simulating preindustrial climate and shows a discrepancy that is less than or equal to the reconstructed mPWP SSTAs for 10 out of 17 models, which suggests that the MCM is a valuable method for exploring multi-model results.