Low-frequency variability of the Pacific Subtropical Cells as reproduced by coupled models and ocean reanalyses

Low-frequency variability of the Pacific Subtropical Cells (STCs) is investigated using outputs from several models included in the two latest phases of Coupled Model Intercomparison Project (CMIP), CMIP5 and CMIP6, as well as ocean reanalysis products. Our analysis focuses on historical simulations and an idealised future scenario integration. Mass and heat transport diagnostics are employed to assess how coupled models and ocean reanalyses reproduce Pacific STCs total and interior transport convergence at the equator and their relationship with equatorial Pacific sea surface temperature (SST). Trends of mass and heat transport are also evaluated, in order to study how the STCs are expected to change in a warming climate. A large spread is obtained across models in simulated mass transports, confirming that coupled models do not agree on reproducing observed Pacific STCs dynamics, with very limited improvement by CMIP6 models. Compared to ocean reanalysis products, coupled models tend to underestimate the STCs interior transport convergence, and are less efficient on propagating the signal generated by the subtropical wind stress towards the equator. Also, mass transport obtained from ocean reanalyses exhibit larger variability, and these products also better reproduce the STCs-SST relationship. Future scenario simulations suggest a weakening (strengthening) of the heat transport by the North (South) Pacific cell under warmer conditions, with a general agreement across models. Equatorward mass transport trends do not confirm this for total and interior components, but they do for the western boundary component.


