In this study, we used the ocean general circulation model NOAA-GFDL Modular Ocean Model version 5 (MOM5) (Griffies, 2012; Todd et al., 2020). MOM5 is a free-surface primitive equation model and uses the z☆ -vertical coordinate. The horizontal resolution is 1° with refined meridional resolution equatorward of 30°. It uses 50 vertical levels, with 10 levels in the top hundred meters. Subgrid mesoscale processes are parameterized with the Gent-McWilliams skew-flux closure scheme (Gent and McWilliams, 1990; Gent et al., 1995; Griffies, 1998) and submesoscale eddy fluxes according to Fox-Kemper et al. (2008) and Fox-Kemper et al. (2011). Vertical mixing is represented with a K-profile parameterization (Large et al., 1994).

For our experiments, we followed the protocol proposed by the ocean-only Flux-Anomaly-Forced Model Intercomparison Project (FAFMIP; Todd et al., 2020) in order to exclude any coupled feedback. A spin-up was performed by integrating the model towards a near-equilibrium state, using as forcing fields the total mass flux of liquid precipitation plus liquid runoff minus evaporation, net heat flux, and surface wind stress (see Todd et al., 2020). An ensemble of pi-Control CMIP5 AOGCMs was used to generate climatologies for all fluxes. Restoring of both sea surface temperature (SST) and sea surface salinity (SSS) employs a time scale of 2 and 3 months, respectively (Todd et al., 2020).

After equilibrium, a restored Control simulation was obtained through an additional 70 years of integration. During this Control, total effective air–sea heat and freshwater fluxes were stored at daily intervals, including restoring-induced SSS fluxes. Finally, a twin second Control was obtained with SSS restoration deactivated and by prescribing the effective heat and freshwater fluxes diagnosed from the first Control. Differences between both controls have been quantified and are negligible. The flux-forced Control is used as the reference state, and will henceforth referred to as the Control simulation.

Sensitivity experiments are forced by FAFMIP surface flux perturbations, derived from a set of 13 CMIP5 AOGCMs, where each surface flux represents the transient response of the climate system to a doubled CO₂ concentration under the idealized scenario 1pctC02 (Gregory et al., 2016). Anomalous fields are computed as the difference between the climatological monthly time-means of years 61–80 of 1pctCO₂ and the corresponding time-means of the non-evolving and unforced variability equilibrium state known as pre-Industrial or piControl run (Taylor et al., 2011). All experiments are global, 70 years long, branching off from the Control simulation and forced by the Control effective surface fluxes plus the FAFMIP momentum (faf-stress) or heat (faf-heat) anomalous fluxes ( Figures 2 , 3 ). In faf-stress, global zonal and meridional wind stress perturbations are added to the momentum balance equation, whereas faf-heat adds heat flux perturbations to the total surface downward heat flux.

The Tropical Pacific Ocean is vertically divided into depth-dependent integrals representative of shallow and deep circulations. The shallow ocean, involving the subtropical overturning cells, spans in both hemispheres from the tropics up to the mid-latitudes where subduction of surface water takes place (McPhaden and Zhang, 2004), approximately 32° in our model. In the vertical, it is bounded at the top by the mixed layer and at the bottom by the deep ocean. In the latitudinal extent, the shallow ocean extend from near the Equator up to the latitude of minimum zonal wind stress (32°N and 30°S, respectively). Here, mixed layer depth estimation consists of the determination of the vertical level at which temperature change from the surface is equal to 0.5°C, as suggested by Levitus (1982). Our definition for the shallow overturn considers all vertical levels away from the influence of both the surface mixed layer and the deep ocean circulation.

Our definition of the Shallow Wind-driven Overturning Circulation differs from that of subtropical cells (STCs) by previous studies (McCreary and Lu, 1994; McPhaden and Zhang, 2002; Zhang and McPhaden, 2006; Farneti et al., 2014; Graffino et al., 2018). Traditionally, the STC is characterized in terms of equatorward mass transport within the upper thermocline. In contrast, our definition encompasses a wider latitudinal extent and includes not only the subtropical overturning cells but also the lower-latitude tropical overturn and equatorial upwelling region, and intentionally avoids the mixed layer region.

