Surging of Global Surface Temperature due to Decadal Legacy of Ocean Heat Uptake

BABLU SINHA

National Oceanography Centre, Southampton, United Kingdom

FLORIAN SÉVELLEC

Laboratoire d'Océanographie Physique et Spatiale, Univ Brest CNRS IRD Ifremer, Brest, France, and Ocean and Earth Science, University of Southampton, Southampton, United Kingdom

JON ROBSON

National Centre for Atmospheric Science, Department of Meteorology, University of Reading, Reading, United Kingdom

A. J. GEORGE NURSER

National Oceanography Centre, Southampton, United Kingdom

(Manuscript received 26 November 2019, in final form 4 April 2020)

ABSTRACT

Global surface warming since 1850 has consisted of a series of slowdowns (hiatus) followed by surges. Knowledge of a mechanism to explain how this occurs would aid development and testing of interannual to decadal climate forecasts. In this paper a global climate model is forced to adopt an ocean state corresponding to a hiatus [with negative interdecadal Pacific oscillation (IPO) and other surface features typical of a hiatus] by artificially increasing the background diffusivity for a decade before restoring it to its normal value and allowing the model to evolve freely. This causes the model to develop a decadal surge that overshoots equilibrium (resulting in a positive IPO state), leaving behind a modified, warmer climate for decades. Water-mass transformation diagnostics indicate that the heat budget of the tropical Pacific Ocean is a balance between large opposite-signed terms: surface heating/cooling resulting from air-sea heat flux is balanced by vertical mixing and ocean heat transport divergence. During the artificial hiatus, excess heat becomes trapped just above the thermocline and there is a weak vertical thermal gradient (due to the high artificial background mixing). When the hiatus is terminated, by returning the background diffusivity to normal, the thermal gradient strengthens to prehiatus values so that the mixing (diffusivity × thermal gradient) remains roughly constant. However, since the base layer just above the thermocline remains anomalously warm, this implies a warming of the entire water column above the trapped heat, which results in a surge followed by a prolonged period of elevated surface temperatures.

1. Introduction

The estimated 0.85-K global surface warming since 1850 (Alexander et al. 2013, p. 5) has been realized as a succession of periods of stronger followed by weaker warming (Risbey 2015). Understanding the causes of these hiatus and surge events is of key importance for interannual-to-decadal prediction (Guemas et al. 2013; Meehl et al. 2014;

Sévellec and Drijfhout 2018) and for constraining uncertainty in long-term climate projections.