Introduction
The Pacific Subtropical Cells (STCs) are shallow meridional overturning circulation structures, connecting the equatorial ocean to the subtropical regions in the Pacific, Atlantic, and Indian Oceans. In particular, the Pacific STCs are extensively studied in literature: theorised in the 1990s (McCreary and Lu 1994;Liu 1994;Lu et al. 1998), their properties have been evaluated through observations Zhang 2002, 2004;Zhang and McPhaden 2006) and model simulations [see Section 1.1 of Graffino (2019) for a review]. Other than STCs, the shallow meridional overturning circulation is also defined by the Tropical Cells (TCs; Lu et al. 1998;Molinari et al. 2003), localised closer to the equator, but forced by different processes.
The structure of the Pacific STCs was first depicted by McCreary and Lu (1994), with a 2 1/2-layer model forced by an idealised wind stress pattern. The circulation depicted by McCreary and Lu (1994) was finally able to reconcile theoretical works (Luyten et al. 1983), observations (Fine et al. 1981(Fine et al. , 1987, and modelling studies (McCreary and Yu 1992) about the need of a connection between tropics and subtropics in the Pacific Ocean. Following studies by Liu (1994) and Liu et al. (1994) identified the deep branches of the STCs, helping to draw a schematic of the circulation. On a time-average perspective, the Pacific STCs are seen as a pair of overturning cells on both sides of the equator, with a subduction branch located at the subtropics, an equatorward (and mostly westward) sub-surface flow, a rising branch along the equatorial thermocline, and a returning flow toward the poles at the surface (Schott et al. 2004;Capotondi et al. 2005).
The equatorward mass transport is split into two contributions: western boundary (WB) and interior (INT). There is an hemispheric asymmetry in the pathways followed by water parcels in the tropical Pacific (Fine et al. 1981(Fine et al. , 1987. This arises from the presence of a high potential vorticity ridge, located at 9 • N (Lu and McCreary 1995). As a consequence, in the Northern Hemisphere the pathway followed by a water parcel goes more westward than in the Southern Hemisphere (Johnson and McPhaden 1999). This causes the WB component to be larger than the INT contribution in the Northern Hemisphere, and also makes the northern STC mass convergence smaller than the southern STC mass convergence (Lu et al. 1998).
The INT flow occurs roughly between 180 • W and 140 • W in the Northern Hemisphere, and between 160 • W and 90 • W in the Southern Hemisphere (Capotondi et al. 2005;Zhang and McPhaden 2006). Despite being smaller, it has been shown that the INT component dominates the Pacific STCs variability on interannual to decadal timescales (Lee and Fukumori 2003;Cheng et al. 2007;Hong et al. 2014), although it is partly compensated by the WB component (Capotondi et al. 2005;Cheng et al. 2007;Lübbecke et al. 2008). Furthermore, interannual variability of the INT component is able to affect ENSO variability (Huang and Wang 2001;Schott et al. 2008;Zilberman et al. 2013). For these reasons, in this study we mostly focus on the INT component. Klinger et al. (2002) and Graffino et al. (2019) focused on the relationship between the subtropical wind stress, the main STCs forcing mechanisms according to McCreary and Lu (1994), and the equatorial sea surface temperature (SST). This teleconnection process was introduced by Kleeman et al. (1999) as the v ′ T mechanism: thermal anomalies can be driven at the equator by changing the strength of the STCs circulation. A different hypothesis, called the v T ′ mechanism (Gu and Philander 1997), is related to the advection of thermal anomalies from the subtropics to the equator via the STCs. Despite being corroborated by Zhang et al. (1997), several observational studies questioned this hypothesis, suggesting that such thermal signal would shortly decay while propagating from the subtropics to the equator (Schneider et al. 1999;Pierce et al. 2000;Hazeleger et al. 2001).
By employing both uncoupled atmospheric and oceanic simulations, Farneti et al. (2014b) schematised a coupled interaction between the tropics and the subtropics involving the STCs, inferring an important influence of such overturning structures on the Pacific Decadal Oscillation (PDO; Mantua et al. 1997) and on the tropical Pacific Ocean state. Despite Nonaka et al. (2002) and McGregor et al. (2008) depicted a secondary role for STCs on driving SST decadal variability in the equatorial Pacific Ocean, their relevance was instead highlighted by several modelling studies (e.g., Solomon et al. 2003;Lohmann and Latif 2005;Wu et al. 2007;Lübbecke et al. 2008;Farneti et al. 2014a;Hong et al. 2014;Graffino et al. 2019) employing both ocean-only and coupled ocean-atmosphere simulations.
With respect to STCs decadal variability, Zhang (2002, 2004) observed a decreasing transport convergence at the equator from the 1970s to the 1990s, and then an increasing convergence during the early 2000s. A similar behaviour was reproduced in many studies (Nonaka et al. 2002;Lee and Fukumori 2003;Capotondi et al. 2005), and STCs variations at decadal timescale are confirmed by ocean-only (Cheng et al. 2007;Farneti et al. 2014a) and coupled ocean-atmosphere simulations (Lohmann and Latif 2005;Zhang and McPhaden 2006), as well as ocean reanalysis products (Schott et al. 2008). In particular, Farneti et al. (2014a) showed that an ocean model forced with atmospheric reanalysis is able to reproduce most of the observed STCs and equatorial SST variability in the Pacific Ocean during the twentieth century.
On the other hand, coupled ocean-atmosphere models struggle on reproducing observed features of the observed Pacific STCs variability in historical simulations, as well as its relationship with equatorial Pacific SST. In fact, Merryfield and Boer (2005) and Lohmann and Latif (2005) showed that coupled models reproduce the relationship between STCs interior transport changes and equatorial SST with different magnitudes, while Zhang and McPhaden (2006) proved that the anti-correlation between INT and WB components can vary significantly in coupled simulations. Comparing their results with Solomon and Zhang (2006), Zhang and McPhaden (2006) also suggested that an underestimated STCs-SST relationship in coupled models may arise from a reduced variability of subtropical wind stress in terms of suppressed Ekman pumping. In particular, Solomon and Zhang (2006) showed that salinity patterns over the 1025 kg m −3 density surface are strikingly different in the equatorial regions between coupled and ocean-only simulations, with the latter more faithful to observations; again, they ascribed this to an inadequate reproduction of the remote forcing of the equatorial thermocline via the STCs.
More recently, Yang et al. (2014) compared a subset of CMIP5 historical simulations with the SODAsi.1 ocean reanalysis product over the 1871-2008 period. Yang et al. (2014) found that SODA provides a cooling (warming) subsurface trend in the eastern/central (western) equatorial Pacific Ocean, along with an increased STCs convergence transport in both interior and western boundary components. On the other hand, most CMIP5 models showed a general subsurface warming trend over the whole Pacific basin, along with a weakened subtropical trade wind circulation and a reduced STCs mass transport convergence. This was confirmed by the analysis performed by Farneti (2017) on a subset of CMIP5 historical runs, also showing a weak correlation between STCs interior convergence transport and equatorial Pacific SST. Farneti (2017) attributed these discrepancies to an underestimation of the remote forcing provided by off-equatorial wind stress on the STCs.
Due to the large number of processes at play, which might impact STCs in the future (England et al. 2020), and the inherent uncertainties in climate projections (see reflections included in Knutti and Sedláček (2013)), trying to infer future Pacific STCs changes is particularly challenging. Moreover, in many cases climate projections do not even agree on the overall response of the tropical Pacific Ocean under global warming conditions (Knutson and Manabe 1998;Timmermann et al. 1999;Collins 2005;Vecchi et al. 2008;Park et al. 2009).
There are a few notable exceptions though. Among those, Merryfield and Boer (2005) and Lohmann and Latif (2005) analysed coupled simulations under global warming conditions until 2100, both finding a weakening of the STCs mass transport convergence at the equator. In particular, Lohmann and Latif (2005, Fig . 13) showed the different behaviour of the cells in the two hemispheres, with a weakening (strengthening) of the northern (southern) STC. This was confirmed by Park et al. (2009) analysing coupled model experiments with linearly-increasing CO 2 concentration, and by Wang and Cane (2011) using a subset of CMIP3 simulations of the twenty-first century. Furthermore, this occurs along with a weakening (strengthening) of the surface wind stress in the subtropical North (South) Pacific (Luo et al. 2009;Park et al. 2009;Wang and Cane 2011). Luo et al. (2009) studied STCrelated changes of the tropical thermocline in the twentysecond century, and found a robust increase (decrease) of the WB (INT) contribution across analysed models, but with no significant change of the total STCs transport. A robust weakening of the interior transport convergence was also found by Wang and Cane (2011), and it is explained in terms of barotropic and baroclinic adjustment of the tropical Pacific Ocean to the changing subtropical wind stress forcing. However, unlike Luo et al. (2009), Wang and Cane (2011) found no indication of a compensation by the western boundary contribution, leading to an overall decrease of the STCs total convergence. Comparing coupled model simulations in pre-industrial and twenty-first century globalwarming conditions, He et al. (2019) showed that the overall poleward heat transport by the STCs circulation in the Indo-Pacific Ocean is expected to decrease.
Realising the importance of Pacific STCs in shaping both local and global climate, our aim is to revisit their behaviour in state-of-the-art climate models and ocean reanalyses. Therefore, one of the purposes of this paper is to assess the relationship between the Pacific STCs strength and the equatorial Pacific SST at decadal timescales in CMIP5 and CMIP6 models simulations. In particular, we are interested in the role of STCs mass transport convergence, that is the sum of the equatorward components from both hemispheres, in driving the low-frequency equatorial Pacific SST variability, especially at decadal timescale. To achieve this, we compare STCs mass and heat transports, and their relationship with equatorial Pacific SST, as reproduced during historical conditions and in a future global warming scenario. Data and methods employed in the present study are described in Sect. 2. Low-frequency variability of the STCs mass transport and its relationship with equatorial Pacific SST are explored in Sect. 3.1, whereas Sect. 3.2 focuses on expected changes of the STCs circulation under an idealised future scenario. Main results, discussed and compared with previous works, and conclusions are provided in Section 4.

Data and methods
This section includes details about all coupled models and ocean reanalyses employed in the present study. A list of models and products, along with their references, is given in Table 1. All data are converted into annual means.

Coupled models
Coupled models are sophisticated tools, combining different numerical components to simulate the Earth system. They are essential for studying many aspects of climate science but, despite their increased complexity over the last decade, they still show deficiencies and biases with respect to observations. The latest phases of the Coupled Model Intercomparison Project, CMIP5 (Taylor et al. 2012) and CMIP6 (Eyring et al. 2016), include a wider range of models and experiments with respect to the previous phases, providing to the climate science community an unprecedentedly large dataset. All CMIP5 and CMIP6 simulations are available through the Earth System Grid Federation (ESGF) portal. All data were remapped on a rectilinear grid, keeping the original horizontal and vertical resolution unchanged unless otherwise specified. We used pre-industrial control (piControl), historical, and 1pctCO2 experiments.
Only one ensemble member was considered for each model. While we acknowledge that the use of many ensemble members might give more insights about model performances, we also think this choice does not undermine the main outcomes of this study. Furthermore, earlier studies about STCs in coupled models only used one ensemble member for most or all employed models (Lohmann and Latif 2005;Zhang and McPhaden 2006;Wang and Cane 2011;Farneti 2017), without diminishing the relevance of their results.
PiControl simulations are obtained by prescribing the CO 2 atmospheric concentration to pre-industrial values, mostly for studying unforced climate variability. They should represent the quasi-equilibrium state of the climate system; in reality, since the deep ocean would take thousands of years to adjust to the imposed forcing, shorter runs are often employed (Eyring et al. 2016). For particular purposes, the model long-term drift should be subtracted from all relevant quantities before the analysis. In our case, our region of interest does not go below the ocean thermocline, so we have not opted for this option.
Historical simulations use prescribed CO 2 atmospheric concentrations to simulate the recent past climate conditions, from the mid nineteenth century to present day. For all CMIP5 historical runs, the historical period ranges from 1850 to 2005 (Taylor et al. 2012). CMIP6 historical runs consider a longer period (1850Eyring et al. 2016). Historical experiments include both natural (such as solar  variability and volcanic eruptions) and anthropogenic (such as greenhouse gases and aerosol) forcing, trying to simulate the behaviour of the climate system over the last 150 years. Finally, 1pctCO2 is an idealised future scenario, in which the CO 2 atmospheric concentration is increased annually by 1% starting from pre-industrial values. This configuration allows to study the response of coupled models under idealised forcing. There is a small difference between CMIP5 and CMIP6 in the 1pctCO2 experimental setup: in the former, the annual 1% increase is limited to the first 140 years of simulation, by that time the CO 2 concentration is quadrupled with respect to pre-industrial values. Instead, in CMIP6 1pctCO2 experiments the CO 2 concentration keeps on increasing (beyond quadrupling) for the whole simulation (Eyring et al. 2016, see Appendix A1.4). For this reason, in the analysis we consider only the first 140 years of simulation.

Ocean reanalyses
Ocean reanalyses are obtained by running a global ocean model simulation constrained with observations (where available), through a data assimilation scheme. They are essential products for studying historical changes of the ocean, and they are valuable tools for exploring the role of the ocean on climate. A summary of recent improvements in ocean reanalyses, along with an evaluation of their uncertainties, is provided by Balmaseda et al. (2015) and Storto et al. (2019).
We compare results from historical simulations with four ocean reanalysis products: SODA 2.2.4 (Carton and Giese 2008), ORAS4 (Balmaseda et al. 2013), ORAP5 (Zuo et al. 2017), and ORAS5 ). Selected reanalysis products span different time periods, overlapping over the last decades of the twentieth century. Despite the rather poor reproduction of observed ocean properties by SODA 2.2.4 before the 1950s (Giese and Ray 2011; Chen and Wu 2012; Ray and Giese 2012), we decided to use the full dataset in order to assess STCs behaviour in the late nineteenth century and in the first half of the twentieth century.

Computation of mass and heat transport
Our analysis is based on mass and heat transport diagnostics, as well as SST and zonal wind stress. The total meridional mass transport (in Sverdrups; 1 Sv = 10 6 m 3 s −1 ) is evaluated as where 1 , 2 is the longitudinal extension over which the integration is performed, h is the ocean depth, is the sea surface, and v includes both resolved and parameterised meridional velocity components. We only considered the zonally and vertically integrated equatorward meridional transports in the Pacific basin, and in the uppermost 1000 m as in Graffino (2019) and Graffino et al. (2019). Other common metrics include computing the meridional transport in the pycnocline (see for example Capotondi et al. (2005) and Zhang and McPhaden (2006)), or tracking the maximum of the meridional streamfunction (as in Park et al. (2009) and He et al. (2019)). In our opinion, all these methods present shortcomings. The pycnocline method needs to be tuned for the relevant density range over which the meridional transport is occurring. This is hard to do for a large number of models, because the meridional transport occurs over different density ranges in different models (not shown). Instead, considering a large depth range, while masking out all non-equatorward meridional transport, allows to easily evaluate the STCs transport for all models. On the other hand, the streamfunction method does not allow for an evaluation of the Pacific STCs, because the integration must be performed over the whole Indo-Pacific basin. Due to the conformation of the Indian Ocean, the Indian STCs dynamics is very different from the Pacific one (Schott et al. 2004). This aspect is particularly relevant for the computation of the meridional transport in the Southern Hemisphere. Also, the streamfunction includes the aforementioned Tropical Cells, and caution must be used when tracking the maximum value of the streamfunction.
We employ the meridional energy transport diagnostic introduced in Graffino et al. (2019) and also used in Graffino (2019), for assessing the meridional heat transport associated with the STCs. The diagnostic relies on the method developed by Klinger and Marotzke (2000), allowing for the computation of the STCs-related meridional energy transport using Ekman dynamics. The main assumptions are that the STCs are wind-driven structures, and that the equatorward flow is isothermal. As Klinger and Marotzke (2000) point out, these conditions are usually not met in realistic numerical simulations. However, the simplicity of this method makes it helpful in qualitatively assessing the STCs meridional transport, as well as its changes over time.
The expression for the STCs meridional energy transport is where C p = 3992.1 J kg −1 • C 1 is the heat capacity for seawater at constant pressure, M E = − (y)∕f (y) is the Ekman mass transport obtained as the ratio between the surface wind stress and the Coriolis parameter, and is SST. The derivation of Eq. 2 can be found in Graffino et al. (2019, see Appendix). The energy transport is integrated meridionally between y 1 = 10 • and the latitude of zero wind stress (around 30 • for all models and datasets), although contributions to the STCs mass transport in the real ocean can also come from more poleward locations (McCreary and Lu 1994;Liu et al. 1994).

STCs mass transport and relationship with equatorial Pacific SST
In this section, we evaluate the low-frequency variability of the Pacific STCs mass transport, and its relationship with the equatorial Pacific SST. We also assess the time-mean circulation associated with the cells. This is done by computing the equatorward mass transport (Eq. 1), and separating western boundary and interior contributions. According to Zhang and McPhaden (2006), the interior component is computed from the South American coastline to 145 • E in the Northern Hemisphere, and to 165 • E in the Southern Hemisphere. All transport are evaluated at 9 • N and 9 • S, as in McPhaden and Zhang (2002Zhang ( , 2004 and Zhang and McPhaden (2006). Linear trends are removed before doing the analysis. Being interested in low-frequency variability, a 6-year low-pass filter is applied on all time series as in Zhang and McPhaden (2006), unless otherwise specified. The relationship between STCs interior transport and equatorial Pacific SST is well established in both observations and models (McPhaden and Zhang 2002Zhang , 2004Zhang and McPhaden 2006;Farneti et al. 2014a;Farneti 2017;Graffino et al. 2019). To explore such a link, we compute the linear correlation coefficient between interior mass transport convergence (evaluated as the sum of the southward mass transport at 9 • N and the northward mass transport at 9 • S) and the SST averaged over the central and eastern equatorial Pacific Ocean (9 • N-9 • S, 90 • -180 • W). In general, the largest correlation is found at lag 0 for most models and reanalysis products (not shown). Therefore, all correlation coefficients are computed at lag 0. Statistical significance of the coefficients is evaluated by using the z-transformation method (Fisher 1992) at the 0.05 level, with effective degrees of freedom computed according to Sun et al. (2019, see Supporting Information).

PiControl simulations
We start by looking at the piControl simulations. A common time period (the first 500 years) is chosen for all models. As shown in Fig. S1a, linear correlation coefficients are ranging between −0.73 (GFDL-CM2.5-FLOR and GFDL-ESM2M) and −0.1 (GFDL-CM3.0). Apart from GFDL-CM3.0, they are all statistically significant. The relationship between STCs transport convergence and equatorial SST is also shown in the time series in Fig. S2. Analysing the timemean transport convergence (Fig. S1b), the total mass convergence is in the 55-60 Sv range for most models, with largest (63.1 Sv) and weakest (41.7 Sv) simulated by GFDL-ESM2M and GFDL-CM2.5, respectively. Values concerning time-mean interior transport show a smaller spread, all lying between 10 and 15 Sv. However, in terms of percentage, total and interior transport show a similar range of variation across models (around 30%).
It is also evident that interior mass convergence from piControl simulations show large variations at decadal timescales. The standard deviation of STCs interior mass convergence and equatorial SST (Fig. S3a) allows to better evaluate this aspect. Notably, there are two "clusters" of models in terms of SST standard deviation: GFDL-CM2.1, GFDL-CM2.5-FLOR, and GFDL-ESM2M all lie around 0.8 • C , while GFDL-CM2.5, GFDL-CM3.0, and GFDL-CM4.0 have a standard deviation of ∼ 0.5 • C . In terms of interior convergence standard deviation, the models are spread between 2.8 and 4.3 Sv, with three models (GFDL-CM2.5, GFDL-CM2.5-FLOR, and GFDL-CM4.0) lying around 3.3 Sv. However, our subset of piControl runs is too small to infer a definite relationship between the standard deviation of the STCs interior transport convergence and STCs-SST anti-correlation (i.e. negative correlation; Fig. S3b).

Historical and 1pctCO2 simulations from coupled models, and ocean reanalysis products
Moving to historical and 1pctCO2 runs from CMIP5 models, Fig. 1 shows the correlation between STCs interior transport convergence and equatorial SST (panels a and c), and the time-mean STCs interior transport convergence (panels b and d). Figure 1a, b also include ocean reanalysis products. Focusing on STCs-SST correlations, there is a large spread of values among CMIP5 models in both historical and 1pctCO2 configurations (respectively in Fig. 1a, c), although they all give negative values of correlation. A smaller spread is shown by ocean reanalyses. It is also worth to note that the coefficients do not change significantly from historical to 1pctCO2 runs for most models, with notable exceptions being ACCESS1-0 (larger negative correlation) and CNRM-CM5 (smaller negative correlation). On average, the mean STCs-SST correlation is −0.72 for the reanalyses, −0.58 for the historical runs, and −0.59 for 1pctCO2 runs. Time-mean STCs transport convergence from CMIP5 models also shows a large spread (Fig. 1b, d) for both total and interior contribution. Here too, values are very similar between historical and 1pctCO2 runs. It is evident that ocean reanalyses give larger interior transport contribution than coupled models, with the only exception being MPI-ESM-MR, and larger anti-correlations between STCs interior transport convergence and equatorial Pacific SST.
CMIP6 models show an even larger intermodel spread, in terms of STCs-SST anti-correlations (Fig. 2a, c), with respect to CMIP5 models. By looking at the linear correlation coefficients, we can say that, on average, the selected CMIP6 models reproduce a weaker STCs-SST relationship than ocean reanalysis products and corresponding CMIP5 models. The mean STCs-SST correlation is −0.46 for historical runs, and −0.39 for 1pctCO2 runs.
Time-mean STCs mass transport convergence for CMIP6 models has a comparable intermodel spread (Fig. 2b, d) with CMIP5 models. Computed values are again very similar from historical and 1pctCO2 configurations, but there is no indication of a relationship between STCs-SST anti-correlation and time-mean transport convergence (neither total nor interior contribution). Again, ocean reanalysis products give the largest STCs-SST anti-correlations.
To summarise, correlation coefficients shown in Figs. 1, 2 and listed in Table 2 give a large spread in the STCs-SST relationship across historical CMIP5 and CMIP6 simulations. Ocean reanalysis products give better results with a mean correlation of −0.72 , close to the observed value of −0.85 (Zhang and McPhaden 2006). In this respect, the Fig. 1 a Linear correlation coefficients at lag 0 of interior mass transport convergence at 9 • of latitude (Sv) and equatorial Pacific SST anomaly ( • C ; 9 • N-9 • S , 90 • -180 • W ), as simulated by ocean reanalysis products and CMIP5 historical simulations. Linear trend is removed before the computation, and a 6-year low-pass filter is applied to both mass transport and sea surface temperature time series. For each model, 95% confidence interval is shown. Apart from ORAP5 and MPI-ESM-P, all coefficients are significant at the 0.05 significance level. b Time-mean total equatorward mass transport convergence at 9 • of latitude (Sv), separated in interior and western boundary contributions, as simulated by ocean reanalysis products and CMIP5 historical simulations. Black bars represent interior transport, grey bars represent western boundary transport. c Same as (a), but for CMIP5 1pctCO2 simulations. Apart from ORAP5 and MPI-ESM-P, all coefficients are significant at the 0.05 significance level. d Same as (b), but for CMIP5 1pctCO2 simulations choice of using the full-length SODA 2.2.4 dataset has a small but not negligible negative impact on its performance, as the pre-1950 STCs-SST anti-correlation is slightly weaker than the more recent one (not shown). Looking at the listed historical simulations, the multimodel-mean correlation is −0.58 for CMIP5 models and −0.46 for CMIP6 models.
Time series of STCs interior convergence transport and equatorial Pacific SST anomalies are shown in Figs. 3, 4 for historical CMIP5 and CMIP6 runs respectively, and are compared with time series obtained from ocean reanalysis products. In both figures, we can see how differently models reproduce the STCs-SST relationship, already highlighted by the correlation analysis in Figs. 1, 2.
Looking at the time series in Figs. 3, 4, large intermodel variations can be seen in terms of STCs interior transport convergence variability, for both CMIP5 and CMIP6 models. The same can be said for 1pctCO2 simulations, shown in Figs. S4, S5. To better evaluate this aspect, we show in Fig. 2 a Linear correlation coefficients at lag 0 of interior mass transport convergence at 9 • of latitude (Sv) and equatorial Pacific SST anomaly ( • C ; 9 • N-9 • S , 90 • -180 • W ), as simulated by ocean reanalysis products and CMIP6 historical simulations. Linear trend is removed before the computation, and a 6-year low-pass filter is applied to both mass transport and sea surface temperature time series. For each model, 95% confidence interval is shown. Apart from ORAP5, ACCESS-CM2, BCC-ESM1, CanESM5, and MPI-ESM1-2-LR, all coefficients are significant at the 0.05 significance level. b Time-mean total equatorward mass transport convergence at 9 • of lat-itude (Sv), separated in interior and western boundary contributions, as simulated by ocean reanalysis products and CMIP6 historical simulations. Black bars represent interior transport, grey bars represent western boundary transport. c Same as (c), but for CMIP6 1pctCO2 simulations. Linear trend is removed before the computation, and a 6-year low-pass filter is applied to both mass transport and sea surface temperature time series. Apart from BCC-CSM2-MR, BCC-ESM1, MPI-ESM1-2-LR, and MRI-ESM2-0, all coefficients are significant at the 0.05 significance level. d Same as (b), but for CMIP6 1pctCO2 simulations panels a, c of Figs. 5, 6 the relationship between standard deviation of STCs interior transport convergence and equatorial Pacific SST, for CMIP5 and CMIP6 models respectively. Values obtained from ocean reanalysis products are shown in historical runs plots for comparison. A large spread of values is shown by both CMIP5 and CMIP6 models, in terms of STCs interior transport convergence. Concerning the equatorial SST anomaly, there is more agreement among models for what concerns its standard deviation, especially for CMIP5 models. In general, models showing a large STCs mass transport variability also exhibit large equatorial Pacific SST variations at decadal timescales; this aspect was also highlighted by Zhang and McPhaden (2006). Focusing on the multimodel-mean standard deviation of the STCs interior transport convergence, CMIP5 models average at 2.8 Sv (2.72 Sv) for historical (1pctCO2) simulations, while the value for CMIP6 models is smaller for both historical (2.47 Sv) and 1pctCO2 simulations (2.39 Sv). However, the largest values are shown by ocean reanalysis products, averaging at 5.03 Sv.
Panels b, d of Figs. 5, 6 focus instead on the relationship between the standard deviation of STCs interior transport convergence and the correlation between STCs interior transport convergence and equatorial Pacific SST. We explored this possibility to find a reason why reanalyses better reproduce the observed link between STCs transport and equatorial SST with respect to models, and why some models reproduce this link better than others. From what we see in Figs. 5, 6, the variability in the STCs interior transport convergence might explain why this is the case. In fact, values obtained from ocean reanalyses (grey markers in Figs. 5b, 6b) stand out on average from those related to coupled models. Furthermore, this is also valid within the same dataset in the case of CMIP6 historical runs. Farneti (2017) already noted this aspect, while analysing CMIP5 historical simulations and comparing them with observations, as well as with ocean-only runs forced with atmospheric reanalysis. Farneti (2017) argued that coupled models underestimate the subtropical wind stress variability. Because of that, Pacific STCs mass transport variability is underestimated as well, thus providing a too weak forcing to the equatorial Pacific SST.
We now investigate the SST pattern related with the Pacific STCs. Many studies highlighted the role of STCs variability on setting Pacific SST decadal variability in both observations and models (e.g.; Solomon et al. 2003;Capotondi et al. 2005;Zhang and McPhaden 2006;Cheng et al. 2007;Farneti et al. 2014b) In Fig. 7 we show whether this is consistent between ocean reanalyses and selected coupled models. Plots obtained from models are shown as multimodel mean for convenience; in order to do so, results are remapped on a 1 • × 1 • horizontal grid before the computation.
All zero-lag correlation patterns shown in Fig. 7 resemble the well-known PDO pattern (Mantua et al. 1997), one of the main driver of low-frequency variability in the Pacific Ocean. Several works already found a connection between STCs variability and PDO (Farneti et al. 2014a;Hong et al. 2014), as well as with the most recent PDO regime shifts Zhang 2002, 2004). Furthermore, SST adjustment in the eastern and central sector of the Pacific Ocean are consistent with theories about the STCs convergence (Zhang and McPhaden 2006;Farneti et al. 2014a;Graffino et al. 2019). However, the extratropical SST pattern may arise from additional SST adjustments, due to atmospheric tropical-extratropical teleconnections and extratropical air-sea interactions (e.g., Newman et al. 2003).
Focusing on individual models, the consistency of the SST pattern among models is striking, as shown in Figs. S6, S7 for historical runs from CMIP5 and CMIP6 models respectively. In general, models characterised by low STCs-SST anti-correlation also show the least "canonical" PDO pattern in their spatial correlation response, as can be seen by comparing Figs. S6, S7 with the linear correlation coefficients listed in Table 2. Correlation patterns from 1pctCO2 runs show no remarkable difference from their historical counterpart (not shown).

STCs mass and energy transport in future scenario simulations
In Sect. 3.1 we showed that the low-frequency variability of the Pacific STCs in 1pctCO2 runs is similar to that obtained under historical conditions. However, we did not mention yet about changes occurring to the STCs transport. By using Eq. 2, we compute the STCs energy transport for 1pctCO2 simulations, and we compare that with historical runs. In this case, time series are not detrended. We are showing in Fig. 8 the diagnostics for CMIP5 and CMIP6 simulations, focusing on the North and South Pacific STC. By looking at the cumulative meridional energy transport (that is, the value reached by the curve at the central vertical line of each plot), we can see how the two cells behave differently in the multimodel mean. According to Fig. 8a, the northern (southern) cell experiences a weakening (strengthening) of its meridional energy transport in 1pctCO2 runs, with larger and more consistent changes in the Southern Hemisphere. This response is also found in CMIP6 models (Fig. 8b). This different hemispheric behaviour of the Pacific STCs in a warming climate was shown by previous coupled modelling efforts (Lohmann and Latif 2005;Park et al. 2009). Our analysis also confirms what found by Graffino (2019), in which STC meridional energy transport from piControl and 1pctCO2 experiments was compared in a subset of CMIP5 models. It is worth to stress that in the present study the comparison is made with respect to historical simulations, including part of the forced response, and not to (unforced) piControl runs as in Graffino (2019).
There are two major players in Eq. 2: the meridional Ekman transport, and the meridional SST gradient. In Figs. S8, S9 the two terms are shown for both CMIP5 and CMIP6 multimodel means. We are showing the absolute value of both quantities, because we are interested in comparing historical and 1pctCO2 values. Starting from the meridional Ekman transport (Fig. S8), we can see that most of the models give a weakening (strengthening) of the North (South) Pacific transport. This is driven by analogous changes of the zonal wind stress, since the meridional Ekman transport is equal to − (y)∕f (y) . However, the pattern of zonal wind stress change is not zonally uniform for some models (not shown), which might complicate the interpretation of the results. Moving to the meridional SST gradient (Fig. S9), we see the same tendency: smaller (larger) gradient over the North (South) Pacific Ocean. Therefore, both terms of Eq. 2 contribute to the response seen in Fig. 8.
We compare these findings with trends of STCs mass transport, employing Eq. 1 as we did in Sect. 3.1. Unlike the STC meridional energy transport, in this case we have no univocal response about future STCs changes. Both CMIP5 and CMIP6 models (Figs. 9a and 10a, respectively) show a Table 2 Linear correlation coefficients at lag 0 relating interior mass convergence and equatorial Pacific SST as simulated by pre-industrial control, historical, and 1pctCO2 coupled models simulations, as well as ocean reanalysis products Linear trend is removed, and a 6-year low-pass filter is applied to both mass transport and sea surface temperature time series. Mean correlation coefficients for all datasets are also shown. Apart from those marked with an asterisk, all coefficients are significant at the 0.05 significance level piControl Historical 1pctCO2 There is no agreement in this case among models, although the majority of them gives a weakening trend. This is at odds with the meridional energy transport changes shown in Fig. 8, but in agreement with what found by Wang and Cane (2011). Looking at the INT component, panels c, d of Figs. 9, 10 show that there is no agreement among models for changes Fig. 3 Pacific STCs interior mass transport convergence anomaly at 9 • of latitude (Sv, thick line), and equatorial Pacific SST anomaly ( • C , thin line), as simulated by ocean reanalysis products and CMIP5 historical simulations. All time series are detrended and a 6-year low-pass filter is applied to both mass transport and sea surface temperature time series. The transport anomaly scale is inverted to highlight its anti-correlation with the SST anomaly Fig. 4 As in Fig. 3, but for CMIP6 historical simulations at 9 • N , while most models give a weakening trend at 9 • S . This negative trend is partially compensated in many cases by the WB component (Figs. 9f, 10f). There is also a robust indication for a weakening trend of the WB component at at 9 • N (Figs. 9e, 10e).
In terms of equatorward mass transport convergence (i.e. the sum of the northern and the southern cells contributions), we find that a general (but not univocal) weakening of the INT component is obtained for both CMIP5 and CMIP6 models (Figs. S10b, S11b respectively). However, there is no agreement across models about a compensation provided by the WB component (Figs. S10c, S11c), leading to an overall weakening of the total mass transport convergence in 1pctCO2 conditions (Figs. S10a, S11a). This overall weakening of the Pacific STCs was obtained by other studies in terms of mass (Merryfield and Boer 2005;Lohmann and Latif 2005;Wang and Cane 2011) and heat transport (Levine and Schneider 2011;Zelinka and Hartmann 2012;He et al. 2019), while others found no reduction in the overall STCs strength in a warming climate (Luo et al. 2009;Park et al. 2009). Furthermore, Wang and Cane (2011) also found a lack of western boundary compensation to the weakening interior transport component, and ascribed the STCs mass transport reduction to changing wind stress forcing in the Fig. 5 a Standard deviation of equatorial Pacific SST anomaly ( • C ) as a function of standard deviation of interior mass transport convergence at 9 • of latitude (Sv), as simulated by ocean reanalysis products (grey markers) and by CMIP5 historical simulations (black markers). Linear trend is removed before the computation. b Linear correlation coefficients at lag 0 computed between STCs interior mass transport convergence at 9 • of latitude and equatorial Pacific SST, as a function of standard deviation of STCs interior mass transport convergence at 9 • of latitude (Sv), as simulated by ocean reanalysis products (grey markers) and by CMIP5 historical simulations (black markers). The vertical axis is inverted. c As in (a), but for CMIP5 1pctCO2 simulations. d As in (b), but for CMIP5 1pctCO2 simulations warming scenario. For what concerns the equatorial Pacific SST, all models give a positive trend (Figs. S10d, S11d), although with different magnitudes. However, disagreeing trends in both equatorward mass transport components and equatorial Pacific SST are obtained for historical conditions as well (not shown). As noted by Zhang and McPhaden (2006), this aspect simply highlights that other processes, other than the STCs variability, contribute to the modelled warming trend in the equatorial Pacific Ocean.

Discussion and conclusions
An analysis of the low-frequency variability of the Pacific STCs is presented in this study. We employed state-of-theart coupled models included in the phase 5 and 6 of the Coupled Model Intercomparison Project (CMIP) listed in Table 1. After checking that pre-industrial control simulations can successfully generate STCs variability at decadal timescales under unforced conditions, we proceeded on diagnosing mass and heat meridional transports, and their statistical relationship with the equatorial SST, for historical simulations. As reference, we used several ocean reanalysis Fig. 6 a Standard deviation of equatorial Pacific SST anomaly ( • C ) as a function of standard deviation of interior mass transport convergence at 9 • of latitude (Sv), as simulated by ocean reanalysis products (grey markers) and by CMIP6 historical simulations (black markers). Linear trend is removed before the computation. b Linear correlation coefficients at lag 0 computed between STCs interior mass transport convergence at 9 • of latitude and equatorial Pacific SST, as a function of standard deviation of STCs interior mass transport convergence at 9 • of latitude (Sv), as simulated by ocean reanalysis products (grey markers) and by CMIP6 historical simulations (black markers). The vertical axis is inverted. c As in (a), but for CMIP6 1pctCO2 simulations. d As in (b), but for CMIP6 1pctCO2 simulations products. We also compared results from historical simulations with outcomes from 1pctCO2 condition, an idealised global warming scenario.
Generally speaking, caution must be used when studying low-frequency variability across model runs with different lengths, as in the case of our piControl simulations. In fact, for some models different time windows of the same simulation give different correlation coefficients (not shown), implying that assessments of the STCs-SST relationship critically depends on the considered time window and on the length of the simulation. At the same time, available historical records may not be a representative sample for studying long-term climate variability, as suggested by Wittenberg (2009) for ENSO statistics.
That said, the linear correlation between the interior mass transport convergence and the equatorial SST is important Fig. 7 Spatial correlation at lag 0 between the Pacific STCs interior mass transport convergence at 9 • of latitude (Sv) and the Pacific SST ( • C , thin line), as simulated by ocean reanalysis products (top two rows), and historical (third row) and 1pctCO2 simulations (bottom row) from CMIP5 and CMIP6 models. All fields are detrended and a 6-year low-pass filter is applied. Results obtained from coupled models are remapped on a 1 • × 1 • grid before computing the multimodel mean to characterise the influence of the STCs on the equatorial ocean state at decadal timescales. It is not trivial to disentangle the connection between forcing mechanism (the subtropical wind stress) and final response (e.g., mass transport convergence at the equator), but there are few interesting hypothesis in literature. Solomon and Zhang (2006) and Zhang and McPhaden (2006) pointed at the subtropical wind stress variability, deemed to be too weak in models to drive a strong STCs response. In particular, they refer to a too-weak Ekman pumping in coupled models. We could not find such feature in our simulations, although ocean reanalyses tend to give larger meridional Ekman transport than coupled models (not shown).
At the equator, the STCs-SST connection is explained by the v ′ T mechanism (Kleeman et al. 1999): a strengthened (weakened) STCs circulation drives a cold (warm) response at the equator; thus, we expect to see an anti-correlation (i.e. negative correlation) between STCs interior transport convergence and equatorial Pacific SST. We computed correlation coefficients on detrended, low-pass filtered time series for all considered models and products; they are shown in Figs. 1, 2 and listed in Table 2. Averaging these values, we The shading shows the intermodel spread. The left-hand side shows the South Pacific, and the right hand side shows the North Pacific. b Same as (a), but for CMIP6 models obtained that ocean reanalyses give larger anti-correlations ( −0.72 ) than both historical CMIP5 ( −0.58 ) and CMIP6 models ( −0.46 ). However, some coefficients computed from CMIP6 models are not statistically significant; discarding models with non significant coefficients, CMIP6 multimodel-mean correlation becomes −0.6 , slightly larger (i.e. more negative) than the CMIP5 value. This confirms that coupled models do not agree on the strength of the STCs-SST relationship, although, in all cases but one, a stronger (weaker) STCs circulation leads to cold (warm) equatorial SST. Furthermore, the disagreement is even larger within CMIP6 models than within CMIP5 models.
Comparing these results with previous studies, our analysis shows little improvement in coupled models with respect to what obtained by Zhang and McPhaden (2006). All the considered CMIP5 models and 8 out of 12 CMIP6 models show statistically significant anti-correlation, against 15 out of 18 runs analysed by Zhang and McPhaden (2006). Further, the multimodel-mean correlation is −0.58 and −0.46 ( −0.6 , if we remove non significant coefficients) for CMIP5 and CMIP6 models respectively, whereas it is −0.58 in Zhang and McPhaden (2006). Coupled models also give a not coherent response in terms of the well-known negative correlation between INT and WB components, an aspect already highlighted by Zhang and McPhaden (2006) and Wang and Cane (2011). This is resumed in Table 3, showing the linear correlation coefficients between INT and WB components at 9 • in the two hemispheres. It is evident that the WB-INT compensation is better reproduced by ocean reanalyses than coupled models. Also, a stronger agreement is found among CMIP5 than among CMIP6 models. However, several coefficients are not statistically significant. In fact, the compensation mechanism works much better at 9 • S than at 9 • N ; this effect was already noted by Lee and Fukumori (2003), and explained in terms of the different magnitude of the off-equatorial wind stress curl in the two hemispheres.
If we consider the time-mean STCs mass transport convergence (that is, the sum of the equatorward mass transport at 9 • N and 9 • S ), we see a large intermodel spread of values across models for both CMIP5 and CMIP6 datasets (see Figs. 1, 2). Comparing CMIP5 values with ocean reanalysis products, we see that some models give similar total transport convergence (like the GFDL models), but few of them match the reanalyses in terms of interior transport convergence. This might explain why reanalyses give larger STCs-SST anti-correlations than models, and why some models perform better than others in this respect (look for example at MPI-ESM-MR in Fig. 1a, b). However, this indication lacks robustness: for example, CanESM2 gives a much smaller interior transport convergence and, at the same time, a large STCs-SST anti-correlation. Also, the time average of the total STCs mass transport convergence for most models is larger than any previous estimate. The different method employed to compute the equatorward mass transport might be the cause. However, values of the interior component are much closer to previous estimates from observations Zhang 2002, 2004;Schott et al. 2004), assimilation products and reanalyses (Schott et al. 2007(Schott et al. , 2008Hong et al. 2014), and coupled models (Zhang and McPhaden 2006;Luo et al. 2009;Farneti 2017).
Moving to the STCs interior transport convergence variability at decadal timescales, we analysed the relationship between the standard deviation of the STCs interior transport convergence with the standard deviation of equatorial Pacific SST, as well as with the STCs-SST anti-correlation (Figs. 5,6). Large variations are seen across models; this was already evident by looking at the time series in Figs. 3, 4. Values of standard deviation of both STCs interior transport convergence and equatorial SST are larger than estimated by Zhang and McPhaden (2006) and Farneti (2017). Furthermore, it is suggested that, averaging over datasets, a larger standard deviation of STCs interior transport convergence gives a larger STCs-SST anti-correlation. This is also true within the same dataset in the case of CMIP6 models. Another interesting hint is given by the relationship between the variability of STCs interior transport convergence and equatorial Pacific SST (already described by Zhang and McPhaden (2006)). By computing regression lines of the standard deviations of these two quantities, we see that the slope does not change between historical and 1pctCO2 runs for CMIP5 models, while it gets more positive from historical to 1pctCO2 conditions for CMIP6 models. In that case, in a global warming scenario, we can expect a larger equatorial Pacific SST decadal variability from models showing larger STCs interior transport convergence variability.
Another interesting aspect comes from the relationship between Pacific STCs and PDO, a coupled ocean-atmosphere mode generating decadal variability in the Pacific Ocean. This was already noted by Zhang and McPhaden (2006), with Farneti et al. (2014b) and Hong et al. (2014) also providing different feedback mechanisms to explain this relationship. The correlation pattern shown in Fig. 7 is remarkably robust across reanalyses and historical coupled runs. By looking at individual models (Figs. S6, S7) Fig. 9 a Decadal trend (Sv/decade) of the total (interior plus western boundary contribution) STCs mass transport at 9 • N , for CMIP5 1pctCO2 runs. Please note that a positive trend corresponds to a weakening of the equatorward mass transport. b Decadal trend (Sv/ decade) of the total (interior plus western boundary contribution) STCs mass transport at 9 • S , for CMIP5 1pctCO2 runs. In this case, a positive trend corresponds to a strengthening of the equatorward mass transport. c Same as (a), but for the interior contribution. d Same as (b), but for the interior contribution. e Same as (a), but for the western boundary contribution. f Same as (b), but for the western boundary contribution ◂ Fig. 10 As in Fig. 9, but for CMIP6 1pctCO2 simulations and comparing them with coefficients listed in Table 2, it is evident that models better reproducing the STCs-SST relationship also show a more "canonical" PDO pattern in historical conditions. This does not change when considering non-detrended fields, while it does in the 1pctCO2 configuration (not shown).
We employed two different approaches to study Pacific STCs future changes in 1pctCO2 conditions. The first approach is based on a method developed by Klinger and Marotzke (2000), and successfully employed by Graffino et al. (2019). Based on Ekman dynamics, it computes the meridional heat transport by STCs using only surface properties (zonal wind stress and SST). Results obtained with this method are visually shown in Fig. 8, and suggest that Pacific STCs are going to change differently in the two hemispheres, with a weakening (strengthening) of the northern (southern) cell, as also found by previous studies (Lohmann and Latif 2005;Park et al. 2009;Wang and Cane 2011). Unfortunately, this diagnostic does not allow for an easy comparison with previous works, nor for an analysis of what Table 3 Linear correlation coefficients at lag 0 relating interior and western boundary mass convergence contributions as simulated by preindustrial control, historical, and 1pctCO2 coupled models simulations, as well as ocean reanalysis products PiCntrl refers to the PiControl runs, Hist to the historical runs, and 1pCO2 to the 1pctCO2 runs. Linear trend is removed, and a 6-year low-pass filter is applied to both mass transport and sea surface temperature time series. Mean correlation coefficients for all datasets are also shown. Apart from those marked with an asterisk, all coefficients are significant at the 0.05 significance level it is causing this hemispheric difference in behaviour. However, we could verify that both meridional Ekman transport and meridional SST gradient collaborate on shaping the STCs future changes in both hemispheres (Figs. S8, S9). A prominent role for the subtropical wind stress in STCs future changes was already predicted by several authors (Luo et al. 2009;Park et al. 2009;Wang and Cane 2011). The second approach involves the computation of equatorward mass transport trends. In this case, by looking at individual contributions coming from each cell (Figs. 9, 10), we confirm a general reduction of the equatorward total transport for the northern cell. There is less agreement about the southern cell, which on average also shows a weakening trend, as opposed to our meridional energy transport metrics. However, by looking at trends, we see that most models agree on a weakening (strengthening) of the northern (southern) contribution from the western boundary component, and on a general weakening of the interior component. About this last point, our results are consistent with the Pacific Ocean pycnocline analysis of Luo et al. (2009) and Wang and Cane (2011). On the other hand, a final answer about the role of the two components (WB and INT) in STCs future changes is hard to find in literature, due to the large range of response from coupled models. In any case, the disagreement between our two approaches requires further attention. Our energy transport diagnostics need to be compared with the total meridional heat transport computed by the models, as well as its sensitivity to different components of the STCs dynamics.
Finally, we confirmed that coupled models do not agree on several aspects of the simulated low-frequency variability of the Pacific STCs, as well as their relationship with equatorial Pacific SST. However, we believe that we made some progress in explaining why this is the case. In fact, the interesting relationship between STCs interior transport convergence and STCs-SST anti-correlation suggests that, underestimating the mass transport variability, coupled models might lose part of the SST-driving forcing.
The STCs are wind-driven overturning structures. Thus, a more accurate study of the subtropical wind stress as reproduced by coupled models is needed. Also, fast variations of the zonal wind stress, and consequent rapid changes of the meridional Ekman transport, might not be captured by using annual means. Future efforts will have to employ data at higher temporal resolution to better clarify this aspect. Furthermore, STCs interior transport convergence and equatorial Pacific SST show significant negative correlations even without low-pass filtering, and the resulting correlation pattern resembles the canonical ENSO shape (not shown). The influence of Pacific STCs on modulating ENSO variability was approached by several works in the past, but it would be useful to revisit some of these concepts with new available models and reanalysis products.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.