In the case of the deep ocean, this region is defined as the region below the deepest vertical limit of the shallow overturn ( Figure 1A ). As suggested by Talley (2003), we defined the base of the shallow overturn as the maximum subducting level at the subtropical gyre. The definition implies using meridional zonally integrated Sverdrup circulation sections and determining either the vertical level at which the transport is nearly zero or simply by determination of the deepest closed streamline. Once the mixed layer depth and the base of the overturn are determined, the separate contribution from the shallow circulation can be estimated, minimizing the influence of both atmosphere and the deep ocean.

Explicitly, the deep ocean encompasses all the vertical levels of the ocean below the shallow overturn up to near the bottom of the ocean. Given the minimization of the mixed layer influence and the determination of the shallow ocean, all the remaining circulation contributions must represent the circulation related to the deep branches of the MOC. However, the above-mentioned partitioning of the water column is only for the purpose of the determination of the budget integrals. In Section 2.3, we will provide a link between the heat budget analysis and an independent analysis of the dominant circulation patterns in the water column.

The following provides a diagnostic accounting of the ocean heat budget for a grid cell in the model. The relative contribution of each physical process to the overall heat budget, concerning the equilibrium state and perturbed climate conditions, is quantified using the three-dimensional process-based tendency diagnostics for prognostic temperature (Griffies et al., 2016; Dias et al., 2020). Such a procedure allows us to determine in physical terms the processes responsible of setting the overall heat storage balance.

The time tendency for the heat content in a grid cell per unit horizontal area, expressed as a semi-discrete finite volume formulation within the grid cell, is shown by Griffies et al. (2016)

where T represents the temperature of a grid cell, C p is the heat capacity coefficient, and ρ 0  d z is the mass of seawater per horizontal area in a grid cell with ρ 0 being the reference density and d z being the grid cell thickness. The tracer fluxes at the surface of the ocean are represented as Q . Generalized three-dimensional velocity is expressed as u while F accounts for the parameterized subgrid-scale processes, including circulation, diffusion, eddies, and vertical mixing.

The tracer formulation’s vertical extent is time-dependent, which implies that each heat budget term scales according to the thickness of the cell. In contrast, the horizontal cross-sectional area is time-independent, which implies area normalization if the tracer grid cell is multiplied by the horizontal area.

The time tendency for the heat content for each grid cell over an ocean column arises as the convergence of advective processes ( u   T ) in addition to subgrid-scale heat fluxes ( F ) (Eq. 1) and the net time tendency for the heat content is represented as the convergence of the general tracer flux F xz and is expressed as

Here, ∂ t T net represents the model diagnosed tracer for the time tendency in a grid cell, and is equivalent to the temperature time derivative in the advective–diffusive equation, or can be expressed in terms of a generalized tracer flux F , which is the convergence of advective or subgrid-scale fluxes crossing grid cell faces. The diagnosed advective processes account for the model prognostic advection ( ∂ t T adv ) and parametrized eddy-induced advection ( ∂ t T eddy ), including both mesoscale and submesoscale flux components (Fox-Kemper et al., 2008; Fox-Kemper et al., 2011; Gent and McWilliams, 1990). Within the diagnosed diffusive mixing processes, we distinguish among isopycnal ( ∂ t T iso ) and diapycnal mixing ( ∂ t T dia ). The first refers to diffusion oriented along neutral or isopycnal surfaces and the second encompasses vertical convection and vertical diffusion, including parameterization of diapycnal diffusion due to tidal mixing, shear instabilities, double diffusion, mixing due to convective instabilities, and the overflow parameterization of dense water flowing down topographic slopes.

Since our formulation does not take into account mixed layer depths, contributions related to the non-advective radiative transfer are neglected (i.e., Q = 0 ). In contrast, we included a diagnostic related to parameterized oceanic absorption of downward solar radiation ( ∂ t T swh ). This term stands for shortwave heat flux at the bottom of a tracer grid cell face, generally expressed as the attenuated downwelling shortwave radiative heat per unit area incident at the ocean surface. Owing to the similar nature of both contributions, diapycnal mix and shortwave heat penetration terms are grouped together and presented as a vertical diffusive processes term (vDIFF). The expression ∂ t T net stands for both the iterated model diagnosed term and the linear sum of all the different components. By construction, temperature tendency diagnostics should close the heat budget. However, both estimations are not necessarily equivalent implying a heat budget residual.

