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

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 papera globalclimate model is forced to adoptan oceanstate corresponding to a hiatus [with negative interdecadal Paciﬁc oscillation (IPO) and other surface features typical of a hiatus] by artiﬁcially 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 modiﬁed, warmer climate for decades. Water- mass transformation diagnostics indicate that the heat budget of the tropical Paciﬁc Ocean is a balance between large opposite-signed terms: surface heating/cooling resulting from air–sea heat ﬂux is balanced by vertical mixing and oceanheat transport divergence. During the artiﬁcial hiatus, excess heat becomestrapped just above the thermocline and there is a weak vertical thermal gradient (due to the high artiﬁcial 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 3 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.


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-todecadal 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 eruptions) versus internal variability  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).
Denotes content that is immediately available upon publication as open access.
During the recent global surface warming hiatus of the 2000s Smith et al. (2015) estimated a reduction in net topof-atmosphere (TOA) radiation of 20.31 6 0.21 W m 22 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 22 ) 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 andTung 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.

a. Climate model and experiment design
We use the HadCM3 climate model (Gordon et al. 2000), which maintains a stable climate simulation without flux adjustments despite relatively coarse grid resolution (1.258 3 1.258 ocean and 3.758 3 2.58 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 25 m 2 s 21 at the surface to ;15 3 10 25 m 2 s 21 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 5 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(u, f) at which water crosses some isotherm with potential temperature u integrated zonally over all longitudes l and meridionally north of some latitude f. This is related by mass (or volume in a Boussinesq fluid) conservation to the meridional influx across f, of water colder than u (where y is meridional velocity and R is the radius of Earth) and the rate of change of the volume V(u, f) north of f and colder than u by This transformation G(u, f) represents water crossing the isotherm from cool to warm, and thus requires heat input; to link G(u, f) to heat inputs we consider (Fig. 1b)  Wild et al. (2013), the TOA outgoing solar comprises all reflected solar radiation including that reflected from the surface.
HadCM3 control simulation Wild et al. (2013)  the (Boussinesq) equation for the heat content of water colder than u (e.g., Nurser et al. 1999): Here (4) is the heat contained within seawater with potential temperature less than or equal to u integrated everywhere north of a given latitude f and ›H/›t is its time derivative. In Eq. (3), D diff is the (cold to warm) diffusive temperature flux across the u isotherm north of latitude f from mixing, and represents the heat fluxes upward and out (first term) through the u 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 u isotherm has outcropped [second term; the net turbulent airsea flux Q NET 5 Q latent 1 Q sensible 1 Q net LW 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 u s denotes the SST]. The term represents the advective heat transport of water colder than u northward across latitude f, and c p r 0 Gu represents the advective heat flux up and out across the u isotherm.
Noting that the transport of water with u , u 0 , u 1 du is dc 5 ðð l,z:u,u 0 (l,f,z),u1du yR cosf dz dl ' ›c ›u du , FIG. 1. Derivation of the water-mass transformation equations: (a) volume conservation for waters with potential temperature less than or equal to u, (b) heat conservation for the same region, and (c) water-mass formation rate between isopycnals u and u 1 du 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 228C are shaded brown. Surface and subsurface waters with temperatures between 228 and 238C are shaded purple. Colored arrows indicate water-mass transformation by air-sea heat exchange (yellow) and ocean circulation/ overturning (green) across the 228 and 238C isotherms. Solid and dashed black arrows indicate diapycnal mixing across the 238 and 228C 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.
we can re-express (integrating by parts for the last result) A(u, f) 5 r 0 c p ð u ›c ›u u du 0 5 r 0 c p cu 2 r 0 c p ð u c du 0 (7a) and similarly obtain Eliminating H and A from Eq. (3) using Eqs. (7a) and (7b) and removing the Gu, cu, and (›V/›t)u terms using Eq.
(2) gives Differentiating this by u and again using Eq.
(2) to replace 2›V/›t 1 c by G gives the key result Multiplying Eq. (8) by du makes clear that Eq. (8) represents the heat budget of the volume between the u and u 1 du 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 f, r 0 c p F 5 2 ›D surf ›u ; it then follows ( Fig. 1c) that where the three terms on the RHS of Eq. (10) represent (i) the solar irradiation down through the u 1 du isotherm (a heat source for the layer u , u 0 , u 1 du), (ii) the solar irradiation down through the u isotherm (a heat sink), and (iii) the surface heat flux down through the outcrop u , u s , u 1 du. In terms of F, the water-mass transformation G is then given as 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 ›V/›t and transport c] and F, which is known from the surface fluxes. We can use Eq. (11) to eliminate G in Eq.
(2), arriving at Now (›V/›u)du represents the volume of water with temperature between u and u 1 du, so differentiating Eq. (12) by u gives (13) which describes the increase in volume of water of temperature u in terms of its southern influx [›c(u, f)/›u]du, and the convergence of flow across the u and u 1 du 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.58-308N, 1408E-1208W and 0-300-m depth). The equatorial warm pool is visible to the southwest (u . 278C). From here, isotherms slope up, eastward and northward. Water between 228 and 238C (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 1608E extending eastward to about 1608W, 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 238 , u # 248C (purple shaded surface region extending to the North Pole). This is the larger yellow arrow (F 238C,19.58N ), a loss of volume, which would be compensated in steady state by flow of warmer water across the 238C 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 228 , u # 238C we obtain a smaller value (smaller yellow arrow 5 F 228C,19.58N ), which must be balanced by a diapycnal flow of water across the 228C 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 u # 238C. The lower green arrow represents the volume flux out of the domain for water with u # 228C. 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 u) and a weaker southward flow in a thicker, deeper layer (lower u). As we integrate up from the lowest temperatures, the streamfunction itself is negative, but its vertical gradient is negative at depth/low u and positive at the surface/ high u.
Similar reasoning leads us to conclude that changes to the volume and mixing in the purple region are respectively given by the difference in ›V/›t ›D diff /›u 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 188N 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 238C,188N ; F 228C,188N ). Evidently, the diapycnal fluxes in the region between 188 and 19.58N (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): Here K y is the diffusion coefficient, u is the threedimensional 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 where Q NET (the net surface heat flux) and I s (the surface irradiance) were introduced earlier. Recall that K y consists of a background value, which increases with depth, and a time-varying part, which depends on stratification and velocity shear.
In the surface mixed layer K y becomes large while ›u/›z is small, and the surface flux Q NET 2 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.