A large uncertainty remains concerning the detailed processes responsible for surge and hiatus periods, with debate centered around the relative roles of external forcing (e.g., anthropogenic greenhouse gas and aerosol emission variations, solar radiation changes, and volcanic

Denotes content that is immediately available upon publication as open access.

Publisher's Note: This article was revised on 8 March 2021 to include the CCBY license that was missing when originally published.

Corresponding author address: Bablu Sinha, National Oceanography Centre, European Way, Southampton SO14 3ZH, United Kingdom. E-mail: bs@noc.ac.uk

<figure><img src="image_1.png" /></figure>

This article is licensed under a Creative Commons Attribution 4.0 license (http://creativecommons.org/licenses/by/4.0/).

DOI: 10.1175/JCLI-D-19-0874.1
eruptions) versus internal variability (Smith et al. 2015) and the dominance or otherwise of the Pacific Ocean over the Atlantic Ocean (England et al. 2014; Chen and Tung 2014; Drijfhout et al. 2014; Meehl et al. 2011).

During the recent global surface warming hiatus of the 2000s Smith et al. (2015) estimated a reduction in net top-of-atmosphere (TOA) radiation of −0.31 ± 0.21 W m⁻² between 1999 and 2005 but could not identify whether this was due to internal processes or external drivers. Even if the source of the TOA reduction were disentangled, Song et al. (2014) demonstrate using the CMIP5 multimodel ensemble that on decadal time scales TOA imbalance is only very weakly related to surface air temperature (SAT) because natural variability is superimposed on global warming (see also Xie et al. 2016). Also, Hedemann et al. (2017) point out that it is the difference between TOA imbalance and ocean heat transfer between the surface (~100 m) layer and the deep ocean that matters for a hiatus, and this could be too small (~0.08 W m⁻²) to resolve with the current observing system.

Leaving aside the ultimate causes of hiatus and surge events (external forcing vs internal variability), the hiatus of the 2000s coincided with a negative phase of the interdecadal Pacific oscillation (IPO), consistent with the results of Meehl et al. (2013), who found the IPO to be one of three main processes responsible for hiatus and surge behavior in the CCSM4 coupled climate model, with the others being Antarctic Bottom Water formation and the Atlantic meridional overturning circulation (AMOC). Observations demonstrated that the hiatus was associated with anomalously cold Pacific sea surface temperatures (SSTs) and warm thermocline anomalies (Nieves et al. 2015), due to strengthened subtropical circulation cells forced by strengthened trade winds (Farneti et al. 2014; England et al. 2014). That the IPO, especially the equatorial Pacific, played a key role in the hiatus was demonstrated by England et al. (2014), Watanabe et al. (2014), and Kosaka and Xie (2013) using pacemaker experiments in which either tropical Pacific winds or SST were forced to follow observations in coupled model simulations, which then reproduced the hiatus. However, increased ocean heat uptake was also seen in the deep North Atlantic and Southern Oceans during this period (Balmaseda et al. 2013; Chen and Tung 2014; Drijfhout et al. 2014) and the role of these regions, especially the North Atlantic where the AMOC is strongly associated with heat content changes (Moat et al. 2019; Williams et al. 2015) and the surface atmospheric temperature is highly variable (Sévellec et al. 2016), remains unclear.

More recent work has documented the long residence time of heat taken up by the ocean and stored below the mixed layer during the recent hiatus of the 2000s. Maher et al. (2018), using an ocean model forced by negative IPO-like surface atmospheric conditions, found drawdown of

heat into the main thermocline (~300-m depth) and surface cooling in the Pacific. The extra subsurface heat was advected by subtropical cells, strengthened as a result of enhanced trade winds, into the Indian Ocean where it remained buried in the thermocline, as seen earlier in analyses of hiatus periods in the CMIP5 multimodel ensemble (Lee et al. 2015). A similar result was found in the pacemaker experiments of Gastineau et al. (2019), although the spatial redistribution of heat showed differences; for example, much of the additional heat entering the Indian Ocean from the Pacific during the hiatus was returned to the atmosphere via surface fluxes, or advected southward into the Southern Ocean. The spreading of this additional heat from the Pacific to the Indian Ocean, seen in observations as well as models (Nieves et al. 2015; Liu et al. 2016; Lee et al. 2015), and its partial reemergence after the hiatus suggest a substantial legacy of, and potentially an active role for, the ocean heat uptake during the hiatus period. This was illustrated in the study of Maher et al. (2018) by the response to a subsequent reversal of IPO phase, in which the sequestered heat was not re-extracted back to the atmosphere.

The consensus then is that the IPO determines the pattern of heat drawdown above the tropical thermocline and subsequent horizontal redistribution via associated changes in the strength of the shallow subtropical overturning cells.

The significance of the increased heat uptake seen in the deep Southern Ocean and North Atlantic Ocean (Meehl et al. 2013; Drijfhout et al. 2014; Chen and Tung 2014) remains ambiguous. Indeed, Drijfhout et al. (2014) suggest that the enhanced Southern Ocean heat uptake may be an ongoing multidecadal process linked to external forcing of the southern annular mode rather than natural variability. The role of the North Atlantic remains a thorny issue with the mechanism relating AMOC changes to ocean heat uptake still debated (Chen and Tung 2014 vs Drijfhout et al. 2014).

The central role of ocean heat uptake in variations of global mean surface temperature and the TOA radiation balance motivates us to explore the effects of its variability on decadal time scales. Tanaka et al. (2012) show that changes in ocean mixing can affect Pacific decadal variability, a key player in the hiatus of the early 2000s (Meehl et al. 2011; England et al. 2014; and others; see above). Hence here we explore the effects of transient enhancement of ocean mixing, and the subsequent climate adjustment, on ocean heat uptake and global mean surface temperature.

## 2. Method

### a. Climate model and experiment design

We use the HadCM3 climate model (Gordon et al. 2000), which maintains a stable climate simulation without
<table><caption>TABLE 1. Summary of annual mean energy flows (W m⁻²) in the HadCM3 control simulation and comparison with observational estimates. Note that, as defined by Wild et al. (2013), the TOA outgoing solar comprises all reflected solar radiation including that reflected from the surface.</caption><thead><tr><th></th><th>HadCM3 control simulation</th><th>Wild et al. (2013) observation-based estimates</th></tr></thead><tbody><tr><td>TOA incident solar</td><td>341.4</td><td>340 (340–341)</td></tr><tr><td>TOA outgoing solar</td><td>101.1</td><td>100 (96–100)</td></tr><tr><td>Atmospheric solar absorption</td><td>76.1</td><td>79 (74–91)</td></tr><tr><td>Surface incident solar</td><td>188.5</td><td>185 (179–189)</td></tr><tr><td>Surface reflected solar</td><td>24.3</td><td>24 (22–26)</td></tr><tr><td>Surface solar absorption</td><td>164.2</td><td>161 (154–166)</td></tr><tr><td>Outgoing longwave</td><td>240.5</td><td>239 (236–242)</td></tr><tr><td>Surface upwelling longwave</td><td>394.6</td><td>398 (394–400)</td></tr><tr><td>Surface downwelling longwave</td><td>332.0</td><td>342 (338–348)</td></tr><tr><td>Surface sensible heat</td><td>17.0</td><td>20 (15–25)</td></tr><tr><td>Surface latent heat</td><td>84.1</td><td>84 (70–85)</td></tr></tbody></table>

flux adjustments despite relatively coarse grid resolution (1.25° × 1.25° ocean and 3.75° × 2.5° atmosphere). The HadCM3 preindustrial control simulation provides a reasonable simulation of the observed energy flows between components of the Earth system (Table 1). We spin experiments off the control by instantaneously doubling the background vertical diffusivity at every gridpoint and integrating the model for one decade before instantaneously returning the diffusivity to its control value and continuing the integration for one or more additional decades. To address internal variability, we create a four-member ensemble by starting from different points of a 140-yr portion of the 2000-yr control simulation. All ensemble members begin on the same day of the year, 1 December.

The in situ vertical diffusivity for tracers consists of a time-independent (background) part, which only varies with depth, and a time-dependent part, which depends on the stratification and velocity shear present and varies spatially and temporally (Gordon et al. 2000, see appendix 2 therein). The background part increases with depth from ~10⁻⁵ m² s⁻¹ at the surface to ~15 × 10⁻⁵ m² s⁻¹ in the abyss (Gordon et al. 2000). We instantaneously double this background part at the beginning of the ensemble experiment and set it back to normal at the beginning of the second decade. Ocean models also include substantial numerical mixing up to several times the explicit mixing (Megann 2018), and hence the mixing perturbation we apply only modifies the actual in situ mixing coefficient by perhaps 10%–20%, not far removed from its observational uncertainty.

Ensemble mean quantities are evaluated for significance against the control (all variables are first filtered using a 5-yr running mean unless otherwise stated) using a one-tailed *t* test, taking into account the number *n* of independent degrees of freedom in the 140-yr control simulation (estimated to be *n* = 26 because of autocorrelation in the time series).

## b. Water-mass transformation framework

Since we modify the diffusion coefficient and hence the mixing in our experiments, a natural framework for diagnosis of the results is water-mass transformation analysis (Walin 1982). The strength of this framework is that it diagnoses mixing from all sources, including numerical mixing, which cannot easily be quantified by alternative methods (Megann 2018; Lee et al. 2002). In addition, the water-mass transformation framework can be applied directly to observations and was developed in that context (Walin 1982).

For convenience of exposition, we describe this framework for a zonally integrated domain; however, the extension to three dimensions is important for the real ocean, and hence we also discuss how to interpret the equations in a 3D context.

Referring to Fig. 1a, we begin with the “transformation”—that is, the rate $G(\theta, \phi)$ at which water crosses some isotherm with potential temperature $\theta$ integrated zonally over all longitudes $\lambda$ and meridionally north of some latitude $\phi$. This is related by mass (or volume in a Boussinesq fluid) conservation to the meridional influx across $\phi$,

$$ \psi(\theta, \phi) = \iint_{\lambda,z:\theta'(\lambda,\phi,z)<\theta} \nu R \cos\phi \, dz \, d\lambda, \quad (1) $$

of water colder than $\theta$ (where $\nu$ is meridional velocity and $R$ is the radius of Earth) and the rate of change of the volume $V(\theta, \phi)$ north of $\phi$ and colder than $\theta$ by

$$ \frac{\partial V(\theta, \phi)}{\partial t} = \psi(\theta, \phi) - G(\theta, \phi). \quad (2) $$

This transformation $G(\theta, \phi)$ represents water crossing the isotherm from cool to warm, and thus requires heat input; to link $G(\theta, \phi)$ to heat inputs we consider (Fig. 1b)
<figure><img src="image_2.png" /><figcaption>FIG. 1. Derivation of the water-mass transformation equations: (a) volume conservation for waters with potential temperature less than or equal to $\theta$, (b) heat conservation for the same region, and (c) water-mass formation rate between isopycnals $\theta$ and $\theta + \delta\theta$ from surface flux forcing. Also shown are illustrations of 3D water-mass transformation diagnostics: (d) ocean temperature in the western tropical Pacific (contours). Surface waters below 22°C are shaded brown. Surface and subsurface waters with temperatures between 22° and 23°C are shaded purple. Colored arrows indicate water-mass transformation by air-sea heat exchange (yellow) and ocean circulation/ overturning (green) across the 22° and 23°C isotherms. Solid and dashed black arrows indicate diapycnal mixing across the 23° and 22°C isotherms, respectively. (e) As in (d), but for a domain extended slightly to the south. The corresponding diapycnal mixing in the southward extension to the domain is indicated by hollow arrows.</figcaption></figure>

the (Boussinesq) equation for the heat content of water colder than $\theta$ (e.g., Nurser et al. 1999):

$$ \frac{\partial \mathbf{H}}{\partial t} = -D_{\text{diff}} - D_{\text{surf}} + A - c_p \rho_0 G \theta . \quad (3) $$

Here

$$ \mathbf{H} = \rho_0 c_p \iint_{\lambda, \phi', z: \theta'(\lambda, \phi', z) < \theta, \phi' < \phi} \theta' R \cos\phi' \, d\lambda \, d\phi \, dz \quad (4) $$

is the heat contained within seawater with potential temperature less than or equal to $\theta$ integrated everywhere north of a given latitude $\phi$ and $\partial \mathbf{H} / \partial t$ is its time derivative. In Eq. (3), $D_{\text{diff}}$ is the (cold to warm) diffusive temperature flux across the $\theta$ isotherm north of latitude $\phi$ from mixing, and

$$ D_{\text{surf}} = - \iint_{\lambda, \phi, \theta_s(\lambda, \phi) \le \theta} I(\lambda, \theta, z_{\theta}) R \cos\phi \, d\lambda \, d\phi \\ - \iint_{\lambda, \phi, \theta_s(\lambda, \phi) \le \theta} [I_s - Q_{\text{net}}(\lambda, \phi)] R \cos\phi \, d\lambda \, d\phi \quad (5) $$

represents the heat fluxes upward and out (first term) through the $\theta$ isotherm from solar irradiance $I$ (sign convention is positive downward) where the isotherm is close to the ocean surface and upward through the surface boundary where the $\theta$ isotherm has outcropped [second term; the net turbulent air-sea flux $Q_{\text{NET}} = Q_{\text{latent}} + Q_{\text{sensible}} + Q_{\text{LW}}^{\text{net}}$ is made up of (positive upward) latent, sensible, and net longwave (LW) radiative heat fluxes and the surface irradiance $I_s$ is positive downward; note that $\theta_s$ denotes the SST]. The term

$$ A(\theta, \phi) = \rho_0 c_p \iint_{\lambda, z: \theta'(\lambda, \phi, z) < \theta} \theta' v R \cos\phi \, dz \, d\lambda \quad (6) $$

represents the advective heat transport of water colder than $\theta$ northward across latitude $\phi$, and $c_p \rho_0 G \theta$ represents the advective heat flux up and out across the $\theta$ isotherm.

Noting that the transport of water with $\theta < \theta' < \theta + \delta\theta$ is

$$ \delta\psi = \iint_{\lambda, z: \theta < \theta'(\lambda, \phi, z) < \theta + \delta\theta} v R \cos\phi \, dz \, d\lambda \approx \frac{\partial\psi}{\partial\theta} \delta\theta, $$
we can re-express (integrating by parts for the last result)

$$A(\theta, \phi) = \rho_0 c_p \int_{\theta}^{\theta} \frac{\partial \psi}{\partial \theta} d\theta' = \rho_0 c_p \psi \theta - \rho_0 c_p \int_{\theta}^{\theta} \psi d\theta' \quad (7a)$$

and similarly obtain

$$\mathbb{H}(\theta, \phi) = \rho_0 c_p \int_{\theta}^{\theta} \frac{\partial V}{\partial \theta} d\theta' = \rho_0 c_p V \theta - \rho_0 c_p \int_{\theta}^{\theta} V d\theta'. \quad (7b)$$

Eliminating $\mathbb{H}$ and $A$ from Eq. (3) using Eqs. (7a) and (7b) and removing the $G\theta$, $\psi\theta$, and $(\partial V/\partial t)\theta$ terms using Eq. (2) gives

$$-\rho_0 c_p \int_{\theta}^{\theta} \frac{\partial V}{\partial t} d\theta' = -D_{\text{diff}} - D_{\text{surf}} - \rho_0 c_p \int_{\theta}^{\theta} \psi d\theta'$$

Differentiating this by $\theta$ and again using Eq. (2) to replace $-\partial V/\partial t + \psi$ by $G$ gives the key result

$$\rho_0 c_p G = -\frac{\partial D_{\text{diff}}}{\partial \theta} - \frac{\partial D_{\text{surf}}}{\partial \theta}. \quad (8)$$

Multiplying Eq. (8) by $\delta\theta$ makes clear that Eq. (8) represents the heat budget of the volume between the $\theta$ and $\theta + \delta\theta$ isotherms, relating the transformation (volume flux from cold to warm) across the isothermal surfaces to the heat acquired from the convergence of air-sea fluxes and mixing. It is customary to introduce $F$, the transformation driven by air-sea fluxes $\phi$,

$$\rho_0 c_p F = -\frac{\partial D_{\text{surf}}}{\partial \theta}; \quad (9)$$

it then follows (Fig. 1c) that

$$\begin{aligned} \rho_0 c_p F \delta\theta &= -\frac{\partial D_{\text{surf}}}{\partial \theta} \delta\theta \\ &= \iint_{\lambda, \phi; \phi' > \phi, \theta_s (\lambda, \phi) > \theta + \delta\theta} I(\lambda, \phi', z_{\theta+\delta\theta}) R \cos\phi \, d\lambda \, d\phi' \\ &\quad - \iint_{\lambda, \phi; \phi' > \phi, \theta_s (\lambda, \phi) > \theta} I(\lambda, \phi', z_{\theta}) R \cos\phi \, d\lambda \, d\phi' \\ &\quad + \iint_{\lambda, \phi'; \phi' > \phi, \theta < \theta_s (\lambda, \phi) \le \theta + \delta\theta} [I_s - Q_{\text{net}}(\lambda, \phi')] R \cos\phi \, d\lambda \, d\phi', \end{aligned} \quad (10)$$

where the three terms on the RHS of Eq. (10) represent (i) the solar irradiation down through the $\theta + \delta\theta$ isotherm (a heat source for the layer $\theta < \theta' < \theta + \delta\theta$), (ii) the solar irradiation down through the $\theta$ isotherm (a heat sink), and (iii) the surface heat flux down through the outcrop $\theta < \theta_s < \theta + \delta\theta$. In terms of $F$, the water-mass transformation $G$ is then given as

$$\rho_0 c_p G = \rho_0 c_p F - \frac{\partial D_{\text{diff}}}{\partial \theta}. \quad (11)$$

While it is, in principle, possible to evaluate the diffusive fluxes, they are difficult to evaluate in observations and many models, so they are often in practice estimated as a residual between $G$, [which can be found by Eq. (1) from changes in volume $\partial V/\partial t$ and transport $\psi$] and $F$, which is known from the surface fluxes.

We can use Eq. (11) to eliminate $G$ in Eq. (2), arriving at

$$\frac{\partial V(\theta, \phi)}{\partial t} = \psi(\theta, \phi) - F(\theta, \phi) + \frac{1}{c_p \rho_0} \frac{\partial D_{\text{diff}}}{\partial \theta}. \quad (12)$$

Now $(\partial V/\partial \theta)\delta\theta$ represents the volume of water with temperature between $\theta$ and $\theta + \delta\theta$, so differentiating Eq. (12) by $\theta$ gives

$$\frac{\partial^2 V(\theta, \phi)}{\partial t \partial \theta} = \frac{\partial \psi(\theta, \phi)}{\partial \theta} - \frac{\partial F(\theta, \phi)}{\partial \theta} + \frac{1}{c_p \rho_0} \frac{\partial^2 D_{\text{diff}}}{\partial \theta^2}, \quad (13)$$

which describes the increase in volume of water of temperature $\theta$ in terms of its southern influx $[\partial\psi(\theta, \phi)/\partial\theta]\delta\theta$, and the convergence of flow across the $\theta$ and $\theta + \delta\theta$ isotherms.

To bring out the physical meaning in three dimensions of the terms in Eqs. (12) and (13), Fig. 1d shows simulated temperature in one particular December in a region of the Pacific dominated by the shallow northern subtropical overturning cell (19.5°–30°N, 140°E–120°W and 0–300-m depth). The equatorial warm pool is visible to the southwest ($\theta > 27^\circ C$). From here, isotherms slope up, eastward and northward. Water between 22° and 23°C (purple shading) outcrops at the surface in a crescent oriented northwest to southeast. In 3D, the water mass resembles a bowl shape truncated by the southern boundary of the region where the water mass persists as a thin near-horizontal feature at about 80-m depth at 160°E extending eastward to about 160°W, where it rises vertically to intersect the surface.

Applying Eq. (12), we integrate the surface heat flux across all surface waters to the north with $23^\circ < \theta \le 24^\circ C$ (purple shaded surface region extending to the North Pole). This is the larger yellow arrow ($F_{23^\circ C, 19.5^\circ N}$), a loss of volume, which would be compensated in steady state by flow of warmer water across the 23°C isotherm (solid black arrows). With the heat flux integrated back to the North Pole, we cannot say anything about the latitudinal distribution of this diapycnal flow. Now, restricting ourselves to waters to the north with $22^\circ < \theta \le 23^\circ C$ we obtain a smaller value (smaller yellow arrow $= F_{22^\circ C, 19.5^\circ N}$), which must be balanced by a diapycnal flow of water across the
22°C isotherm (dashed black arrows). Hence the *difference* [i.e., applying Eq. (13)] between these two transformation rates (known as the formation rate) results in a tendency for the volume of the purple region to increase.

The formation rate may be partially balanced by ocean circulation/overturning (green arrows). The upper green arrow represents the volume flux out of the domain for water with $\theta \le 23^\circ C$. The lower green arrow represents the volume flux out of the domain for water with $\theta \le 22^\circ C$. Hence, the volume flux out of the purple region is equal to the difference between the green arrows. Both arrows are southward because in this part of the ocean the circulation consists of a strong northward flow in a thin layer at the surface (higher $\theta$) and a weaker southward flow in a thicker, deeper layer (lower $\theta$). As we integrate up from the lowest temperatures, the streamfunction itself is negative, but its vertical gradient is negative at depth/low $\theta$ and positive at the surface/ high $\theta$.

Similar reasoning leads us to conclude that changes to the volume and mixing in the purple region are respectively given by the difference in $\partial V/\partial t \partial D_{\text{diff}}/\partial \theta$ with respect to temperature.

To localize these quantities spatially, consider two nearby latitudes. Figure 1e is identical to Fig. 1d except an additional slab of ocean to 18°N is added to the south, and the yellow arrows represent the surface heat flux integrated from the North Pole down to the new southern boundary ($F_{23^\circ C,18^\circ N}$; $F_{22^\circ C,18^\circ N}$). Evidently, the diapycnal fluxes in the region between 18° and 19.5°N (hollow arrows) are given by the difference between values of $F$ at the two different latitudes. Similar considerations apply to the overturning (green arrows) and to the time derivative and mixing terms.

Note the resemblance of Eq. (12) to the equation for the evolution of ocean temperature (neglecting horizontal diffusion):

$$ \frac{\partial T}{\partial t} = \frac{\partial}{\partial z} \left( K_v \frac{\partial T}{\partial z} \right) + \frac{I}{\rho_0 C_p} - \nabla \cdot (\mathbf{u}T). \quad (14) $$

Here $K_v$ is the diffusion coefficient, $\mathbf{u}$ is the three-dimensional ocean velocity, and $I$ is the solar irradiance encountered earlier. The correspondence between the terms in Eqs. (12) and (14) becomes more obvious when we consider the surface boundary condition

$$ K_v \left. \frac{\partial \theta}{\partial z} \right|_{z=0} = \frac{Q_{\text{NET}} - I_s}{\rho_0 C_p}, \quad (15) $$

where $Q_{\text{NET}}$ (the net surface heat flux) and $I_s$ (the surface irradiance) were introduced earlier. Recall that $K_v$ consists of a background value, which increases with depth,

<figure><img src="image_3.png" /><figcaption>FIG. 2. Time series of ensemble mean global mean area-weighted temperature anomaly (K) with respect to the 140-yr control. Annual means are plotted, and a 4-yr boxcar filter is applied to smooth the data. Color shading indicates values that are significantly different from the control simulation at the 95% confidence level. The dotted vertical line separates the hiatus period (years 1–10; doubled vertical diffusivity) from the surge period (years 11–20; standard vertical diffusivity).</figcaption></figure>

and a time-varying part, which depends on stratification and velocity shear.

In the surface mixed layer $K_v$ becomes large while $\partial \theta/\partial z$ is small, and the surface flux $Q_{\text{NET}} - I_s$ is effectively spread over a vertical distance on the order of the mixed layer depth. We shall use this property to develop a conceptual model in section 3f.

## 3. Results

### a. Ocean temperature response

The global temperature response in the perturbed ensemble experiment is shown in Fig. 2, which is a time-depth plot of horizontal averaged ensemble mean ocean temperature anomaly relative to the mean of the 140-yr control. In the first decade after the diffusivity is increased, near-surface ocean temperatures (0–100 m) decrease by ~0.15 K, reaching a minimum at year 8 [2 years before the diffusivity is returned to its original value; this is seen in the unfiltered data (not shown) and is not an artifact of filtering]. Simultaneously, the subsurface ocean (below 100 m) warms, with maximum warming at ~200 m around year 9. Thus, the change in mixing results in transfer of heat from the surface layer to the thermocline, reminiscent of the observed hiatus of the 2000s. When the diffusivity is returned to its original value at the beginning of year 11, the surface ocean warms, with peak surface temperature occurring at year 15. In the subsurface ocean, there is slight cooling in the 200–300-m depth interval, but below this the temperature remains roughly constant for the remaining 10 years of the experiment. Hereinafter, we refer to years 1–10 of the ensemble
<figure><img src="image_4.png" /><figcaption>FIG. 3. Ensemble mean zonal mean temperature anomaly (K) with respect to the 140-yr control mean during (a) the hiatus period (years 6–10) and (b) the surge period (years 11–15). Thick black contours denote values that are significantly different from the control simulation at the 95% confidence level. Thin black contours denote the control mean temperature (°C).</figcaption></figure>

experiment as the hiatus period and years 11–20 as the surge period, although strictly this should be defined in terms of anomalous global mean temperature tendency. The heat taken up during the simulated hiatus period remains in the ocean for at least a decade after the diffusivity has been restored to its original value. Although anomalously cool during the hiatus, global mean SST becomes anomalously warm over the surge period.

Ensemble mean zonal average ocean temperature anomalies during years 6–10 (Fig. 3a) indicate cooling in the tropics to subtropics (30°S–30°N) from the surface to about 100 m and warming in the permanent thermocline throughout the tropics and midlatitudes, but particularly marked in the tropics (20°S–20°N; 100–400 m). This pattern of surface cooling and deep warming is strongly reminiscent of the observed hiatus of the 2000s (e.g., Drijfhout et al. 2014).

During years 11–15 there is warming throughout the tropics and midlatitudes down to 1000-m depth (Fig. 3b). However, weak cooling of the high latitudes that was present during years 1–10 remains, particularly in the Northern Hemisphere to depths of 800 m, but is counterbalanced by a

slight warming from 800 to 2500 m (not shown). Once again, the heat taken up by the ocean into the main thermocline during the hiatus period persists into the surge period.

Figure 4a shows the ensemble mean SST anomaly averaged over years 6–10. Surface cooling is particularly marked in the tropics. The western boundary current regions, by contrast, exhibit warming, albeit over a smaller area. A negative IPO-like pattern is apparent with a wedge-shaped region of cooling in the tropical Pacific and warming in the western North Pacific and the southeastern Pacific.

Following restoration of the diffusivity to its original value, the surface warms, with Pacific SST showing a strongly positive IPO pattern averaged over years 11–15 (Fig. 4b). Marked warming remains over the western boundary currents and their extensions in the Northern Hemisphere, and the Indian Ocean shows substantial surface warming. The Southern Ocean, southwestern Pacific, and South Atlantic show a weak cooling. Areas of statistically significant warming during the surge are less extensive than areas of cooling during the hiatus. However, a clear global-scale pattern is evident.

The SST anomaly patterns can be compared with the 0–1000-m ocean heat content anomaly. During years 6–10 there are substantial rises in ocean heat content in the tropics and subtropics (40°S–40°N), which remain largely in place in years 11–15; that is, the uptake of heat during the hiatus is not reversed during the surge (Figs. 4c,d). Unlike the hiatus, the surge SST anomaly pattern in years 11–15 echoes that of the ocean heat content anomaly over much of the ocean: in the tropical and South Pacific, the Indian Ocean, and the North Atlantic, but not in the South Atlantic, the North Pacific, or the southwest Pacific.

With regard to the ocean surface heat uptake anomaly (annual mean anomaly of net surface heat flux into the ocean; Figs. 4e and 4f), there are large increases in years 6–10 throughout the tropics and subtropics of all ocean basins, and in some subpolar regions. A zonally banded structure is seen, seemingly related to oceanic fronts and mode water formation regions with strong anomalies over the ITCZ and the western boundary currents (particularly the Gulf Stream) and adjacent to Antarctica. When the diffusivity is returned to its standard value in years 11–15, the ocean heat uptake anomaly reduces substantially across the globe but remains positive in many regions, including the tropical Pacific and over the Gulf Stream. The North Atlantic subpolar gyre by contrast shows a strong reduction in ocean heat uptake during the surge period.

### b. Energy budget

Figure 5 indicates the impact of the perturbation on Earth's energy flows (Wild et al. 2013). There is no effect on incident solar radiation at the TOA and hence this is not plotted. Figure 5a shows the impact on the remaining
<figure><img src="image_5.png" /><figcaption>FIG. 4. Ensemble mean SST minus 140-yr mean control SST (K) during (a) the hiatus period (years 6–10) and (b) the surge period (years 11–15). (c), (d) As in (a) and (b), but for 0–1000-m heat content (TJ). (e), (f) As in (a) and (b), but for net surface heat flux (W m⁻²). Heat flux is positive downward. The thin black contour lines denote values that are significantly different from the control simulation at the 95% confidence level.</figcaption></figure>

shortwave (SW) flows. The reflected SW radiation at the TOA (black) is initially reduced by 0.25 W m⁻² (a 5-yr filter is applied to the data) and then gradually returns to its control value over the next 20 years. The atmospheric absorption (red) is also initially reduced, but returns to the control value earlier, after about 12 years, and overshoots slightly in the surge period. The reduction in both reflection and absorption by the atmosphere means that the incident surface SW radiation (green and cyan) is increased by about 0.6 W m⁻², with the change in surface reflection (blue) being negligible. The incident surface SW radiation returns to its unperturbed value after about 12 years (2 years after the diffusivity returned to normal).

The impact on nonsolar energy flows is depicted in Fig. 5b. Outgoing LW radiation at the TOA (black) displays a symmetrical pattern with a reduction in the first 10 years (consistent with reduced SST) and an increase in years 11–20 consistent with temperatures warmer than the control. The biggest changes occur in the surface upwelling (red) and downwelling (green) LW, which

decrease by about 1 W m⁻² relative to the control during the hiatus but increase rapidly by 2 W m⁻² once the diffusivity returns to its control value. Thus, unlike the shortwave, the LW terms strongly overshoot their original values in the surge regime. Of the turbulent heat flux terms, sensible heat (blue) shows negligible variation, while latent heat (cyan) bears a strong resemblance to atmospheric absorption (weaker latent heat loss and reduced SW absorption during the hiatus period may be linked via reduced atmospheric humidity). Both the latter terms remain slightly higher than their control values in the last 10 years of the perturbed simulations.

Global mean ocean heat uptake, net TOA heat flux, and global mean heat content tendency show very similar variation as required by energy conservation, since most of the heat capacity resides in the ocean; hence, only the ocean heat uptake is illustrated (Fig. 6). The control displays large decadal fluctuations superimposed on a global net heat loss (TOA imbalance) of about -0.2 W m⁻² (vertical black line), indicative of long-term drift. This is
<figure><img src="image_6.png" /></figure>

FIG. 5. Ensemble annual mean anomalous global energy flows (W m⁻²): (a) Shortwave terms (black indicates TOA upward shortwave, red is atmospheric shortwave absorption, green is incident surface shortwave, blue is reflected shortwave, and cyan is oceanic shortwave absorption) and (b) longwave and turbulent heat flux terms. In (b), black indicates outgoing longwave radiation at the TOA, red is surface upwelling longwave radiation, green is surface downwelling longwave radiation, blue is sensible heat flux, and cyan is latent heat flux. Colored horizontal lines denote zero anomaly. To allow the individual anomalies to be seen clearly, each quantity has had an offset subtracted: black is no offset, red is -1 W m⁻², green is -2.5 W m⁻², blue is -3 W m⁻², and cyan is -4 W m⁻². Annual means are plotted, and a 4-yr boxcar filter is applied to smooth the data. Thick portions of the curves indicate values that are significantly different from the control simulation at the 95% confidence level. The dotted vertical line separates the hiatus period (years 1–10; doubled vertical diffusivity) from the surge period (years 11–20; standard vertical diffusivity).

associated with the model Antarctic Bottom Water (see Fig. 10a) in the deep Southern Ocean and is unlikely to affect the signals in the top 1000 m seen in Figs. 2–4.

The diffusivity perturbation causes substantially increased ocean heat uptake of about +0.6 W m⁻² relative to the control. This elevated heat uptake declines steeply between years 8 and 12 and becomes indistinguishable from the control. The ocean heat uptake thus displays an asymmetry: heat gained during the hiatus is not fully discharged from the ocean in the subsequent surge (within the time scale of our experiments).

### c. Ocean and atmosphere surface temperature response

The ensemble mean global average SST is plotted in Fig. 7 (red) together with the global average SST from the control (vertical black line) and one ensemble member (green), which was integrated for an extended period after the diffusivity was restored to its original value. The mean SST over the 140-yr portion of the

control was 17.95°C. In response to the reduced diffusivity, the ensemble mean global average SST reduces by about 0.15 K to 17.80°C. When the diffusivity is returned to its original value at the beginning of year 11, global average SST increases steeply by ~0.25 K in 2–3 years, although in strict terms the hiatus ended and the surge began when global temperature stopped declining and began to rise around year 8 [this is seen in the unfiltered data (green) and is not an artifact of time filtering]. Subsequently the ensemble mean SST stabilizes at about 18.05°C. The green curve indicates the evolution of a single ensemble member rather than the ensemble mean, which was integrated for a much longer period than the others. The surge in temperatures in years 8–12 leaves the global mean SST higher than the control, possibly for several decades or more. During the surge and hiatus periods the ensemble mean SST is significantly (95% confidence) below/above the control, but the ensemble spread remains similar to the temporal variability of the control. In contrast in the transition
<figure><img src="image_7.png" /></figure>

FIG. 6. Ensemble mean global annual mean net ocean heat uptake (net ocean surface downward heat flux) (W m⁻²; thick red line) and ±1 ensemble standard deviation (thin red lines). Annual means are plotted, and a 4-yr boxcar filter is applied to smooth the data. The vertical black line at the left is centered on the 140-yr mean value from the control simulation and indicates ±1 standard deviations (or ±2 standard errors). The dotted vertical line separates the hiatus period (years 1–10; doubled vertical diffusivity) from the surge period (years 11–15; standard vertical diffusivity).

period (years 10–12), the SST anomaly is not significantly different from the control, but the ensemble spread is smaller than the control variability. This suggests a potentially high predictability of this transition phase.

Figure 8 provides a scatterplot of annual mean SAT versus SST for the control (red circles) and the individual ensemble members for the first 10 years of each simulation (black) and the last 10 years (blue). SAT correlates very well with SST so a hiatus in SST is equivalent to one in SAT on an annual mean and global average basis. The same linear relationship between SST and SAT seen in the control holds for both surge and hiatus periods. Spatial and seasonal differences in this relationship are of interest but are not pursued further in this paper.

<figure><img src="image_8.png" /></figure>

FIG. 7. Ensemble mean global mean SST (°C; thick red line) and ±1 ensemble standard deviation (thin red lines). The green line shows global mean SST for an individual ensemble member run out for 50 years. Monthly means are plotted, and a 36-month boxcar filter is applied to smooth the data. The vertical black line at the left is centered on the 140-yr mean value from the control simulation and indicates ±1 standard deviations (or ±2 standard errors). The dotted vertical line separates the hiatus period (years 1–10; doubled vertical diffusivity) from the surge period (years 11–20; standard vertical diffusivity).

### d. Analysis in temperature space

A consistent finding is that following the relaxation of conditions that force a hiatus, the surface warms rapidly and remains warmer than the previous long-term average state. We now investigate this behavior using water-mass transformation analysis.

#### 1) WATER-MASS TRANSFORMATION

We first apply the water-mass transformation framework, Eqs. (12) and (13), to the control simulation. The net trend $\partial V/\partial t$ (Fig. 9a) is very small in most temperature classes and latitudes, but $\partial^2 V/\partial t \partial \theta$ is negative between 0° and 2°C (80°S–20°N), representing a loss of volume. This is counterbalanced by positive $\partial^2 V/\partial t \partial \theta$ at the same latitude range, between -2° and 0°C, representing a gain in volume. Thus, there is cooling in the deep Southern Ocean (~0°C) consistent with the overall negative TOA balance.

From Eq. (12) surface fluxes, ocean circulation and mixing compete to create a net trend in the volume of seawater in temperature classes. In Figs. 9b–d, each panel shows one term from Eq. (12). Values plotted for each latitude $\phi$ and temperature $\theta$ represent all the fluid with temperature $\le \theta$ and north of $\phi$. Therefore, as with streamfunction plots in depth–latitude space, a convergence (divergence) of isolines represents accumulation (loss) of volume at a particular $\theta$ and $\phi$. Along isolines, input and output cancel: water transformed from the class (and/or latitude) below is counterbalanced by the same amount of water being transferred to the class (and/or latitude).
<figure><img src="image_9.png" /><figcaption>FIG. 8. Scatterplot of global annual mean SAT vs SST (°C) for control (red circles), hiatus period (years 1–10 of the ensemble experiment; black symbols) and surge period (years 11–20 of the ensemble experiment; blue symbols). Individual ensemble members are differentiated by different symbols: circles, stars, squares, and triangles.</figcaption></figure>

latitude) above. Hence, we can draw arrows connecting sources of water with particular properties to sinks, or closed loops representing continuous transformation pathways. The reader should remember that the 2D temperature–latitude plots in Fig. 9 represent complex water pathways in 3D physical space (Fig. 1).

Figure 9b shows the overturning streamfunction $\psi$ in temperature space [Eq. (1)]. It is calculated by integrating up lateral transports of fluid of various temperatures, so in principle the diagnosed apparent “vertical” diathermal motion may be associated with volume trends $\partial V/\partial t$ rather than actual warming or cooling. However, since the trends are weak (Fig. 9a), the diagnosed diathermal motion should be accurate. Negative and positive values respectively indicate clockwise and anticlockwise transformation cells. Shallow subtropical cells (STCs), 25–30 Sv (1 Sv $\equiv 10^6 m^3 s^{-1}$) in strength, are apparent, centered at about 22°C, 10°S in the Southern Hemisphere and, slightly cooler, 20°C, 10°N in the north. These deliver cold water upwelled near the equator (the sharp meridional gradient between 15° and 25°C at 0° latitude implies transfer from cold to warm temperatures as explained in section 2b) to the downwelling/subduction regions in the subtropics (Fig. 1 gives a 3D view of this process). The model North Atlantic Deep Water (NADW) cell (i.e., the AMOC) carries ~25 Sv of surface waters northward into the North Atlantic and returns cold (3°–7°C), deep waters southward. These upwell at high middle southern latitudes (~40°S) and rejoin the surface circulation. The Antarctic Bottom Water cell downwells ~25 Sv of cold surface waters in the Southern Ocean and carries them northward as near-bottom currents (at ~0°C). These currents rise and warm slightly at middle high latitudes in the Northern

<figure><img src="image_10.png" /><figcaption>FIG. 9. Transformation streamfunctions (Sv) for the 140-yr control simulation in temperature space: (a) time derivative, (b) ocean circulation, (c) surface heat flux, and (d) mixing. Blue shading with arrows indicates clockwise transformation, and red shading with arrows indicates counterclockwise transformation. Black lines show the average SST at each latitude and $\pm 1$ standard deviation.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>FIG. 10. Latitudinal divergence of the volume flux (Sv per degree of latitude) across isopycnals (transformation rate) (a) at the equator (3°S–3°N) and (b) for the “off equatorial” regions (average of 18°–12°S and 12°–18°N) as a function of temperature associated with surface heat flux (green), overturning (red), and mixing (blue). Green circles demarcate the boundaries between the prevailing thermal regimes.</figcaption></figure>

Hemisphere and join the NADW, making their way back southward.

Air-sea fluxes directly transform water shifting it from one temperature class to another. The associated streamfunction $F$ [Eq. (10)] is analogous to $\psi$ (Fig. 9c) but differs in that it is defined from the (known) diathermal flux; it gives diagnosed (unrealistic) horizontal transports consistent with weak volume trends and no mixing. Air-sea fluxes only act on a temperature class if the bounding surface of that class reaches to the surface at the latitude in question (or close enough that penetrative solar radiation passes through its boundary) as illustrated (section 2b) by the purple region of Fig. 1d. Thus, there is a blank region in Fig. 9c between 40°S and 30°N below 10°C where air-sea fluxes cause no temperature change. Air-sea fluxes can make warm waters even warmer (tropical heating) and cool waters even cooler (subtropical/high-latitude cooling). Water in the 23°C range is transformed by air-sea fluxes into warmer waters (up to 30°C) mainly in the tropics. Waters in this class are also cooled by air-sea fluxes between the subtropics and midlatitudes to temperatures around 15°C. Farther poleward, waters around 5°C are transformed into cold deep water below 0°C, largely in the Southern Ocean and the North Atlantic.

Mixing, obtained as a residual $(\partial D_{\text{diff}}/\partial \theta)$ in Eq. (12), opposes air-sea fluxes, removing extremes of temperature and adding volume to intermediate temperature classes (Fig. 9d). However, mixing also acts within the ocean interior, where air-sea fluxes cannot penetrate; for example, mixing is active between the equator and 40°S, transforming water from the 15° to the 20°C range, and balances the NADW overturning cell at around 50°N where waters between 5° and 10°C are transformed into the range 0°–5°C.

## 2) TRANSFORMATION AND FORMATION RATES

We now focus on the tropics and subtropics where a large part of the surge signal originates and examine the horizontal derivative of the transformation rate. Following Eq. (12), in Fig. 10a we display the latitudinal divergence of the volume flux across isotherms (the transformation rate) at the equator (averaged over 3°S–3°N) as a function of temperature due to air-sea fluxes, mixing, overturning, and volume change. This is essentially a section through Figs. 9a–d at the equator, differentiated with respect to latitude. We plot curves for the control simulation to elucidate the fundamental balances in the model and subsequently investigate how these balances change during the hiatus and surge in the perturbed ensemble.
<figure><img src="image_12.png" /><figcaption>FIG. 11. Formation rate (Sv per degree of latitude) (a) at the equator (3°S–3°N) and (b) for the off-equatorial regions (average of 18°–12°S and 12°–18°N) in the control simulation associated with overturning (red), surface heat flux (green), and mixing (blue) for the four thermal regimes (i–iv) defined in Fig. 10.</figcaption></figure>

At the equator, air–sea fluxes (green) cause transformation in temperature classes 17°–32°C (there is no direct effect on isotherms below 17°C because these do not outcrop at the equator). The lines have a negative gradient (formation rate) between 26° and 31°C; hence, air–sea fluxes tend to generate water masses in this temperature range [F has a negative sign in Eq. (12)]. Between 17° and 25°C, the lines have a positive gradient and air–sea fluxes remove water in these classes. The overturning (red) acts to remove water between 23° and 30°C and add water between 10° and 23°C. Mixing (blue) differs from the other terms as the associated formation rate has three regimes: between 28° and 32°C mixing removes water, between 21° and 28°C it adds water, and between 10° and 21°C it removes water. This is consistent with mixing tending to eliminate extremes. The result is that at the equator there are four thermal regimes. In each, two out of the three terms (in formation rate) are in approximate balance and the third is relatively unimportant. Thus between 28° and 32°C the predominant balance is mixing versus air–sea fluxes (equatorial regime iv), between 25° and 28°C it is mixing versus overturning (regime iii), between 21° and 25°C it is mixing versus air–sea fluxes (regime ii), and between 10° and 21°C it is mixing versus overturning (regime i). Below 15°C there is little coherent signal. In all four regimes, mixing is one of the dominant terms. Volume change is a very small term over most temperature classes, averaged over the 140-yr control simulation.

Figure 10b shows transformation rates for the off-equatorial regions (average over 18°–12°S and 12°–18°N; the strongest subsurface anomalies in Fig. 4 are between 20°S and 20°N). Transformation rates are much smaller relative to those at the equator: 1–2 Sv per degree of latitude as compared with 8–10 Sv per degree of latitude, and the roles of the individual processes differ. Air–sea fluxes add volume between 28° and 32°C and between 20° and 25°C, and they remove volume between 25° and 28°C. The overturning behaves in an opposite manner than at the equator, adding volume between 25° and

28°C and removing it between 20° and 25°C. Mixing removes volume between 28° and 32°C (off-equatorial regime iv) and adds it between 20° and 25°C (regime iii). Between 10° and 20°C, overturning and mixing balance, mixing adding volume and overturning removing it (regime ii). The situation is reversed between 5° and 10°C with overturning adding volume and mixing removing it (regime i). Details vary between the northern and southern off-equatorial regions (not shown), but the regimes are similar. Again, net volume change is a very small term. Although the thermal regimes are different from those at the equator, the boundaries between the regimes occur at about the same temperatures.

In both equatorial and off-equatorial regions, the corresponding formation rate curves for the perturbed ensemble (not shown) are similar in shape to the control and the boundaries between the thermal regimes occur at the same temperatures.

We further characterize the thermal regimes by calculating the formation rate terms (section 2b and Fig. 1) for the control simulation. Figure 11a shows the net volume flux convergence for each regime i–iv, associated with overturning (red), air–sea fluxes (green), and mixing (blue) in the equatorial region. The net volume inflation in each of these regimes is very small, and there is a near-exact balance between the three terms. In regime iv, air–sea fluxes create about 6 Sv per degree of latitude of water (36 Sv in total for the 3°S–3°N region) while mixing and overturning remove 4 and 2 Sv, respectively, per degree of latitude. Mixing delivers 3 Sv per degree of latitude of water into regime iii, largely balanced by overturning. Regime ii is effectively a reversed version of regime iv (loss of ~6 Sv per degree of latitude by air–sea fluxes, and gains of 4 and 2 Sv per degree of latitude by mixing and overturning, respectively). Similarly, regime i is a reversed version of regime iii. This suggests that regimes i and iii and ii and iv are causally connected. The situation is subtly different in the off-equatorial regions (Fig. 11b) where regimes iii and iv
both have mixing balancing overturning and air-sea fluxes, but with opposite signs. Similarly, regimes i and ii both have mixing largely balancing overturning, but with opposite signs.

In both equatorial and off-equatorial regions, there are regimes in which air-sea fluxes add volume and others in which they remove it, suggesting that the regimes are not all stacked vertically, but more likely horizontally (e.g., warming in the western subtropics and cooling in the eastern). This would explain the causal connection between off-equatorial regimes iii and iv, for example, if the spatial regions corresponding to these regimes are connected by atmospheric or oceanic teleconnections. Figure 1d, depicting a strong east-west temperature gradient in the near-surface of the northern tropical Pacific, supports this conjecture. Albeit only for a single month from one ensemble member, it suggests that above the 20°C isotherm off-equatorial regimes ii-iv are stacked east-west and are linked by the northern STC.

### 3) CHANGES OF VOLUME IN TEMPERATURE CLASSES

Figure 12 shows the logarithm of the ratio (ensemble mean/control mean) of the globally integrated volume in temperature space, calculated in 1-K intervals and based on monthly means, with a 12-month boxcar filter applied to smooth the data (before taking the logarithm). Negative and positive values respectively indicate reduced and increased volume in temperature classes in the ensemble mean versus the control. In years 1–10 of the perturbed simulations, there is a depletion of water in the range 27°–30°C. Conversely there is a small increase in volume at cooler temperatures, confined mainly to 25°–27°C for the first five years but subsequently extending to cooler temperatures (down to 15°C). From year 11 onward there is an abrupt increase in the volume of the warmest waters, 27°–30°C, and a reduction of volume in the 25°–27°C range. The volume of water in the 18°–25°C range decreases more gradually, and in the 15°–16°C range there is no reduction at all; indeed there is a long-term increase in the volume of these waters with little sign of impending decrease a decade after the surge began. There are also notable increases of volume at around 10°C, still present at year 20. The increase in volume in the range from –1° to 0°C is a reflection of the drift in the control associated with the long-term negative TOA imbalance of –0.2 W m⁻².

Extending the plots to cover latitude, the logarithm of the absolute volume per unit latitude for temperature classes in the control simulation (Fig. 13a) shows the largest volume of seawater in the coldest temperature classes, and a comparatively low volume in the warmest

classes (black lines indicate the zonal and time average SST ± 1 standard deviation at each latitude).

Figure 13b shows the logarithm of the ratio (ensemble mean/control mean) of globally integrated volume for the first (blue) and second (red) decades of the perturbed simulations. Negative and positive values respectively indicate a loss and gain of volume relative to the control. In the first decade there is a strong decrease in volume above 25°C and a compensating but more modest increase in the 15°–25°C layer. In years 11–20 the warmer layers regain some of the volume that was lost in the hiatus decade at the expense of loss in volume of waters between 20° and 25°C; however, the volume gained between 15° and 20°C is retained in the subsequent surge and new equilibrium period. The 15° and 10°C volume classes, in particular, remain much thicker than in the control in years 11–20.

Next we examine the change in volume per unit latitude in temperature classes as a function of latitude for the hiatus (years 1–10) and the surge (years 11–20). In years 1–10 (Fig. 13c) there is considerable loss of volume of the warmest tropical waters (>25°C) between 20°S and 20°N. Below these waters we note an increased volume of waters between 15° and 25°C, extending poleward to 30°S/30°N, slightly greater in latitudinal extent than the surface cooling/loss of volume. There is also increased volume of the warmest waters in the Southern Ocean (80°–40°S) and the northern midlatitudes (20°–50°N).

In years 11–20 (Fig. 13d), we see a widespread increase in volume of the warmest waters, including a large increase in the tropics where the volume was previously depleted. The Southern Ocean remains as warm during the surge period as it was during the hiatus, but the northern midlatitudes see a very strong increase

<figure><img src="image_13.png" /></figure>

FIG. 12. Logarithm to base 10 of the ratio of ensemble mean to control seawater volume in temperature classes (1-K bins centered on integral values), evaluated over the global ocean. The dotted vertical line separates the hiatus period (years 1–10; doubled vertical diffusivity) from the surge period (years 11–15; standard vertical diffusivity).
<figure><img src="image_14.png" /><figcaption>FIG. 13. (a) Logarithm to base 10 of the control simulation seawater volume in temperature classes (1-K bins centered on integral values) evaluated at each latitude. (b) Logarithm to base 10 of the ratio of ensemble mean to control global seawater volume in temperature classes (1-K bins centered on integral values). Blue indicates the average over the hiatus period; red is the average over the surge period. The vertical black line indicates a zero value. Also shown is logarithm to base 10 of the ratio of ensemble mean to control seawater volume in temperature classes (1-K bins centered on integral values), evaluated at each latitude and averaged over (c) the hiatus period, years 1–10, and (d) the surge period, years 11–20. Black lines show the average SST at each latitude and ±1 standard deviation.</figcaption></figure>

in the volume of warmer waters. At temperatures between 20° and 25°C, the extra volume deposited during the hiatus is lost, but below 20°C the increased volume of water gained at the expense of the warmest waters remains throughout the subsequent surge and the new equilibrium period.

In Figs. 13c and 13d we have superimposed the circulation cells from Fig. 9b. These demonstrate that the warming of the waters warmer than about 15°C during the hiatus in years 1–10 coincides with the positions of the subtropical cells. The warming and cooling pattern north of 40°N seems associated with the AMOC.

We also investigate pentadal changes to the formation rate in the thermal regimes in the perturbed ensemble (Fig. 14) to explain the changes in volume seen in Figs. 12 and 13. At the equator, across all regimes, during the hiatus the changes to the mixing (blue) create anomalies in the overturning and surface heat fluxes (Figs. 14a,b). The overturning anomalies die down in the early part of the surge (Fig. 14c) and completely disappear in the second pentad (Fig. 14d), leaving a modified balance in the long term where the mixing remains changed, but is balanced by a modified surface heat flux. There is a suggestion that regimes ii and iv and regimes i and iii are opposite to each other in the final pentad.

A similar pattern occurs in the off-equatorial regions where the changes in mixing during the hiatus induce changes in the ocean circulation (Figs. 14e,f), which however die down and disappear during the surge (Figs. 14g,h). The final state once again pitches long-term changes in mixing against changes in air-sea fluxes and regimes iii and iv mirror each other, reinforcing the earlier suggestion that there is a causal connection between them (Fig. 11).

We defer detailed analysis of the mechanism governing this transient behavior to a future study; however, we link the behavior in temperature space in Fig. 14 to depth space and previous studies (e.g., England et al. 2014) by plotting the pentadal evolution of the STCs (Figs. 15a–d). The thin contours show the meridional overturning streamfunction from the control simulation in the tropical (30°S–30°N) Indo-Pacific. There is an ~30-Sv clockwise cell in the northern tropics/subtropics and a somewhat stronger (~40 Sv) anticlockwise cell in the Southern Hemisphere. Both have a meridional extent of ~25° of latitude, with the southern cell extending slightly deeper (~300 m as compared with ~200 m for the northern cell).

In the first pentad (hiatus; Fig. 15a) we see a modest strengthening of the southern cell, but in the second there is a much bigger strengthening (on the order of 5 Sv) of the
<figure><img src="image_15.png" /><figcaption>FIG. 14. Ensemble mean anomalous formation rate (Sv per degree of latitude) (left) at the equator (3°S–3°N) and (right) for the off-equatorial regions (average of 18°–12°S and 12°–18°N) associated with overturning (red), surface heat flux (green), mixing (blue), and volume change (cyan) for the four thermodynamic regimes (i–iv) defined in Fig. 10, for (a),(e) years 1–5, (b),(f) years 6–10, (c),(g) years 11–15, and (d),(h) years 16–20.</figcaption></figure>

northern cell (Fig. 15b). In the third pentad (surge; Fig. 15c) the northern cell returns to normal strength, while the southern cell weakens by ~5 Sv. In the fourth (Fig. 15d), and consistent with Fig. 14, we find both STCs have returned to normal. Note that when only one of the overturning cells increases, we are likely to see opposing changes in the equatorial versus the off-equatorial regions (e.g., in regimes iii and iv the red bars are of opposite sign between the equatorial and off-equatorial regions), since increasing

divergence at the equator is associated with increasing convergence off-equator and vice versa.

Consistent with previous studies (Meehl et al. 2013), increased heat drawdown in the ocean is associated with strengthening STCs, and the reverse is true for surge conditions. Thus, while the perturbed background diffusivity artificially induces a hiatus in our ensemble experiment, heat uptake in the 300-m layer increases due to a realistic physical mechanism. The stronger and

<figure><img src="image_16.png" /><figcaption>FIG. 15. Ensemble mean anomalous (top) overturning streamfunction (Sv) and (bottom) zonal wind (ms⁻¹) in the tropics (color shading) for (a),(e) years 1–5, (b),(f) years 6–10, (c),(g) years 11–15, and (d),(h) years 16–20. The thick black contours and dotted regions denote anomalies that are significantly different from the control simulation at the 95% confidence level. The thin black contours in (a)–(d) denote the control mean streamfunction (Sv).</figcaption></figure>
weaker STCs in Figs. 15b and 15c are in addition associated with strengthened and weakened easterly trade winds, respectively (Figs. 15f,g). In contrast there is no marked signal in the wind stress when the STCs are near normal (Figs. 15e,h). Therefore, ocean heat uptake and climate variability feed back on each other, and artificial changes to the background diffusivity induce a physically realistic circulation response in line with previous studies.

### e. Conceptual model of transition

The previous section demonstrated that, in the model, mixing is always of first-order importance in the thermal balance in the tropics and subtropics (Figs. 10 and 11). Combining this result with the fact that the temperature tendency due to mixing is proportional to the second derivative of the temperature [Eq. (14)], we hypothesize that any modification to one of the forcing terms will modify the temperature leading to changes in stratification, and bringing the mixing term back into balance with its antagonist.

Consider an idealized upper ocean consisting of a 200-m layer in which surface heating $I$ is balanced by diffusion $[\partial(K_v\partial\theta/\partial z)/\partial z]$ similar to regime iv (Fig. 10a), with the further simplification that ocean heat transport divergence $\nabla \cdot (\mathbf{u}T)$ is identically zero. We apply a uniform heating of $20 \text{ W m}^{-2}$ distributed evenly throughout the upper 80 m representing solar input. Initially, we specify zero heat flux through the surface boundary, relying on the upper ocean heating to drive the model to equilibrium. We use the diffusive heat flux at the lower boundary to drive a deeper layer of thickness 100 m balanced by a fixed ocean heat transport divergence of $20 \text{ W m}^{-2}$. This deeper layer is thus similar to regime iii in Fig. 10a.

The upper 200 m is discretized with a forward timestepping scheme and standard second-order accurate finite differences in the vertical (time step $\Delta t = 864$ s; vertical grid spacing $\Delta z = 1$ m). The lower 100 m is treated numerically as a single layer. We assume a constant background diffusion coefficient $K_v = 5 \times 10^{-5} \text{ m}^2 \text{s}^{-1}$, set the initial temperature to $20^\circ\text{C}$ in all the layers, then integrate the model forward for 120 years. At the beginning of year 121 the background diffusivity is increased by 20% to $6 \times 10^{-5} \text{ m}^2 \text{s}^{-1}$; in addition, heat flux feedback to the atmosphere is switched on [Eq. (15)], damping the SST back to the long-term mean (based on years 40–80) with a feedback parameter of $1 \text{ W m}^{-2} \text{ K}^{-1}$. At the beginning of year 131, the background diffusivity is returned to its standard value but the surface heat flux feedback is retained, and the model is integrated for a further 30 years.

The system evolves to equilibrium in about 20 years such that diffusion transports heat downward from the

surface (where it is input by heating) to the deepest layer where it is removed by heat transport divergence and continues without significant change until year 120. When the diffusivity is increased at the beginning of year 121, there is an immediate elevation of the isotherms above 140 m and a depression of isotherms below, resulting in a cold anomaly above 140-m depth and a warm anomaly below 140-m depth (Fig. 16a)—that is, a hiatus.

When the diffusivity is returned to normal at the beginning of year 131, the system experiences a temporary perturbation such that the export of heat by mixing from the upper ocean through the lower boundary is reduced. There is thus a temporary imbalance in favor of the surface heat flux, which modifies the temperature gradient to restore a balanced regime (note that, although we perturbed the system by changing the diffusion coefficient, we could have modified one of the other processes, advection or air-sea fluxes, to create the hiatus).

This behavior is very reminiscent of Fig. 2, albeit with larger anomalies because it is intended to represent behavior near the equator rather than an average over the globe. In Fig. 16b temperature anomaly profiles with depth are plotted before and after the hiatus-surge transition. The hiatus state (red) shows a ~0.7-K reduction in temperature at the surface and a corresponding increase at 200 m. The cyan line shows how the temperature profile changes at year 145 when the diffusivity has been reduced to its original value. A surface temperature increase of 0.2 K over control values, with little depth variation, is seen. This final equilibrium state thus has a similar temperature gradient to the prehiatus state in the top 200 m but is systematically warmer at all depths.

We now compare the zonally averaged response in the full climate model with this conceptual model. Figure 16c shows the ensemble mean zonal average (confined to the Pacific basin) of the temperature in the top 200 m of the equatorial Pacific (3°S–3°N) in the full climate model. We have noted previously that the major temperature anomaly during the simulated hiatus sits below 200 m and remains there during the surge. During the hiatus (red line, representing an average of years 6–10 of the ensemble experiment), temperature above 120 m reduces by up to 0.5 K while below this there is a temperature increase of a similar magnitude. In the first year of the surge (cyan line; average over year 11), temperature at 200 m stays at the same elevated level attained during the hiatus. Above this depth there is an increase in temperature of 0.4 K, approximately independent of depth. Hence in the surge simulated by the climate model, as in the conceptual model, the temperature gradient returns to prehiatus values, but the temperatures themselves are elevated as a result of the increased temperature of the deep layer.
<figure><img src="image_17.png" /><figcaption>FIG. 16. (a) Temperature anomaly (K) predicted by the conceptual model. The dotted vertical lines mark the hiatus period (years 121–130; higher vertical diffusivity) before the surge period (years 131–160; standard vertical diffusivity). (b) Anomaly-depth profile predicted by the conceptual model at hiatus year 125 (red) and surge year 145 (cyan). (c) Ensemble mean zonally averaged equatorial (3°S–3°N) temperature anomaly profile from the coupled climate model experiment as a function of depth in the Pacific Ocean averaged over the hiatus period (red) and the first year of the surge (cyan). (d) As in (c), but for the off-equatorial regions (average of 18°–12°S and 12°–18°N). The cyan line shows the average over year 3 of the surge.</figcaption></figure>

The northern off-equatorial region, 12°–18°N (Fig. 16d), displays similar behavior to that at the equator. Here, the cyan lines are plotted for year 3 of the surge (and not year 1) because it takes longer for the system to attain equilibrium away from the equator (varying Rossby wave speed is the likely explanation). The southern off-equatorial region displays the same behavior except that the warming during the hiatus occurs deeper down at 400 m (not shown).

Interestingly, in the conceptual model (Fig. 16b) we only had to increase the diffusion coefficient by 20% to achieve a realistic response, whereas in HadCM3 we applied a 100% increase. Beyond the simplification inherent to the conceptual model, it is possible that this mismatch comes from other sources of mixing in HadCM3, such as numerical mixing and shear-driven instability (e.g., Megann 2018).

4. Conclusions and discussion

We performed sensitivity experiments using the HadCM3 climate model to understand decadal changes in the rate of

global warming. Specifically, we increased the ocean heat uptake by increasing the vertical diffusivity in the model for 10 years, before returning vertical diffusivity to normal values. The main conclusions are as follows:

* Experiments with increased ocean diffusivity show Pacific-dominated hiatus-like behavior: the Pacific undergoes a transition to a negative interdecadal Pacific oscillation like state. Outgoing longwave radiation at the top of the atmosphere, oceanic latent heat loss, and surface upwelling and downwelling longwave radiation are all strongly reduced due to reduced surface temperatures; however, the ocean absorbs more shortwave radiation due to reduced atmospheric absorption.

* When ocean vertical mixing is increased global surface temperature decreases (i.e., the climate is forced into a hiatus-like state). When ocean vertical mixing returns to normal values the surface temperature increases, and overshoots the climatological values (i.e., a surge occurs).

* The climatological thermal balance in the warmest (near surface) waters of the tropical Pacific Ocean is maintained by mixing balancing surface heat flux.
Hiatus to surge mechanism

<figure><img src="image_18.png" /><figcaption>FIG. 17. Schematic of the hiatus-to-surge mechanism described in this paper: (a) conditions during hiatus, (b) conditions during surge, and (c) conditions after surge.</figcaption></figure>

In contrast, in the cooler layer below, mixing is balanced by ocean heat transport divergence. Related but slightly different balances prevail in the off-equatorial regions.

* The overshoot into surge conditions is due to long-term depression of deep (~15°C) isotherms following the hiatus (put simply, the ocean has no mechanism to quickly remove the heat deposited at depth). To maintain thermal balance, the upper-layer temperature, the most responsive variable, must rise to higher than normal values, via increased surface heat input, in order to restore the thermal gradient, which supports mixing (see schematic in Fig. 17).

* Other things being equal, the changes in stratification caused by a hiatus increases the probability that a subsequent surge will occur.

Because of our use of fundamental thermal balances in this study, we are confident that our conclusions are correct in their essentials. However, there are some caveats to consider:

* We have used one particular model of relatively coarse resolution (although it maintains a very stable climate similar to the present day).

* The thermal balances may differ quantitatively from model to model.

* Our model does not resolve ocean mesoscale eddies, which play an important role in the thermal balance of the ocean (Griffies et al. 2015). In mitigation our model employs the Gent-McWilliams eddy parameterization scheme.

Our results are consistent with previous observational and modeling studies, in particular Maher et al. (2018). The authors modified Pacific trade winds to induce a hiatus in a similar way to our modification of the mixing, and show similar results in terms of the changes in subtropical cells, heat content, the transfer of excess heat to the Indian Ocean, and the legacy of anomalously high heat content even under a reversal of the perturbation. This agreement, as well as our Figs. 14 and 15, suggests a tight coupling between the three competing processes of surface fluxes, wind-driven ocean circulation, and mixing, identified in the present study. In a similar way, observations (Nieves et al. 2015) confirm the transfer of heat during the recent hiatus period to the Indian Ocean.

Last, we consider the wider implications of our findings. As is widely appreciated, when the deep ocean is heated on decadal time scales, the additional heat cannot easily and quickly be returned to the atmosphere. We have shown that this additional heat has a legacy
impact on the global mean surface temperature on decadal and longer time scales. The additional heat uptake during hiatus periods (which is well attested) conditions the Earth system to warm even more strongly when the hiatus ends. While we have engineered a hiatus by modifying the vertical mixing, in principle we would expect to see the same behavior if the heat enters the deep ocean by other means (e.g., via modification of the shallow overturning cells from wind stress perturbations, or changes to the surface heating and subsequent subduction, mechanisms that have been hypothesized to have driven the hiatus in the early 2000s). Idealized/conceptual one-dimensional models of Earth's energy balance do not currently include the effects described in this paper.

Our results also have important implications for decadal and longer predictions using model-based forecast systems. In particular, we suggest that accurate initialization of the subsurface ocean may be a key element of such systems as inaccuracies could result in erroneous predictions of hiatus or surge conditions [consistent with the findings of Sévellec and Fedorov (2017) and Germe et al. (2018)]. Similarly, inaccurate simulation of forced hiatus periods could lead to systematic biases or additional uncertainty in centennial climate projections.

Further work should establish the robustness of our result in a wider variety of models (including at higher ocean resolution). The water-mass transformation method we used could also in principle be applied to observational data, although in practice great care would be needed to take into account observational uncertainties and internal variability.

*Acknowledgments.* Authors Sinha, Robson, and Sévellec were supported by the U.K. Natural Environment Research Council (NERC) under the Securing Multidisciplinary Understanding and Prediction of Hiatus and Surge events (SMURPHS) project, Grants NE/N005686/1 (Sinha), NE/N006054/1 (Robson), and NE/N005767/1 (Sévellec). Author Nurser was supported by the NERC Transient Tracer-Based Investigation of Circulation and Thermal Ocean Change (TICTOC) project (Grant NE/P019293/1). We thank Jonathan Gregory for useful discussions. All model output used in this paper will be made available upon request to the corresponding author.

## REFERENCES

Alexander, L. V., and Coauthors, 2013: Summary for policymakers. *Climate Change 2013: The Physical Science Basis*, T. F. Stocker et al., Eds., Cambridge University Press, 3–29, https://www.ipcc.ch/site/assets/uploads/2018/03/WG1AR5_SummaryVolume_FINAL.pdf.

Balmaseda, M. A., K. E. Trenberth, and E. Källen, 2013: Distinctive climate signals in reanalysis of global ocean heat content. *Geophys. Res. Lett.*, 40, 1754–1759, https://doi.org/10.1002/grl.50382.

Chen, X., and K.-K. Tung, 2014: Varying planetary heat sink led to global-warming slowdown and acceleration. *Science*, **345**, 897–903, https://doi.org/10.1126/science.1254937.

Drijfhout, S. S., A. T. Blaker, S. A. Josey, A. J. G. Nurser, B. Sinha, and M. A. Balmaseda, 2014: Surface warming hiatus caused by increased heat uptake across multiple ocean basins. *Geophys. Res. Lett.*, **41**, 7868–7874, https://doi.org/10.1002/2014GL061456.

England, M. H., and Coauthors, 2014: Recent intensification of wind-driven circulation in the Pacific and the ongoing warming hiatus. *Nat. Climate Change*, **4**, 222–227, https://doi.org/10.1038/nclimate2106.

Farneti, R., A. Dwivedi, F. Kucharski, F. Molteni, and S. M. Griffies, 2014: On Pacific subtropical cell variability over the second half of the twentieth century. *J. Climate*, **27**, 7102–7112, https://doi.org/10.1175/JCLI-D-13-0070.1.

Gastineau, G., A. Friedman, M. Khodri, and J. Vialard, 2019: Global ocean heat content redistribution during the 1998–2012 Interdecadal Pacific Oscillation negative phase. *Climate Dyn.*, **53**, 1187–1208, https://doi.org/10.1007/s00382-018-4387-9.

Germe, A., F. Sévellec, J. Mignot, A. V. Fedorov, S. Nguyen, and D. Swingedouw, 2018: The impacts of oceanic deep temperature perturbations in the North Atlantic on decadal climate variability and predictability. *Climate Dyn.*, **51**, 2341–2357, https://doi.org/10.1007/s00382-017-4016-z.

Gordon, C., C. Cooper, C. A. Senior, H. Banks, J. M. Gregory, T. C. Johns, J. F. B. Mitchell, and R. A. Wood, 2000: The simulation of SST, sea ice extents and ocean heat transports in a version of the Hadley Centre coupled model without flux adjustments. *Climate Dyn.*, **16**, 147–168, https://doi.org/10.1007/s003820050010.

Griffies, S. M., and Coauthors, 2015: Impacts on ocean heat from transient mesoscale eddies in a hierarchy of climate models. *J. Climate*, **28**, 952–977, https://doi.org/10.1175/JCLI-D-14-00353.1.

Guemas, V., F. J. Doblas-Reyes, I. Andreu-Burillo, and M. Asif, 2013: Retrospective prediction of the global warming slowdown in the past decade. *Nat. Climate Change*, **3**, 649–653, https://doi.org/10.1038/nclimate1863.

Hedemann, C., T. Mauritsen, J. Jungclaus, and J. Marotzke, 2017: The subtle origins of surface-warming hiatuses. *Nat. Climate Change*, **7**, 336–339, https://doi.org/10.1038/nclimate3274.

Kosaka, Y., and S.-P. Xie, 2013: Recent global-warming hiatus tied to equatorial Pacific surface cooling. *Nature*, **501**, 403–407, https://doi.org/10.1038/nature12534.

Lee, M.-M., A. C. Coward, and A. J. G. Nurser, 2002: Spurious diapycnal mixing of the deep waters in an eddy-permitting global ocean model. *J. Phys. Oceanogr.*, **32**, 1522–1535, https://doi.org/10.1175/1520-0485(2002)032<1522:SDMOTD>2.0:CO2.

Lee, S.-K., W. Park, M. O. Baringer, A. L. Gordon, B. Huber, and Y. Liu, 2015: Pacific origin of the abrupt increase in Indian Ocean heat content during the warming hiatus. *Nat. Geosci.*, **8**, 445–449, https://doi.org/10.1038/ngeo2438.

Liu, W., S.-P. Xie, and J. Lu, 2016: Tracking ocean heat uptake during the surface warming hiatus. *Nat. Commun.*, **7**, 10926, https://doi.org/10.1038/ncomms10926.

Maher, N., M. H. England, A. Sen Gupta, and P. Spence, 2018: Role of Pacific trade winds in driving ocean temperature during the recent slowdown and projections under a wind trend reversal. *Climate Dyn.*, **51**, 321–336, https://doi.org/10.1007/s00382-017-3923-3.

Meehl, G., J. Arblaster, J. Fasullo, H. Aixue, and K. Trenberth, 2011: Model-based evidence of deep-ocean heat uptake during
surface-temperature hiatus periods. *Nat. Climate Change*, **1**, 360–364, https://doi.org/10.1038/nclimate1229.

—, A. Hu, J. M. Arblaster, J. Fasullo, and K. E. Trenberth, 2013: Externally forced and internally generated decadal climate variability associated with the interdecadal Pacific oscillation. *J. Climate*, **26**, 7298–7310, https://doi.org/10.1175/JCLI-D-12-00548.1.

—, H. Teng, and J. M. Arblaster, 2014: Climate model simulations of the observed early-2000s hiatus of global warming. *Nat. Climate Change*, **4**, 898–902, https://doi.org/10.1038/nclimate2357.

Megann, A., 2018: Estimating the numerical diapycnal mixing in an eddy-permitting ocean model. *Ocean Modell.*, **121**, 19–33, https://doi.org/10.1016/j.ocemod.2017.11.001.

Moat, B. I., and Coauthors, 2019: Insights into decadal North Atlantic sea surface temperature and ocean heat content variability from an eddy-permitting coupled climate model. *J. Climate*, **32**, 6137–6161, https://doi.org/10.1175/JCLI-D-18-0709.1.

Nieves, V., J. K. Willis, and W. C. Patzert, 2015: Recent hiatus caused by decadal shift in Indo-Pacific heating. *Science*, **349**, 532–535, https://doi.org/10.1126/science.aaa4521.

Nurser, A. J. G., R. Marsh, and R. G. Williams, 1999: Diagnosing water mass formation from air-sea fluxes and surface mixing. *J. Phys. Oceanogr.*, **29**, 1468–1487, https://doi.org/10.1175/1520-0485(1999)029<1468:DWMFFA>2.0.CO;2.

Risbey, J., 2015: Free and forced climate variations. *Nature*, **517**, 562–563, https://doi.org/10.1038/517562a.

Sévellec, F., and A. V. Fedorov, 2017: Predictability and decadal variability of the North Atlantic Ocean state evaluated from a realistic ocean model. *J. Climate*, **30**, 477–498, https://doi.org/10.1175/JCLI-D-16-0323.1.

—, and S. S. Drijfhout, 2018: A novel probabilistic forecast system predicting anomalously warm 2018–2022 reinforcing the long-term global warming trend. *Nat. Commun.*, **9**, 3024, https://doi.org/10.1038/s41467-018-05442-8.

—, B. Sinha, and N. Skliris, 2016: The rogue nature of hiatuses in a global warming climate. *Geophys. Res. Lett.*, **43**, 8169–8177, https://doi.org/10.1002/2016GL068950.

Smith, D., and Coauthors, 2015: Earth's energy imbalance since 1960 in observations and CMIP5 models. *Geophys. Res. Lett.*, **42**, 1205–1213, https://doi.org/10.1002/2014GL062669.

Song, Y., Y. Q. Yu, and P. F. Lin, 2014: The hiatus and accelerated warming decades in CMIP5 simulations. *Adv. Atmos. Sci.*, **31**, 1316–1330, https://doi.org/10.1007/s00376-014-3265-6.

Tanaka, Y., I. Yasuda, and H. Hasumi, 2012: Effects of the 18.6-yr modulation of tidal mixing on the North Pacific bidecadal climate variability in a coupled climate model. *J. Climate*, **25**, 7625–7642, https://doi.org/10.1175/JCLI-D-12-00051.1.

Walin, G., 1982: On the relation between sea-surface heat-flow and thermal circulation in the ocean. *Tellus*, **34A**, 187–195, https://doi.org/10.1111/j.2153-3490.1982.tb01806.x.

Watanabe, M., H. Shiogama, H. Tatebe, M. Hayashi, M. Ishii, and M. Kimoto, 2014: Contribution of natural decadal variability to global warming acceleration and hiatus. *Nat. Climate Change*, **4**, 893–897, https://doi.org/10.1038/nclimate2355.

Wild, M., D. Folini, C. Schär, N. Loeb, E. G. Dutton, and G. König-Langlo, 2013: The global energy balance from a surface perspective. *Climate Dyn.*, **40**, 3107–3134, https://doi.org/10.1007/s00382-012-1569-8.

Williams, R. G., V. Roussenov, M. S. Lozier, and D. Smith, 2015: Mechanisms of heat content and thermocline change in the subtropical and subpolar North Atlantic. *J. Climate*, **28**, 9803–9815, https://doi.org/10.1175/JCLI-D-15-0097.1.

Xie, S.-P., Y. Kosaka, and Y. M. Okumura, 2016: Distinct energy budgets for anthropogenic and natural changes during global warming hiatus. *Nat. Geosci.*, **9**, 29–33, https://doi.org/10.1038/ngeo2581.