Imbalances can be explained in terms of two non-quantified contributions, namely, non-diagnosed processes (diffusive mixing) and non-resolved high-frequency variability (i.e., interannual variability). In particular, non-local turbulent vertical flux boundary layer parameterization, overflow schemes, and mixing across unresolved straits are contributors to such imbalance (Griffies et al., 2016). Despite the fact that such components are indeed diagnosed by the model, those quantities are not considered in our analysis, since they have been quantified to be small. In contrast, the component related to the so-called interannual variability is explained in terms of an inaccurate closure budget due to the time scale on which tendencies are saved. All diagnostics used here are saved as annual averages despite the fact that they are diagnosed in the model at each time step. In order to close the budgets, such contribution is added to the vertical diffusive component given the diffusive mixing nature of the non-diagnosed terms.

It is important to note that imbalances do not arise because of model drift. In the equilibrium state, the global mean net tendency is near-zero but not exactly zero, a common feature among global climate models due to the typically longer time scales required for the deep ocean to reach equilibrium (Exarchou et al., 2014). Such contributions are rather small with no significant impact on the determination of the budgets. In contrast, we refer to imbalances as any transient change (warming or cooling) in the ocean heat content time rate of change (heat storage) derived from the perturbed climate conditions.

Here, we are interested in grouping all of the advective terms together, so that the superresidual advection encompasses the contribution from the parameterized flow, and the super-residual advection (srADV) is defined as

Lastly, previous studies have adopted weighting functions in the time domain to assess changes between the first and last decade of the experimental period (Gregory et al., 2016; Todd et al., 2020). In our study, we are interested in the transient response of the ocean and none of the results make use of any type of weighting and the transient response of each heat budget is evaluated relative to the equilibrium state solution.

The Pacific Ocean shallow wind-driven overturning circulation is described in terms of the meridional mass transport. Our estimations consider mass balance adjusted to zero across 24°S and 24°N. As in Talley (2003), the selection of such latitudes is based on the detection of three circulation features: a strong overturning circulation (i.e., surface poleward and subsurface equatorward transports), a well-defined poleward WBC flow, and a mid-latitude subducting area. In addition, the selected latitude is meant to minimize the effects of water recirculation by the tropical subsurface overturn (between 5° and 10° latitudes).

Our results refer to the annually averaged meridional transport under both sensitivity cases and the equilibrium-state solution, with time-averaged spanning over the last decade of the simulation period. The zonally integrated transport for different density classes is also included. Each water mass is classified according to its ventilation level on the basis of climatological annual means (see Table 1 ). Procedure implies partitioning the water column in terms of the isopycnal distribution and their level of outcropping. The density limits between the different water masses are determined as the outcropping density level at the given latitude when possible, or simply based on the density field distribution.

As suggested by Talley (2003), four types of water classes are defined: upper, intermediate, deep, and bottom waters. The upper waters refer to low latitudes waters bounded by the outcropping isopycnal at 30°S–30°N. Despite the motion of this layer being arguably primarily wind-driven, we did not explicitly isolate such contribution, so transport inside this layer includes both Ekman poleward transport and the geostrophic transport related to the gyre circulation. The intermediate waters refer to mid-latitude waters bounded by outcropping isopycnal surfaces between the region 30° and 60°. The deep water is immediately located below the intermediate water and is only ventilated through the Southern Ocean. Bottom waters are located in the densest layers.

The overall volume transport is adjusted to zero mass balance, which means compensation of volume between layers. We did not estimate any water mass conversion rates, and our calculations are simply based on the assumption of effective subduction rates. The northward Ekman transport is zero at approximately 30°N and all northward Ekman transport at 24°N is at a much lower density than the maximum subducted density; hence, we assumed that all Ekman flow at 24°N loses buoyancy and joins the southward geostrophic flow of the subtropical gyre.