a. Ocean temperature response
The global temperature response in the perturbed ensemble experiment is shown in Fig. 2, which is a timedepth 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 2000 s. 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 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). 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 (308S-308N) from the surface to about 100 m and warming in the permanent thermocline throughout the tropics and midlatitudes, but particularly marked in the tropics (208S-208N; 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 wedgeshaped 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 (408S-408N), 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. 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 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 (8C).

b. Energy budget
shortwave (SW) flows. The reflected SW radiation at the TOA (black) is initially reduced by 0.25 W m 22 (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 22 , 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 22 relative to the control during the hiatus but increase rapidly by 2 W m 22 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 diffusivity perturbation causes substantially increased ocean heat uptake of about 10.6 W m 22 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.958C. In response to the reduced diffusivity, the ensemble mean global average SST reduces by about 0.15 K to 17.808C. 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.058C. 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 FIG. 5. Ensemble annual mean anomalous global energy flows (W m 22 ): (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 21 W m 22 , green is 22.5 W m 22 , blue is 23 W m 22 , and cyan is 24 W m 22 . 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). 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.

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 watermass 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 ›V/›t (Fig. 9a) is very small in most temperature classes and latitudes, but › 2 V/›t›u is negative between 08 and 28C (808S-208N), representing a loss of volume. This is counterbalanced by positive › 2 V/›t›u at the same latitude range, between 228 and 08C, representing a gain in volume. Thus, there is cooling in the deep Southern Ocean (;08C) 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 f and temperature u represent all the fluid with temperature # u and north of f. Therefore, as with streamfunction plots in depth-latitude space, a convergence (divergence) of isolines represents accumulation (loss) of volume at a particular u and f. 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) 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 temperaturelatitude plots in Fig. 9 represent complex water pathways in 3D physical space (Fig. 1) Figure 9b shows the overturning streamfunction c 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 ›V/›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 [ 10 6 m 3 s 21 ) in strength, are apparent, centered at about 228C, 108S in the Southern Hemisphere and, slightly cooler, 208C, 108N in the north. These deliver cold water upwelled near the equator (the sharp meridional gradient between 158 and 258C at 08 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 (38-78C), deep waters southward. These upwell at high middle southern latitudes (;408S) 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 ;08C). These currents rise and warm slightly at middle high latitudes in the Northern 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 c (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 408S and 308N below 108C 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 238C range is transformed by air-sea fluxes into warmer waters (up to 308C) mainly in the tropics. Waters in this class are also cooled by air-sea fluxes between the subtropics and midlatitudes to temperatures around 158C. Farther poleward, waters around 58C are transformed into cold deep water below 08C, largely in the Southern Ocean and the North Atlantic.
Mixing, obtained as a residual (›D diff /›u) 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 408S, transforming water from the 158 to the 208C range, and balances the NADW overturning cell at around 508N where waters between 58 and 108C are transformed into the range 08-58C.

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 38S-38N) 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. At the equator, air-sea fluxes (green) cause transformation in temperature classes 178-328C (there is no direct effect on isotherms below 178C because these do not outcrop at the equator). The lines have a negative gradient (formation rate) between 268 and 318C; hence, air-sea fluxes tend to generate water masses in this temperature range [F has a negative sign in Eq. (12)]. Between 178 and 258C, the lines have a positive gradient and air-sea fluxes remove water in these classes. The overturning (red) acts to remove water between 238 and 308C and add water between 108 and 238C. Mixing (blue) differs from the other terms as the associated formation rate has three regimes: between 288 and 328C mixing removes water, between 218 and 288C it adds water, and between 108 and 218C 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 288 and 328C the predominant balance is mixing versus air-sea fluxes (equatorial regime iv), between 258 and 288C it is mixing versus overturning (regime iii), between 218 and 258C it is mixing versus air-sea fluxes (regime ii), and between 108 and 218C it is mixing versus overturning (regime i). Below 158C 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 offequatorial regions (average over 188-128S and 128-188N; the strongest subsurface anomalies in Fig. 4 are between 208S and 208N). 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 288 and 328C and between 208 and 258C, and they remove volume between 258 and 288C. The overturning behaves in an opposite manner than at the equator, adding volume between 258 and 288C and removing it between 208 and 258C. Mixing removes volume between 288 and 328C (off-equatorial regime iv) and adds it between 208 and 258C (regime iii). Between 108 and 208C, overturning and mixing balance, mixing adding volume and overturning removing it (regime ii). The situation is reversed between 58 and 108C 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 38S-38N 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 offequatorial regions (Fig. 11b) where regimes iii and iv FIG. 11. Formation rate (Sv per degree of latitude) (a) at the equator (38S-38N) and (b) for the off-equatorial regions (average of 188-128S and 128-188N) 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. 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 208C isotherm off-equatorial regimes ii-iv are stacked east-west and are linked by the northern STC. 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 278-308C. Conversely there is a small increase in volume at cooler temperatures, confined mainly to 258-278C for the first five years but subsequently extending to cooler temperatures (down to 158C). From year 11 onward there is an abrupt increase in the volume of the warmest waters, 278-308C, and a reduction of volume in the 258-278C range. The volume of water in the 188-258C range decreases more gradually, and in the 158-168C 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 108C, still present at year 20. The increase in volume in the range from 218 to 08C is a reflection of the drift in the control associated with the long-term negative TOA imbalance of 20.2 W m 22 .

CLASSES
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 6 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 258C and a compensating but more modest increase in the 158-258C 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 208 and 258C; however, the volume gained between 158 and 208C is retained in the subsequent surge and new equilibrium period. The 158 and 108C 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 (.258C) between 208S and 208N. Below these waters we note an increased volume of waters between 158 and 258C, extending poleward to 308S/308N, 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 (808-408S) and the northern midlatitudes (208-508N).
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 in the volume of warmer waters. At temperatures between 208 and 258C, the extra volume deposited during the hiatus is lost, but below 208C 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 158C during the hiatus in years 1-10 coincides with the positions of the subtropical cells. The warming and cooling pattern north of 408N 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 (308S-308N) 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 ;258 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 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 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 [›(K y ›u/›z)/›z] similar to regime iv (Fig. 10a), with the further simplification that ocean heat transport divergence = Á (uT) is identically zero. We apply a uniform heating of 20 W m 22 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 W m 22 . 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 Dt 5 864 s; vertical grid spacing Dz 5 1 m). The lower 100 m is treated numerically as a single layer. We assume a constant background diffusion coefficient K y 5 5 3 10 25 m 2 s 21 , set the initial temperature to 208C 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 3 10 25 m 2 s 21 ; 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 W m 22 K 21 . 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 (38S-38N) 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.
The northern off-equatorial region, 128-188N (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).

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: d Experiments with increased ocean diffusivity show Pacificdominated 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.
d 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).
d The climatological thermal balance in the warmest (near surface) waters of the tropical Pacific Ocean is maintained by mixing balancing surface heat flux. In contrast, in the cooler layer below, mixing is balanced by ocean heat transport divergence. Related but slightly different balances prevail in the offequatorial regions. d The overshoot into surge conditions is due to longterm depression of deep (;158C) 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).
d 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: 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.