Additionally, within the shallow overturning cell, volume transport is divided into the Western Boundary Current (WBC) transport and interior equatorward transport ( Figures 1 , 4 ). The outcrop between equatorward interior flow and WBCs is meant to avoid regions of recirculation water masses (McPhaden and Zhang, 2004; Zhang and McPhaden, 2006). All of the Ekman transport is assumed to remain in the subducting subtropical gyre. To balance the returning flow, we arbitrarily use only the lightest and warmest part of the WBC upper layer. Remaining volume transport is assumed to flow outside of the gyre towards the subpolar region and to return as intermediate or deep water.

Heat transport (HT) is calculated from the zonally integrated transports in isopycnal layers. Expression for the lateral heat transport through a vertical cross-sectional area requires velocity and temperature and is expressed as

where ρ , v , C p , and θ are the density, velocity normal to the area, the specific heat of seawater, and potential temperature, respectively.

The volume-integrated heat budgets for the Northern Hemisphere (NH) and the Southern Hemisphere (SH) Tropical Pacific Ocean are shown in Figure 5 , with vertical integrations accounting for the total column, the Shallow Wind-driven Overturn, and the deep ocean (mixed layer depths not included). In the Control simulation ( Figures 5A, B ) the tropical Pacific Ocean heat tendency imbalance is near-zero but not exactly zero (NET), a common feature in global climate models due to the very long time scales required for the deep ocean to reach equilibrium (Exarchou et al., 2014).

In general, heat balance arises as the compensation between vertical diffusive warming (vDIFF) and the residual advection due to resolved (ADV) and parameterized eddy processes (EDDY), along with eddy mixing terms (ISO). As in Kuhlbrodt et al. (2015) and Dias et al. (2020), here we will refer to super-residual advection (srADV) as the sum of the advective terms and eddy mixing terms. In fact, the net balance can be simplified as a compensation between super-residual advection (srADV) cooling and vertical diffusive warming (vDIFF). The use of the super-residual advection framework allows us to simplify the description of the heat balance analysis. All the individual contributions to the heat balance (ADV, EDDY, rADV, and ISO) are presented in Figure 5 , but heat budget changes will be discussed only in terms of the super-residual advection framework (srADV and vDIFF).

The observed balance seems to disagree with the two opposite vertical regimes proposed by Dias et al. (2020), but it arises from excluding the influence from the mixed layer. Here, both regimes are included in the ocean interior defined by Dias et al. (2020). Our results are thus in accordance with the ocean interior regime proposed by Dias et al. (2020) but also suggest that most of the balance is achieved within the shallow wind-driven overturn. This is confirmed by the zonally integrated heat budget within the ocean column ( Figures 6A–C ), which shows a balance between downward warming by a vDIFF flux and upward cooling by a srADV flux. The balance holds over the whole column, but it is particularly intense within the upper levels of the shallow ocean, near the mixed layer, and within the subducting areas, suggesting a strong association with the convergence and subduction of surface tropical waters. This is in accordance with the circulation pattern held by the cell, which promotes the removal of upper ocean water towards the interior ocean (Dias et al., 2020).

According to our definition and partitioning of the water column, although we have excluded the surface processes, the warming of the rest of the column comes from the mixed layer. Transfer of heat comes from the mixed layer heat transfer from the mixed layer toward the interior ocean given the circulation pattern held by the cell. The last is by removing mixed layer water toward deeper levels, mainly along the isopycnal field toward the interior ocean (Dias et al., 2020). The last is actually in accordance with vDIFF warming and srADV cooling taking place mainly in the subducting area (30°N and 30°S) (see Figures 6B, C ).

In faf-stress, net imbalance is the result of an intensification of the srADV cooling compensated by vDIFF warming in both hemispheres ( Figures 5C, D ). In the NH, the negative imbalance is given by the srADV in the shallow ocean, but in the SH, the positive imbalance is produced by enhanced vDIFF in the deep ocean. In contrast, in faf-heat, total warming imbalances are set primarily by strengthened vDIFF ( Figures 5E, F ). In both hemispheres, the source of heating is the strengthened vDIFF warming and reduced srADV cooling. However, while the shallow ocean in the NH holds nearly 70% of the total net warming imbalance, it only accounts for approximately 10% in the SH ( Table 2 ).

According to the volume-integrated budgets, in faf-stress, net imbalance is the result of intensified srADV and vDIFF balance, whereas in faf-heat, warming in both hemispheres is set by strengthened vDIFF and weakened srADV mainly in the deep ocean.

Both in faf-stress and faf-heat, imbalances are located nearly at the same vertical level and, given the different vertical penetration of the cell in the two hemispheres, the imbalance is promoted by the shallow overturn in the NH and the bottom circulation in the SH ( Figures 6D–I ). This reasserts the importance of separating the ocean into dynamical regimes, rather than using vertical budget integrations based solely on arbitrary selection of boundaries between layers. Our results highlight a local redistribution of heat from the surface to the deep ocean, since imbalances in the unventilated tropical Pacific appear to be linked to subtropical surface anomalies, advected to the deep ocean along mean isopycnal pathways.

Because ocean heat transport is linked to water mass transformation, it is useful to estimate the relative contribution of the different water masses in the overall heat transport (Talley, 2003). Meridional transports are computed for the shallow wind-driven overturning cells and the zonally integrated transport at different density classes, with each water mass classified according to its ventilation level. Heat budget calculation is presented for each experiment. The meridional transport inside the cell is divided into surface poleward transport, WBC contribution, and the interior equatorward transport, computed both at 24°S and 24°N. Balance inside the overturning cell excludes transports within deeper levels of the WBC, which flows outside of the gyre towards the subpolar region (Talley, 2003). Estimations consider mass balance adjusted to zero, but there is a positive net heat transport for this mass-balanced shallow overturn ( Tables 3 , 4 ).

The zonally integrated mass and heat transport within different water masses is shown in Figure 4 for the Control. In this case, the transport associated with the upper layer represents the residual between the total poleward WBC flow and the equatorward interior flow. All remaining layers refer to meridional mass transports below the cell. In the Control simulation, in the NH, the shallow overturning cell recirculates nearly 32 Sv and 2.34 PW according to the interior flow ( Figure 4A ). The net heat transport by the shallow overturning contributes nearly 65% of the total heat transport. Poleward transport by the cell is given not only by the surface Ekman flow but also by the WBC in both hemispheres ( ∼ 60%), a contribution related to the geostrophic circulation within the subtropical gyre.

In contrast, according to the zonally integrated transports, most of the surface northward transport returns within the cell ( > 65%), while the rest returns as intermediate water mixed with upwelled deep water. Transfer between water mass classes, assuming the transfer occurs preferentially between most adjacent layers, is schematically represented with vertical arrows in Figure 4B . In fact, in terms of the equatorward intermediate water, this is mostly composed of surface water with a small contribution from upwelled deep water.

In the SH, the shallow overturning cell recirculates nearly 16~Sv and 1.7~PW according to the interior flow. As in the NH, WBC transport is responsible for nearly 50\% of both poleward mass and heat transports. In the zonally-integrated mass and heat transports, the total and shallow overturning transports are fully accounted for by the surface layer. Most of the surface poleward transport returns within the cell as upper waters ( ∼ 70%), while the rest returns as converted deep waters mixed with poleward downwelling intermediate waters, and the contribution associated with intermediate waters is rather secondary (20%). Deep waters promote a diffusive transfer of heat towards the bottom of 0.08 PW, which nearly represents 25% of the heat convergence within deep waters. The shallow ocean determines nearly the entire meridional mass and heat transports in the NH ( > 90%) and most of the transports in the case of the SH ( ∼ 70%). Moreover, the WBC plays an important role in the overall poleward shallow overturn transport of both mass and heat ( > 50%).

Next, we focus on the response under surface flux perturbations. The relative change of each meridional mass transport component is computed with respect to the equilibrium state and shown in brackets (for details, see Tables 3 , 4 ). In faf-stress, the NH shallow overturning cell reveals an overall weakening ( Figure 7A ). The weakened signal is explained in terms of a weakened interior equatorward transport (−11%) and WBC transport (−20%). In contrast, the production of dense water by the cell is significantly increased given by the increased surface poleward transport. This is in accordance with the strengthened equatorward intermediate and deep water mass production and also the weakened poleward deep water mass production (−20%). These changes can be explained in terms of reduced mid-latitude subduction, according to the zonally integrated mass and heat transports ( Figure 7B ). Although faf-stress imposes relatively weak reduction in both tropical and subtropical wind stress ( Figure 2 ), these changes are able to promote weakening of the mid-latitude subduction and slowing down of the overturning cell. Moreover, circulation changes are in agreement with observed net heat imbalances observed in both the shallow and deep ocean, by increased srADV cooling and vDIFF warming. While strengthened circulation produces cooling of the shallow ocean, intensified intermediate water production produces warming of the deep ocean ( Figure 8A ).

The SH shallow overturning cell shows an overall strengthening. This is based on the strengthened surface mass and heat transport and the strengthened interior equatorward transport. This is in line with the observed strengthening wind stress over the SH Pacific Ocean (see Figure 2A ). However, intensification of the overturn seems to reduce the poleward transport outside of the cell, inhibiting the formation of dense water. Moreover, given the increased mid-latitude subduction of surface water by the cell, the production of dense water is reduced, as suggested by the reduced WBC residual transport and the reduced intermediate and deep water masses transport below the cell. In contrast to the NH, wind stress anomalies promote an increase in mid-latitude subduction and poleward transport outside of the cell. The strengthening of the shallow overturning cell and the significantly decreased formation of dense water at higher latitudes is in agreement with observed warming imbalances. In both cases, the shallow and deep oceans warm by increased srADV cooling and vDIFF warming ( Figure 8A ).

Results in faf-heat are shown in Figure 9 . Based on the large surface poleward flow strengthening in mass and heat transports, the NH shallow overturning cell reveals an overall strengthening (approximately 20%) compared to the control run ( Figure 9A ). This is achieved by the surface layer alone since the WBC is weakened (−6%). In fact, according to the zonally integrated transports ( Figure 9B ), both poleward surface transport and equatorward upper transport are significantly intensified (20% and ∼ 50%, respectively). Hence, despite a cell circulation strengthening, the net poleward surface transport outside of the cell is weakened, as well as the subsequent downwelling of this water into deeper levels. This is shown by a significant reduction of the residual transport (−35%) and the equatorward intermediate water mass production (−31%).

In contrast, the SH shallow overturning cell reveals a non-significant change in circulation. Once we consider the zonally integrated transport, while the poleward surface transport reveals no significant change, the upper equatorward transport shows a significant increase (25%) compared to the control run. This is explained in terms of the significant reduction of residual transport by the WBC (−55%). This reduction in mass export by the cell towards deeper levels seems to be responsible for the observed deep water production weakening (−35%) and the poleward intermediate-to-deep water strengthening by recirculation of the subducted water. The strengthening of the NH overturning cell, accompanied by the significant weakening of dense water production in both hemispheres ( ∼ 50%), is in agreement with the shallow and deep ocean net warming (by reduced srADV cooling and increased vDIFF warming). The strengthened cell circulation seems to promote the reduction of the transport of mass and heat out of the cell, thus increasing the retention of warm water by the cell and potentially explaining the increased heat content ( Figure 8B ). The same intensification of cell circulation seems to promote the reduction of dense water production, thus weakening deep-to-abyssal branches of the overturning circulation (weakened srADV cooling). In the SH, a similar mechanism (reduction of dense water production) may explain the overall warming of the deep ocean. In other words, by reducing dense water production, ventilation of deep ocean is reduced, increasing the retention time for dense water and the heat content as well.

The shallow overturning cell presents significant changes in faf-heat. However, in contrast to faf-stress, where the observed changes are explained in terms of the anomalous wind-driven circulation, in faf-heat, changes in the shallow overturn can be explained in terms of isopycnal dependency of the meridional flow across the cell. The strengthening of the interior pycnocline transport is associated with the intensification of the zonal pycnocline tilt, particularly in the WBC region, due to the coarse-resolution setup (Huber and Zanna, 2017), and by reducing such a strong pycnocline dependency, circulation changes in faf-heat might be significantly reduced. Further analyses are required to assess the impact of the model resolution, parametrizations, and the role of ocean–atmosphere coupled feedbacks.