## Abstract

An ocean mixed layer heat budget methodology is used to investigate the physical processes determining subpolar North Atlantic (SPNA) sea surface temperature (SST) and ocean heat content (OHC) variability on decadal to multidecadal time scales using the state-of-the-art climate model HadGEM3-GC2. New elements include development of an equation for evolution of anomalous SST for interannual and longer time scales in a form analogous to that for OHC, parameterization of the diffusive heat flux at the base of the mixed layer, and analysis of a composite Atlantic meridional overturning circulation (AMOC) event. Contributions to OHC and SST variability from two sources are evaluated: 1) net ocean–atmosphere heat flux and 2) all other processes, including advection, diffusion, and entrainment for SST. Anomalies in OHC tendency propagate anticlockwise around the SPNA on multidecadal time scales with a clear relationship to the phase of the AMOC. AMOC anomalies lead SST tendencies, which in turn lead OHC tendencies in both the eastern and western SPNA. OHC and SST variations in the SPNA on decadal time scales are dominated by AMOC variability because it controls variability of advection, which is shown to be the dominant term in the OHC budget. Lags between OHC and SST are traced to differences between the advection term for OHC and the advection–entrainment term for SST. The new results have implications for interpretation of variations in Atlantic heat uptake in the CMIP6 climate model assessment.

## 1. Introduction

The North Atlantic undergoes variations in sea surface temperature (SST) on multidecadal time scales (e.g., Kerr 2000; Frankcombe et al. 2008; Chylek et al. 2011; Vianna and Menezes 2013), with impacts on the climate of adjacent land areas (e.g., Enfield et al. 2001; Knight et al. 2006; Msadek and Frankignoul 2009; Sutton and Dong 2012; Sutton et al. 2018) and beyond (Lu et al. 2006; Zhang and Delworth 2006). These SST variations are widely referred to as Atlantic multidecadal variability (AMV).

A variety of mechanisms have been proposed to drive AMV, including external forcing by anthropogenic aerosols (Booth et al. 2012), and/or volcanoes (Otterå et al. 2010; Swingedouw et al. 2017), atmospheric forcing (Clement et al. 2015), internal oceanic variability (Sévellec and Fedorov 2013; Gastineau et al. 2018), and coupled ocean–atmosphere modes of variability involving the Atlantic meridional overturning circulation (AMOC; Knight et al. 2005; Ortega et al. 2015). Atmospheric feedbacks are also likely to play a crucial role in setting the AMV pattern (Xie 2009). There is as yet little consensus on the precise mechanism as AMV simulation varies from model to model in both phenomenology and driving processes (Drews and Greatbatch 2017; Muir and Fedorov 2017; Sévellec and Sinha 2018; Sutton et al. 2018).

Observational studies are hindered by the relatively short instrumental record that captures only one or two AMV cycles and lacks information on other variables such as the AMOC. Recent studies have instead utilized AMOC proxies; for example, McCarthy et al. (2015) use a sea level–based indirect proxy of the AMOC to demonstrate a link between the AMOC, OHC in the top 500 m, and AMV from the 1920s to the 2000s.

The link between AMOC and upper ocean OHC is well established in modeling studies (Robson et al. 2012; Zhang 2008; Zhang and Zhang 2015). There is a strong correlation between subtropical AMOC and meridional heat transport (MHT) found in models (Sévellec and Huck 2015; Moat et al. 2016) and observations (Johns et al. 2011). Grist et al. (2010) found, in a model based analysis for 1958–2002, that the subpolar gyre OHC anomaly was more strongly correlated with the ocean heat transport convergence (*r* = 0.75) than with surface fluxes (*r* = 0.5). Similarly, Robson et al. (2014, 2018) and Hodson et al. (2014) found the AMOC and its associated ocean heat transport was the dominant process in the 1990s warming and the 1960s cooling of the subpolar gyre. Likewise, Williams et al. (Williams et al. 2014; R. Williams et al. 2015), using a model that was strongly relaxed to observed temperature and salinity, attributed decadal changes in subpolar gyre OHC to changes in the AMOC.

Whatever the detailed mechanisms and drivers of the AMV, it seems likely that horizontal ocean heat transport convergence and surface fluxes of heat will both play important roles. However the key relationship between changes in oceanic heat transport, OHC, and SST is not well understood, particularly on multidecadal time scales, which is the focus of this paper.

A number of studies have attempted to identify fingerprints of changing AMOC directly on the SST, thus bypassing the need to examine OHC. However, the results from climate models are variable (Roberts et al. 2013; Zhang 2008) and although there is now evidence of a similar pattern associated with the limited duration observational record (Smeed et al. 2018), without a verified mechanism it is difficult to be confident in these fingerprints.

A more rigorous approach was adopted by Buckley et al. (2014) where interannual heat content was evaluated over the depth of the monthly maximum climatological mixed layer (i.e., the portion of the upper ocean in contact with the atmosphere). Using the ECCO state estimate for the period 1992–2010 they estimate that 70% of the variability in mixed layer heat content is explained by local forcing (i.e., air sea heat fluxes and Ekman convergence) and only 30% due to advection over large parts of the North Atlantic. Their use of the monthly maximum mixed layer was an improvement over previous studies which employ a spatially constant depth horizon. However, because of the length of the simulation they were unable to address multidecadal time scales.

The Buckley et al. (2014) approach was extended to the global domain by Roberts et al. (2017), who used a similar theoretical framework, including a spatially variable maximum mixed layer depth (MLD) to differentiate the near-surface layer in contact with the atmosphere from the rest of the ocean. Unlike Buckley et al. (2014) they used observationally based gridded OHC products and surface fluxes from atmospheric reanalyses with a Kalman filter–based method to obtain an estimated heat budget with error bounds for both the mixed layer and the rest of the ocean, evaluating ocean heat transport convergence as a residual. Their results indicated that on interannual time scales there are extensive regions (equator and western boundary currents and the Antarctic Circumpolar Current) where ocean heat transport convergence dominates the OHC variability of the mixed layer and over large parts of the rest of the ocean both ocean heat transport convergence and surface heat fluxes are important. This contrasts with the full-depth OHC, which on these time scales is dominated by ocean heat transport convergence.

In this paper we consider temperature changes in the mixed layer, taking account of its time-varying depth using the SST evolution equation described by Stevenson and Niiler (1983), paying particular attention to the diffusive flux at the base of the mixed layer.

We address the following questions using a state-of-the-art coupled climate model that we demonstrate has realistic Atlantic multidecadal variability:

What controls the multidecadal evolution of full depth OHC in the subpolar North Atlantic? What are the respective roles of ocean surface fluxes versus internal ocean processes? Is there a difference between the deep convection regions to the west and the region farther east that is more influenced by the North Atlantic Current?

What controls the multidecadal evolution of SST in the subpolar North Atlantic? Are the respective roles of surface fluxes and internal ocean processes similar and if not how do they differ?

What is the relationship between changes in the deep (sub mixed layer) OHC and SST? How and why are the forcing terms different?

How do both deep OHC and SST depend on the AMOC?

We focus on the subpolar North Atlantic (SPNA) as the AMV spatial pattern is strongly concentrated in this region (e.g., Sutton et al. 2018). In contrast, subtropical AMV is thought to be caused by relatively rapid (months to a few years) adjustment of the subtropical ocean to changes in the subpolar gyre via boundary waves (Johnson and Marshall 2002), or by atmospheric feedbacks to SPNA variability (Sutton et al. 2018).

The paper is organized as follows. We use a rigorous theoretical framework for comparing OHC and SST variability in section 2. Details of the model configurations and methodology are given in section 3. The results are presented in section 4 and conclusions in section 5.

## 2. Theory

### a. Full depth ocean heat content

We define the full depth ocean heat content per unit area (Θ_{FD}) as

where *λ* and *φ* are longitude and latitude, respectively; *t* is time, *z* is depth (increasing upward), *θ* is potential temperature, *H* is local water depth; and *ρ*_{0} and *C*_{p} are a reference seawater density and specific heat capacity respectively.

Changes in Θ_{FD} at any given location can be caused by heating/cooling at the air–sea interface (*Q*_{NET}) or by horizontal advection and/or diffusion (considered here as one term, *R*_{FD}) resulting in a simple evolution equation:

Observationally, ∂Θ_{FD}/∂*t* can be estimated from Eq. (1) using ocean temperature measurements, *Q*_{NET} using atmospheric reanalysis, and hence *R*_{FD} as a residual, although for climate-relevant time and space scales each term would carry considerable uncertainty. Alternatively, a heat-conserving climate model simulation can provide exact ∂Θ_{FD}/∂*t* and *Q*_{NET}, with *R*_{FD} again evaluated as a residual, for comparison with observed estimates. In principle, in a climate model *R*_{FD} could be calculated directly rather than as a residual, but in practice this is rather difficult because diffusive as well as advective lateral transport convergences would be required, and these were not stored for the simulation used in this study.

Equations (1) and (2) could equally be evaluated over different depth horizons if *H*(*λ*, *φ*) in Eq. (1) is replaced by a fixed depth taking bottom topography into account. We examine the sensitivity of our results to choice of depth horizon in section 4e.

### b. Sea surface temperature

We employ the mixed layer temperature evolution equation derived by Stevenson and Niiler (1983)

where *h* is the mixed layer depth, **v**_{a} is the vertical average within the mixed layer of the horizontal velocity vector, $v^$ and $T^$ respectively are deviations of the horizontal velocity and temperature from their vertically averaged values, and **v**_{−h}, *T*_{−h}, *w*_{−h}, and *Q*_{−h} are the horizontal velocity, temperature, vertical velocity, and diffusive heat flux at the base of the mixed layer.

Neglecting horizontal diffusion, changes in the ocean temperature averaged over the mixed layer, *T*_{a}, at any given location can be caused either by heating/cooling at the air–sea interface, *Q*_{NET}, or by horizontal advection, vertical advection/diffusion of heat, and entrainment/detrainment of fluid into or out of the mixed layer.

Defining *ξ* to be the sea surface temperature and substituting *T*_{a} = *ξ* − (*ξ* − *T*_{a}), Eq. (3) can be recast into a simpler form analogous to Eq. (2):

where

represents the aggregated effect of all internal ocean processes plus an error term, ∂(*ξ* − *T*_{a})/∂*t*, which indicates how well the SST tendency, ∂*ξ*/∂*t*, approximates the depth averaged temperature tendency, ∂*T*_{a}/∂*t*. We focus on SST because the AMV index, the main motivation of our study, is defined in terms of SST. Also, we would like to apply our method to observations in the future and *T*_{a} is not routinely available from observations, partly because mixed layer depth is not known with sufficient accuracy whereas there are many high-quality SST datasets available. As for Eq. (2), each of the terms in Eq. (4), with the exception of *R*_{ML}, can be diagnosed from climate model output, or from observations as long as the MLD is available as a function of time. However, observational datasets of the MLD are limited to monthly mean climatologies (e.g., de Boyer Montégut et al. 2004). The rate of change ∂*ξ*/∂*t* can be estimated from observed sea surface temperature.

Once ∂*ξ*/∂*t* and *Q*_{NET}/*h* have been calculated, *R*_{ML}/*h* and *R*_{ML} can be evaluated as a residual from observations (with associated observational uncertainty), or exactly from a heat-conserving climate model.

### c. Anomaly formulation

As we are interested in decadal variations of heat content and SST, we recast Eqs. (1) and (4) in terms of anomalies from long-term mean quantities. For the heat content, this is straightforward; we average Eq. (2) over sufficiently long time scales, the time derivative tends to zero, and we obtain

where the asterisk denotes a deviation from the long-term mean value. We will refer to $RFD*$ as the “advection” term (because lateral diffusion can be assumed to be small) and to $QNET*$ as the “surface flux” term. Note that it is not always true that averaging over longer periods will result in an exact balance between $QNET\xaf$ and $RFD\xaf$. However, for our analysis it is not a necessary condition. The only requirement is that mean values are removed from ∂Θ_{FD}/∂*t*, *Q*_{NET}, and *R*_{FD}. This is the equivalent to detrending *T* and centering *Q*_{NET} and *R*_{FD} on zero for the period of interest.

A similar procedure can be adopted for the SST using Eq. (4):

leading to

We will refer to [*Q*_{NET}/*ρ*_{0}*C*_{P}*h*]* as the “unadjusted surface flux term” for SST and to [*R*_{ML}/*ρ*_{0}*C*_{P}*h*]* as the “unadjusted advection-diffusion-entrainment” term for reasons that will become clear shortly; however, for comparison with Eq. (7) this formulation is not very convenient. Instead we return to Eq. (4) and taking correlations between *h** and *ξ** into account (see appendix A) we obtain the following equation for the SST anomaly *ξ**:

Note that the denominator of the terms on the right-hand side of Eq. (10) is the mean mixed layer depth $h\xaf$, not the instantaneous value *h*, as in Eq. (4). Also $\u211cML$ is a different residual to *R*_{FD}.

The first term on the RHS of Eq. (10) represents “external” forcing of the SST by surface fluxes, while the second term represents trends due to “internal” processes in the ocean. These are analogous to $QNET*$ and $RFD*$ in Eq. (7). The reasons for differing temporal evolution of SST and OHC are contained in Eqs. (7) and (10), in particular the difference between $RFD*$ and $\u211cML*$. At any given point, if $RFD*$ and $\u211cML*$ were identical, $\rho 0CPh\xaf\xi *$ (and hence *ξ**) would have the same temporal evolution as $\Theta FD*$. Hence, we analyze the relationships between these terms later in order to understand differences between the time evolution of SST and OHC in the SPNA. We will refer to $QNET*/\rho 0CPh\xaf$ as the “surface flux” term for SST and to $\u211cML*/\rho 0CPh\xaf$ as the “advection-entrainment” term.

### d. Parameterization of diffusive vertical heat flux

We will find that the terms on the right-hand side of Eq. (10) are generally of opposite sign and ∂*ξ**/∂*t* is much smaller in magnitude than either, which makes it difficult to identify which term is most important. This is because the diffusive vertical heat flux, *Q*_{−h}, can be of the same order of magnitude as *Q*_{NET} in Eq. (3). We can therefore reformulate Eq. (9) as

where $\mathbb{R}ML$ is yet another residual representing advection and entrainment, but excluding vertical diffusion. The term *Q*_{−h} is generally parameterized in models as a function of vertical temperature gradient *K*(*∂T*/*∂z*), where *K* is a time-variable diffusion coefficient. Here we adopt an even simpler approach and crudely parameterize it as a constant proportion of the surface heat flux anomaly $Q\u2212h*=\lambda QNET*$, where *λ* is a constant. This gives an alternative formulation for the SST tendency:

Our motivation in this paper is to relate the SST variation to the full-depth OHC variations, so we select a measure that will maximize the relationship between them. Therefore we determine the value of *λ* by requiring the strongest correlation between $RFD*$ and $\mathbb{R}ML*$ (see section 4d and appendix B). We will refer to $\u2061(1\u2212\lambda )QNET*/\rho 0CPh\xaf$ as the “adjusted surface fluxes” term for SST and to $\mathbb{R}ML*/\rho 0CPh\xaf$ as the “adjusted advection-entrainment” term. We note that our use of the coefficient *λ* is an empirical approach: we do find large correlations between $RFD*$ and $\mathbb{R}ML*$ in the SPNA (up to 0.87 in the eastern SPNA and 0.63 in the western SPNA; appendix B). However further investigation, beyond this paper, is required to understand the full significance of *λ*.

## 3. Model description and analysis procedure

### a. HadGEM3-GC2 coupled climate model

We analyze output from a 300-yr preindustrial control simulation HadGEM3-GC2 (K. Williams et al. 2015), a high-resolution version of the Met Office HadGEM3 climate model, including ocean, atmosphere, sea ice, and land surface components. The ocean configuration is the Global Ocean version 5.0 (Megann et al. 2014) of the v3.4 NEMO model (Madec 2008), which uses the ORCA025 tripolar grid (~0.25° horizontal resolution) and 75 vertical levels. The sea ice component, also on the ORCA025 grid, is version 4.1 of the Los Alamos Sea Ice Model (CICE; Hunke and Lipscomb 2010), which includes five sea ice thickness categories and has improved representation of Arctic sea ice concentration with respect to previous configurations (Rae et al. 2015).

The atmosphere component is the Global Atmosphere version 6.0 of the Met Office Unified Model (UM; Walters et al. 2011), with a horizontal resolution of N216 (~60 km at midlatitudes) and 85 levels in the vertical. The land surface model is the Global Land version 6.0 of the Joint U.K. Land Environment Simulator (JULES; Best et al. 2011), which shares the same grid as the atmospheric component.

This control simulation has been employed in many studies to examine a variety of climate system processes. For example, the model has been used to examine mechanisms of decadal variability in the Labrador Sea (Ortega et al. 2017), predictions of the winter NAO (Scaife et al. 2014; Dunstone et al. 2016), and climatic trends in the North Atlantic (Robson et al. 2016).

### b. Analysis procedure

Equations (1)–(12) were evaluated from the GC2 climate model simulation (K. Williams et al. 2015) using monthly mean potential temperature, MLD (defined as the depth at which the potential density referenced to the surface differs from the surface density by 0.01 kg m^{−3}), and mean net surface heat flux.

For each model year we take each month and calculate the average tendency terms for SST and OHC for the 1-yr period from that month to the same month in the next year (January 2294–January 2295, February 2294–February 2295, … , December 2294–December 2295). We then calculate the mean of these 12 averaged tendency terms to obtain a consolidated tendency term representative of the entire year. With this approach an exact heat budget for the annual mean OHC or SST anomaly is obtained. A constant value of *ρ*_{0}*C*_{P} = 4.1 × 10^{6} J m^{−3} K^{−1} was used throughout.

The AMOC at 26° and 50°N was taken from the annual mean overturning streamfunction output as a standard model diagnostic. The AMV index was calculated as the annual mean SST averaged over the North Atlantic (0° to 65°N, 75° to 7.5°W) minus the annual mean global SST normalized by the standard deviation (after Sutton et al. 2018):

where the overbar represents a spatial average, angled brackets represent a time average, and the standard deviation *σ* is taken over the 300-yr simulation.

All variables are filtered to retain periods of 10 years and longer using an 11-point Parzen filter for annual means, or a 121-point filter for monthly means (Press 1986). The results were essentially the same using a running mean filter.

## 4. Results

### a. Mean OHC and SST tendency terms

Over long time scales, the mean OHC tendency is very small and surface fluxes balance advection as in Eq. (6), hence it is sufficient to examine just one of these latter terms in order to understand the mean state. The HadGEM3-GC2 300-yr mean *Q*_{NET} is shown in Fig. 1a. The net heat flux term shows cooling in the Gulf Stream region and SPNA (north of a line connecting Florida with the Bay of Biscay) and warming in the subtropics (south of this line). The cooling is considerably weaker in the central SPNA, and there is a strong region of warming on the shelf region of the Grand Banks. The warming in the subtropics is enhanced toward the shelf-slope regions bordering Africa and South America. Equation (6) indicates that advection has a mean pattern opposite to the surface heat flux term with cooling in the subtropics and warming in the subpolar regions.

Thus, as expected, the model shows warming in the subtropics and cooling in the subpolar regions due to differential radiative heating. The ocean circulation (mainly the AMOC in the North Atlantic) redistributes the excess heat in the tropics toward the pole.

The HadGEM3-GC2 300-yr mean SST tendency due to surface fluxes in the North Atlantic, the first term in Eq. (8), is shown in Fig. 1b. Surface fluxes introduce a warming SST tendency everywhere with the exception of the western boundary regions and some small isolated regions in the tropics, and the Greenland and Labrador Seas. In the Gulf Stream Extension, North Atlantic Current, and subpolar gyre regions the sign is opposite to the effect of surface heat fluxes on the OHC (cf. Fig. 1a). Also the pattern is different, with maximum values over the Grand Banks shelf region, in the subtropical gyre, and in the western subpolar region. The prevailing positive tendencies occur because of the MLD factor *h* in the denominator of the *Q*_{NET}/*ρ*_{0}*C*_{P}*h* term in Eq. (8), which weights the annual mean toward the summer months when the MLD is shallowest and the ocean experiences heat gain from surface fluxes. Advection-diffusion-entrainment opposes the warming effect of surface fluxes and hence is negative in most locations.

The result that in most of the North Atlantic north of 30°N, surface fluxes impose a negative trend on the annual mean full depth heat content while also imposing a positive trend on the annual mean SST is somewhat counterintuitive and bears further explanation. As an illustration, Fig. 1c displays the MLD over the model year 2295 at 55°N, 28°W. In January, the MLD is 300 m. It deepens to a maximum of 400 m in February before shallowing over spring (March–May) to a minimum of about 20 m in June. Over summer (June–August) the mixed layer remains very shallow but during the autumn it deepens, reaching in excess of 100 m in December. For comparison the maximum winter MLD over the 300-yr simulation (482.5 m) is shown as a solid line. Also marked are the 100- and 200-m depth levels. Evidently, use of a temporally fixed depth to characterize the mixed layer (e.g., Buckley et al. 2014; Roberts et al. 2017), while mathematically simpler, is problematic. Heat content in such a fixed layer is not simply related to SST in any season.

Surface heat flux, *Q*_{NET}, for each month of the year is plotted in Fig. 1d (blue). There is strong (turbulent) heat loss from ocean to atmosphere between January and March and again between October and December. In summer, between May and August, the ocean gains heat due to increased insolation. At this example location, the seasonal variation of the net surface heat flux is ±200 W m^{−2}.

The red line in Fig. 1d represents the accumulated net surface heat flux (i.e., the accumulated sum of the values plotted in blue). The accumulated heat flux remains negative over the whole year, indicating that winter heat loss strongly outweighs summer heat gain. Hence in the annual mean, surface heat flux tends to reduce OHC and a negative value is found in Fig. 1a at this location.

The surface flux–related forcing term for SST, *Q*_{NET}/*h*, is plotted in red in Fig. 1e. The high values of *h* in winter, spring and autumn compared to summer (up to 20 times higher) result in much smaller values of *Q*_{NET}/*h* in these seasons so the accumulated value of *Q*_{NET}/*h* (red) is strongly positive from May to the end of the year and a positive value is found in the corresponding location in Fig. 1b.

### b. Simulated AMV variability

A common hypothesis for the observed temporal AMV variability is heat redistribution by the AMOC. While changes in the AMOC and associated changes in horizontal heat transport divergence can potentially affect full depth OHC, whether and by what mechanism changes in full-depth OHC are translated to changes in SST are not clear. In this section, we first examine the relationship between the AMOC and AMV in the HadGEM2-GC2 simulation, then use our theoretical framework to obtain insights into the mechanisms.

Figure 2a shows the AMV index calculated from annual mean model output, together with the AMOC anomaly at 26°N (Fig. 2b), and 50°N (Fig. 2c), with respect to its long-term mean, low-pass filtered with a cutoff period of 10 yr. The AMV index shows multidecadal variability reminiscent of the observations and the time scale of the variability (~50 yr) is within the range estimated from observations and multimodel studies (20–70 yr). There are four large AMOC excursions in the simulation period (Figs. 2b,c) and these are matched with large AMV fluctuations. The spatial pattern associated with the AMV (regression coefficient of the linear correlation of SST with the AMV index) is shown in Fig. 2d. The pattern approximately matches that obtained from observations (e.g., Sutton et al. 2018; see also Kushnir 1994) but the region of low regression in the western subtropics (between Florida and Cape Hatteras) is larger than that seen in observations and in addition the Greenland Sea shows the opposite sign regression coefficient. However, the HadGEM3-GC2 control simulation has fixed atmospheric aerosol and CO_{2} concentrations, whereas the real-world AMV may be influenced by changing concentrations of anthropogenic aerosols or greenhouse gases. Hence, even if the model was perfect, we might not expect or demand exact agreement. On the other hand, the current generation of climate models shows a range of AMV time scales and spatial patterns, hence some of the results presented in this study may be model dependent and it will be important to compare them across a range of models in future.

Correlation analysis shows a lagged relationship between AMOC and AMV with a maximum correlation coefficient of 0.56 (26°N) and 0.52 (50°N) with the AMOC leading by ~5 and ~9 years (Fig. 3a). The thicker black and red lines indicate significance at the 95% level. Both time series were detrended and autocorrelations were considered in determining the degrees of freedom for significance testing (Emery and Thomson 1997). Although significant correlations are found, they do not account for all the AMV variance and many other processes could contribute to the AMV variability including subpolar gyre variability independent of the AMOC, atmospheric teleconnections from the tropics, and variability of the Northern Hemisphere cryosphere including sea ice and snow cover.

The time series of the AMOC at 26°N (Fig. 3b) and 50°N (Fig. 3c) are divided into events (labeled A to D) where each event spans a full AMOC cycle. We subdivide each event into four phases corresponding to decreasing and increasing AMOC during periods of negative and positive AMOC anomalies respectively. Thus, each event has a full cycle of the AMOC during which the AMOC anomaly reduces to a minimum (phase 1; red), increases from the minimum to zero (phase 2; blue), increases to a maximum value (phase 3; cyan), and then decreases to zero (phase 4; magenta). The year numbers of these events (based on the AMOC excursions, not the matching AMV excursions) are listed in Tables 1 and 2. The duration of the events varies between 12 and 65 years and individual phases vary from 3 to 26 years.

In the next section, we will investigate the processes controlling the OHC and SST trends in the different phases of the AMOC cycle. We will concentrate on SST and OHC changes during the four events, focusing both on the full time series, and on a composite of all four events.

### c. OHC trends during different phases of the AMOC

Figures 4a–d shows OHC trend composites based on the AMOC at 26°N (upper panels) in the North Atlantic for each phase in turn of a composite of all four events A–D (annual mean trends averaged over the duration of each phase of all four events; Table 1 lists the model years included in each phase). During phase 1 (AMOC anomaly < 0 and reducing) there is a negative heat content trend in the SPNA coupled with increasing OHC in the subtropical gyre (STG) and in the Nordic seas (Fig. 4a). There is a dipole pattern in the intergyre region (Cape Hatteras to the Bay of Biscay) with positive trends in the west and negative in the east. In phases 2–4 we see positive OHC trends spreading first into the eastern and northern SPNA (Fig. 4b) and later into the western SPNA (Figs. 4c,d). Negative trends appear in the western part of the intergyre region in phase 2, but there is a return to positive trends in phases 3 and 4. The STG (south of 40°N) shows somewhat complicated behavior, with mainly positive trends in phases 1 and 4, and opposite signed north–south dipoles in phases 2 and 3. The Nordic seas as a whole vary coherently, with OHC increasing in phases 1 and 4, and decreasing in phases 2 and 3. The OHC trend composites based on the AMOC at 50°N (Figs. 4e–h) show strong similarity with those based on the AMOC at 26°N. Phase 1 at 26°N (Fig. 4a) and at 50°N (Fig. 4e) look very similar, for example, as do phases 2–4. Of particular note is the fact that in the SPNA, in phase 4, when the AMOC is reducing, the OHC shows a warming tendency.

The SST trend composites based on the AMOC at 26°N show some similarities to the corresponding OHC composites in the SPNA, intergyre, and STG regions, but also substantial differences (cf. Figs. 4a–d with Figs. 4i–l). For example, in phase 3 the patterns are broadly similar with extensive warming over the whole SPNA and the eastern SPNA (Figs. 4c,k), whereas in phase 4 the western SPNA has strongly increasing OHC (Fig. 4d) but the SST is weakly increasing (Fig. 4l). In contrast to the OHC trends, the peak SST warming occurs in phases 2 and 3, not in phase 4. In the SPNA, SST trend composites based on the AMOC at 50°N show a phase lag compared to those based on the AMOC at 26°N. Phases 2, 3, and 4 of the 50°N-based composites (Figs. 4n–p) are very similar to phases 1, 2, and 3 respectively of the 26°N-based composites (Figs. 4i–k). The SST phases at 26°N bear some similarity to the observation based normalized SST trends presented in Fig. 2 of Caesar et al. (2018). The sequence established for the composite event is essentially seen in the individual events (not shown).

Figure 5 shows the two terms, $QNET*$ and $RFD*$, that determine the full-depth heat content tendency. As the AMOC increases from a minimum (Fig. 5a), there is a positive heat flux anomaly in the northwestern subtropical gyre with the exception of the Gulf Stream, which has a negative surface heat flux signature. Elsewhere in the subtropical and subpolar gyres, the heat fluxes are rather weak except over the East and West Greenland boundary current and in the Norwegian Sea where there is anomalous heat input. Subsequently there is widespread heat uptake in both SPNA and subtropical regions (Fig. 5b). Phases 3 and 4 (Figs. 5c,d) then reverse the sequence, with phase 3 being a negative version of phase 1 and phase 4 a negative version of phase 2. It is remarkable that there is strong heat gain (loss) due to surface fluxes in the SPNA in phase 2 (phase 4) when the AMOC is increasing (decreasing). The composites based on the AMOC at 50°N (Figs. 5e–h) are very similar in pattern to those based on the AMOC at 26°N, the main differences being in the magnitude of the anomalous fluxes.

In phases 1 and 2, the advection anomaly term, $RFD*$, is very similar to the surface heat flux anomalies but opposite in sign, and the net tendency is a small residual between the terms (cf. Figs. 5i,j with Figs. 5a,b). Hence it is difficult to pick out by eye which term is the larger. However, in phases 3 and 4 we see larger differences in the patterns and in the SPNA; in particular, it is possible to discern which term is dominant. In phase 3 in the western SPNA, surface fluxes appear to be the larger term whereas in the eastern SPNA advection dominates (cf. Fig. 5k with Fig. 5c). In phase 4 on the other hand it appears that advection is the dominant term throughout the SPNA. As with the surface fluxes, composites of $RFD*$ based on the AMOC at 50°N differ slightly in magnitude from those based on the AMOC at 26°N, but the spatial patterns obtained are very similar.

### d. SST trends during different phases of the AMOC

We now examine the contributions to the net SST tendency shown in Figs. 4i–p, focusing on the advection-diffusion-entrainment related term [*R*_{ML}/*ρ*_{0}*C*_{P}*h*]*. From Eq. (9) we plot [*R*_{ML}/*ρ*_{0}*C*_{P}*h*]* for each phase in Figs. 6a–d. The term shows interesting spatial structure particularly around the Labrador and Irminger Seas (areas of deep convection in the model; Ortega et al. 2017). The Gulf Stream and its extension in particular shows systematic changes in sign and magnitude with a warming signal in phases 1 and 4 and a cooling signal in phases 2 and 3 reminiscent of the advection term in the OHC equation (Figs. 5e–h) The surface flux related term [*Q*_{NET}/*ρ*_{0}*C*_{P}*h*]* (not shown) is essentially similar in pattern, but of negative amplitude. These two terms are much larger than the net tendency, which is the residual between two very large and opposing terms, hence this decomposition yields little insight into the relative role of each term.

Turning to the anomaly formulation in Eq. (10) we now plot $\u211cML*/\rho 0CPh\xaf$ for each phase (Figs. 6e–h). We discern a different temporal evolution, without a strong signal in the convection regions but not so clearly reminiscent of the OHC advection especially in the Arctic and the East and West Greenland currents where the shallow mixed layer results in uniformly large values. The magnitude of the term (±0.5 W m^{−3}) is still much larger than the net tendency term (Figs. 4i–p) and hence $\u211cML*/\rho 0CPh\xaf$ and $QNET*/\rho 0CPh\xaf$ still nearly cancel. Thus we still obtain little insight into the controlling process.

Finally we turn to Eq. (12) and evaluate $\mathbb{R}ML*/\rho 0CPh\xaf$ using *λ* = 0.99 (this choice of *λ* is justified in appendix B; we note that it is obtained by searching for the maximum correlation between $\mathbb{R}ML*/\rho 0CPh\xaf$ and *R*_{FD}). By doing this we obtain magnitudes that are of the same order of magnitude as the net tendency. We therefore adopt this decomposition in the following analysis.

### e. Eastern and western subpolar North Atlantic

There is a tendency for both OHC and SST to show different responses in the western compared to the eastern SPNA (see Fig. 4 for region definitions). For example, for the composites based on the AMOC at 26°N (Fig. 4), in phase 1 there is a more negative OHC tendency in the eastern SPNA than in the western SPNA; in phase 2 the tendencies are of opposite sign; and in phase 3 there is a stronger positive tendency in the eastern SPNA than the west. Accordingly, we investigate the spatially averaged response in each region separately.

OHC and SST spatially averaged over the eastern and western SPNA and the AMOC at 26° and 50°N are plotted in Fig. 7. Both regions show a lagged relationship between the AMOC at both latitudes (black, magenta; Fig. 7c) and the OHC (red; Figs. 7a,b). At 26°N the AMOC leads western SPNA OHC by 15 years and at 50°N by 18 years (see Table 3). The corresponding lead times for the eastern SPNA are 10 and 12 years. This is consistent with our earlier finding that OHC tendencies tend to propagate anticlockwise around the SPNA (Fig. 4). The SST (blue, Fig. 7) is also related to the AMOC, but the lag is smaller compared to the OHC (5 and 7 years at 26° and 50°N respectively) and it does not vary between the eastern and western SPNA.

To explore the possibility that the OHC variability may depend on the depth horizon over which it is evaluated, we evaluate the OHC [Eq. (1)] from the surface to depth horizons of 100, 200, 500, and 1000 m and the full ocean depth [Fig. 8; for this figure (only) we use monthly data in order to accurately characterize the lags between different depth horizons]. In both eastern and western SPNA (Figs. 8a,b), the variability of the OHC is qualitatively similar no matter which depth horizon is employed: correlations of the OHC at 1000 m with the OHC at shallower depths yields *r*^{2} values between 0.63 (100 m) to 0.94 (500 m) in the west and 0.83 (100 m) to 0.93 east (500 m). However, the variability for deeper depth horizons lags with respect to shallower ones (Figs. 8c,d); this is particularly marked in the western SPNA where the correlation also drops more rapidly with depth. Nevertheless, a robust result is that SST leads heat content, irrespective of the depth horizon to which it is evaluated (Figs. 8c,d). This is quantified in Table 4, which shows the maximum correlation, and the lag at which the maximum occurs, of each heat content evaluation with the SST. In the western SPNA, the lag narrows from 45 months for full depth to only 3 months for 100-m depth, but a lag always remains. In the east, the SST and the OHC become almost simultaneous for depths shallower than 200 m and the correlation becomes very close to unity. In summary, for depth horizons greater than 200 m the maximum correlation between OHC and SST never reaches unity and a substantial lag (1.5 to 3.5 yr) occurs.

### f. Balance between surface fluxes and advection

#### 1) OHC

Having discussed the variation in OHC we now examine the processes controlling its rate of change via Eq. (7). Figures 9a and 9b show $QNET*$, and $RFD*$ averaged over the western and eastern SPNA respectively while Figs. 9c and 9d show similar plots for $\u2202\Theta FD*/\u2202t$. The rate of change of OHC (red) displays decadal time scale shifts from positive to negative values, during which OHC rises and falls respectively. The events noted earlier (Table 1) are visible as longer than average periods of increasing OHC (e.g., years 2120–60, 2290–2330, and 2390–2410). This rate of change is caused by the interplay between the surface fluxes (*Q*_{NET} in black) and advection (*R*_{FD} in blue), which tend to oppose to each other, but not always.

The term with the larger absolute magnitude will drive the sign of the OHC tendency. If the other term is of opposite sign then it will act as a brake whereas if it is of the same sign then the two terms act in concert. For example in year 2240 the absolute magnitude of the surface heat flux ($QNET*$) is larger than that of advection, the two terms act in concert, the net rate of change is positive, and the anomalous heat content rises. In contrast, in year 2400, the absolute magnitude of the surface heat flux is less than that of advection and it is of opposite sign, the rate of change is positive, and heat content rises with advection driving and heat fluxes acting as a brake.

In Fig. 9a the surface flux term often leads advection by a few years, *r* = 0.5 at 6.5 years, implying that in the western SPNA surface fluxes control the evolution of the full-depth OHC. However, *Q*_{NET} and *R*_{FD} are significantly anticorrelated and correlated respectively with the AMOC at 26°N with the AMOC leading or simultaneous (Table 3). This implies that it is the AMOC that is the main driver of the heat content. Further support for this conclusion will be given in section 4f, which considers the time evolution of a composite AMOC cycle.

In the eastern SPNA, the opposite pattern occurs (Fig. 9b). First, the decadal variability of advection (*R*_{FD} in blue), 6.1 W m^{−2}, is much larger in magnitude than that of surface fluxes, 3.9 W m^{−2}, unlike the western SPNA where the variability is of roughly equal amplitude (both ~4.3 W m^{−2}). In addition, advection tends to lead surface fluxes by a few years, *r* = 0.3 at 11 years (disregarding a peak at 2.5 years, which is statistically insignificant), suggesting that advection is the controlling process in this region. Once again, the AMOC is significantly related to both terms (Table 3).

#### 2) SST

Moving on to the processes controlling the SST, we have already noted that in order to make progress we need to use Eq. (12) with a parameterized heat flux (*Q*_{−h}) at the base of the mixed layer. This is further illustrated by Fig. 10a, which shows the relative contributions of surface fluxes $\u2061([QNET/\rho 0CPh]*)$ and other processes $\u2061([RML/\rho 0CPh]*)$ in the western SPNA from Eq. (9): these are very different compared to the OHC terms *Q*_{NET} and *R*_{FD} (note that the net tendency terms for SST are plotted in Figs. 10g and 10h). There is no discernible lag between the two terms; they are coincident in time and are of opposite sign, and very small differences in magnitude between them determine the sign of the rate of change of SST. Similar considerations apply to the eastern SPNA (Fig. 10b).

Using Eq. (10), we again find a very high degree of compensation between the surface flux and advection terms (Figs. 10c,d), although now the surface flux–related term for the SST, $QNET*/\rho 0CPh\xaf$, has almost exactly the same variation as for the surface flux term for the OHC, *Q*_{NET} (cf. the black line in Fig. 10c with the black line in Fig. 9a). The small differences arise because here we are applying a spatial average and $h\xaf$ varies spatially, though not with time.

Finally we parameterize the diffusive heat flux Eqs. (11) and (12) at the base of the mixed layer. The simple parameterization results in a separation of the surface heat flux and advection-entrainment term (Figs. 10e,f). This decomposition allows us to draw similar conclusions for the SST as we drew for the OHC, namely that in the western SPNA both surface heat flux and advection-diffusion-entrainment play major roles in setting the net SST tendency. By contrast in the eastern SPNA, the advection-entrainment term is the clear driver of SST variations on decadal time scales.

#### 3) Relationship between OHC and SST tendency terms

A strong relationship emerges between the rates of change of full depth OHC (∂OHC/∂*t*) and rates of change of SST (∂SST/∂*t*) (Fig. 11a). Maximum positive correlations are found at lags of 18 months (west) and 3 months (east). As well as these positive correlations, negative correlations are found when the ∂OHC/∂*t* leads ∂SST/∂*t* by 63 months (west) and 67 months (east).

As mentioned in the previous subsection, the adjusted surface flux–related term for the SST, $\u2061(1\u2212\lambda )QNET*/\rho 0CPh\xaf$, has very similar variation as for the surface flux term for the OHC, *Q*_{NET}, with small differences arising because $h\xaf$ varies spatially. This similarity is illustrated in Fig. 11b. The two surface flux–related terms vary simultaneously and the maximum correlation is unity.

By parameterizing the heat flux at the base of the mixed layer, we obtain strong lagged correlations of the advection-entrainment term, $\mathbb{R}ML*/\rho 0CPh\xaf$ with $RFD*$. In the western SPNA $\mathbb{R}ML*/\rho 0CPh\xaf$ leads $RFD*$ by ~3 years (*r* = 0.62) while in the eastern SPNA $\mathbb{R}ML*/\rho 0CPh\xaf$ is almost simultaneous with $RFD*$ (*r* = 0.78). Additionally, both terms, $\mathbb{R}ML*/\rho 0CPh\xaf$ and $RFD*$, have a significant correlation with the AMOC at 26°N in both regions of the SPNA (Table 3).

### g. Drivers of net tendencies in OHC and SST

We obtain further insights into the controls on SST and OHC variation by forming a composite AMOC anomaly cycle based on all four individual events. To do this we take each phase in turn and assign identical timings for the start and end points of the phase. Thus for phase 1 the start year of each event is set to time zero and the end year is set to 2*π*. For example, phase 3 of event C spans years 2275–88, including 14 years, whereas event D spans 2385–95 for phase 1, a total of 11 years. Thus both 2275 and 2385 are assigned a time of *π* and 2288 and 2395 are assigned a time of 3*π*/2 and all intermediate values are interpolated onto a regular grid with spacing *π*/50. In this way all four events and all four phases can be stretched onto a common timeframe and averaged to form a composite AMOC anomaly at 26°N and associated anomalies of SST (*ξ**) and OHC ($\theta FD*$) in the western and eastern SPNA (plotted as a function of phase, *φ*, in Figs. 12a and 12b). By our definition, the composite AMOC anomaly (black line) is zero at phase values *φ* = 0 and *φ* = 2*π*. In between these values the AMOC is negative between *φ* = 0 and *φ* = *π* and positive between *φ* = *π* and *φ* = 2*π*. Local extrema occur near *φ* = *π*/2 and *φ* = 3*π*/2 and the anomaly is near zero at *φ* = *π*. The minimum value is ~−0.9 Sv (1 Sv ≡ 10^{6} m^{3} s^{−1}) and the maximum slightly larger at ~1.0 Sv. SST anomaly (dark blue) closely follows the AMOC anomaly in both western and eastern SPNA. The minima coincide in phase at *φ* = *π*/2, but there is a slight phase lag between the respective maxima close to *φ* = 3*π*/2. The minimum (maximum) SST anomaly is −0.28 K (+0.23 K) in the western SPNA and −0.37 K (+0.35 K) in the eastern SPNA. The big contrast occurs with OHC (red), which is shifted by a quarter cycle in the western SPNA and a little less (~1/8 of a cycle) in the eastern SPNA, consistent with the lagged correlations presented in Table 3.

Going further, we can form composites of all the quantities in Eqs. (9)–(12). Figure 12c shows composites of the rate of change of heat content $\u2061(\u2202\Theta FD*/\u2202t)$ in the western SPNA (light blue line) together with the surface heat flux ($QNET*$, red) and advection ($RFD*$, dark blue) terms, with the AMOC anomaly (black) at 26°N superimposed for reference. The rate of change of heat content is negative (i.e., heat content is falling) in phases 1 and 2, rises steeply to positive values in phase 3, and declines more slowly in phase 4. Advection closely matches the net tendency (dark and light blue curves) during phases 1–3, but is significantly higher in phase 4. The surface flux term is positive in phases 1–3, weakly opposing the advection term, and rises slightly. In the middle of phase 3, as the heat content peaks, the surface flux term declines steeply, transitioning to negative values in phase 4. Overall it can be seen that the net tendency is largely driven by advection, but in phase 4 there is strong damping by surface fluxes. A similar conclusion can be drawn for the ocean heat content in the eastern SPNA (Fig. 11b). In the western SPNA, the advection term is very clearly related to the AMOC anomaly with a lag of approximately Δ*φ* = *π*/4 (Fig. 12c) whereas in the eastern SPNA the advection covaries with the AMOC anomaly (Fig. 12d). We thus conclude that the AMOC is main driver of large-amplitude decadal variations in OHC.

The SST tendency behaves in a broadly similar way [Figs. 12e and 12f, where the net tendency, ∂*ξ**/∂*t*, is in light blue, the surface flux–related term, $\u2061(1\u2212\lambda )QNET*/\rho 0CPh\xaf$, is red and the advection-entrainment term, $\mathbb{R}ML*/\rho 0CPh\xaf$, is dark blue; the AMOC anomaly at 26°N (black) is again overplotted for reference] but there are some subtle differences. In the western SPNA we see a larger contrast compared to OHC (Fig. 12e). The net tendency (light blue) peaks earlier than the net OHC tendency in the same region (Fig. 12c) and because both quantities have essentially the same surface flux forcing (red) it must be the advection-entrainment term in the mixed layer which is responsible (dark blue). Of interest is the fact that both the net tendency and the advection-entrainment term lead the AMOC and the surface flux term leads the net tendency term. Thus surface fluxes seem to exert some control on the SST in the western SPNA. In the eastern SPNA, the SST and OHC tendency behave very similarly (Fig. 12f) and in particular in both cases, the surface flux term is of opposite sign to the SST suggesting the surface flux term is chiefly having a damping effect. The results strongly suggest that advection is the dominant process controlling the evolution of the OHC in the both the western and the eastern SPNA and additionally, advection-entrainment is the process controlling SST in the eastern SPNA. In the western SPNA, there is a disconnect between the full depth advection and the advection-entrainment in the mixed layer, resulting in an SST peak substantially before the heat content peak. In the eastern SPNA, by contrast the full depth and mixed layer tendencies work in tandem and there is little difference in the timing of the peaks. This explains why there is a lag between OHC peaks in the western and eastern SPNA, but no lag between the SST peaks.

The OHC advection term follows the AMOC according to expectations but surface fluxes release the extra heat input to the atmosphere when the AMOC is rising but the AMOC anomaly is still negative (i.e., in phase 2 of the composite event). It is only when the AMOC anomaly becomes positive that the heat content begins to rise. When the AMOC is falling in phase 4, advection falls too, but OHC increases because the opposing contribution of surface fluxes falls faster. A period of decreasing OHC follows when surface fluxes begin to rise at about the time that the AMOC is halfway between its peak and zero (particularly marked in the western SPNA).

As already noted, the net SST tendency in Fig. 12e leads the AMOC anomaly at 26°N in the western SPNA. Since the advection-entrainment term also lags the SST tendency, but the surface flux term leads all three, this suggests that surface fluxes in the western SPNA are at least partly responsible for the large AMOC variations seen in the model. But the surface fluxes are partially set by the AMOC through its (eventual) control of the SST (via subtropics and the eastern SPNA) emphasizing the coupled nature of the AMOC variability.

### h. Drivers of OHC and SST variations

In this section we summarize the driving terms which characterize the AMOC cycle (Fig. 13). Recalling from section 4e that the term with the larger absolute magnitude [either surface flux related or advection (-entrainment) related] drives the sign of the OHC tendency. If the other term is of opposite sign then it will act as a brake; if it is of the same sign then the two terms act in concert. In Fig. 13a the net OHC tendency for the composite AMOC event is shown in black. We then divide the cycle into regimes depending on which term is dominant (i.e., either |*R*_{FD}| > |*Q*_{NET}| or more rarely |*Q*_{NET}| > |*R*_{FD}|). For each regime the corresponding terms are averaged over the duration of the regime and a constant value plotted in order to quantitatively depict the interplay between the forcing terms during each regime. These regimes do not in general line up with the AMOC phases (p1–p4); for example, midway between phase 1 at *φ* ~ *π*/4 to partway through phase 3 (*φ* ~ 1.1*π*) advection (blue) is the driving term with an average value of approximately −3.0 W m^{−2} and it is opposed by surface fluxes (red) with an average value of approximately +1.0 W m^{−2}. In the subsequent regime, for a brief period surface fluxes dominate as the advection term transitions from negative to positive values as does the net tendency itself. From here to the peak net tendency (*φ* ~ 1.2*π* to *φ* ~ 1.4*π*) the two terms act in concert after which surface fluxes transition to negative values. Advection remains the dominant term in this regime, but receives substantial opposition from surface fluxes.

The situation in the eastern SPNA (Fig. 13b) is similar, but the cycle is shifted to earlier times with respect to the west. As with the west, there is a shift from negative to positive forcing by advection halfway along the period when the net tendency increases (*φ* ~ 0.9*π*) and a shift from positive to negative surface flux forcing close to the time of peak net tendency (*φ* ~ 1.3*π*). In addition, there is an extended regime where surface fluxes are the dominant term (*φ* ~ 1.7*π* to *φ* ~ 2.0*π*) which is not seen in the west. Despite this, advection is clearly the dominant term over most of the cycle for both eastern (66% of the time) and western SPNA (88% of the time). The equivalent plots for the SST are shown in Figs. 13c and 13d. These are quite similar to the OHC plots, especially for the eastern SPNA, but it is noteworthy that surface fluxes play a more important role especially in the west, where there is a long period from *φ* ~ 0.6*π* to *φ* ~ 1.4*π* during which surface fluxes dominate, albeit sometimes narrowly. In both east and west, surface fluxes dominate from *φ* ~ 1.7*π* to *φ* ~ 2.0*π*. Overall advection dominates only 53% of the time in the western SPNA and 61% of the time in the eastern SPNA. Unlike the composite terms in Figs. 12e and 12f the results shown in Figs. 13c and 13d are not very sensitive to whether or not we use the unadjusted [Eq. (10)] or adjusted tendency terms [Eq. (12)].

## 5. Conclusions

We have developed a novel combined approach to the mixed layer and full-depth ocean heat budgets and used it to investigate sea surface temperature (SST) and ocean heat content (OHC) variability on decadal to multidecadal time scales in the subpolar North Atlantic (SPNA), the main center of action of the Atlantic multidecadal variability (AMV). Our analysis has employed a state-of-the-art coupled climate model, HadGEM3-GC2, in which the simulated AMV index and spatial pattern is very similar to observed estimates. The new elements of the approach are development of an equation for evolution of anomalous SST and a parameterization of the diffusive heat flux at the base of the mixed layer.

The results of our analysis show that both OHC and SST tendencies are the result of a competition between two terms representing the effects of surface fluxes and advection for OHC (advection-entrainment for SST). These terms have different forms in the OHC and SST equations, because additional terms related to entrainment appear in the SST equation but not in the OHC equation. Hence, the relationship between OHC and SST becomes an investigation into how and why the surface fluxes and advection-related terms differ between the OHC and SST equations.

The main conclusions are listed below:

Anomalies in the OHC tendency propagate around the SPNA on decadal time scales with a clear relationship to the phase of the AMOC.

In the SPNA, AMOC anomalies lead SST anomalies, which in turn lead OHC anomalies. This result does not depend on the depth used for calculation of OHC and is common to both eastern and western SPNA.

OHC variations in the SPNA on decadal time scales are largely dominated by AMOC variability because it controls variability of advection which is shown to be the dominant term in the OHC budget. Surface heat fluxes modulate the OHC variability, particularly as OHC peaks and declines. Surface heat flux plays a larger role in SST variability.

The advection term covaries with the AMOC in the eastern SPNA but lags the AMOC in the western SPNA, leading to the anticlockwise propagation of OHC anomalies around the SPNA.

The lag between OHC and SST is traced to differences between the advection term for OHC and the advection-entrainment term for SST. The latter leads the former particularly in the western SPNA.

In the western SPNA, surface fluxes and SST appear to precede and cause AMOC changes, whereas in the east AMOC changes cause the changes in SST and surface fluxes.

The main implication of our study is that deep OHC changes are not associated with immediate changes in SST in HadGEM3-GC2; indeed, changes in SST precede OHC deep changes. There is also a very clear difference in the dominant process between the eastern and western SPNA. In the former region, advection is dominant, whereas in the latter surface fluxes dominate. While our study confirms the important role of the AMOC in the decadal variability of the North Atlantic SST, this role cannot be simplified as an increasing AMOC leading to increasing heat content leading to increasing SST, which is a common assumption underlying numerous studies of contemporary and paleo variability of the North Atlantic (e.g., Chen and Tung 2018). For example, in this study using HadGEM3-GC2 the SPNA OHC rarely immediately increases as AMOC increases (phase 2 in Fig. 12), because the advection term must first switch sign from negative to positive (Figs. 13a,b). On the other hand the SST can and does begin rising quite soon after the AMOC starts increasing, because the surface flux term is already driving an increasing SST at this time and reduced opposition to this term from advection reinforces this trend.

In the western SPNA in particular it seems that surface fluxes drive both the subsequent evolution of the advection-entrainment term, and ultimately the AMOC. The detailed mechanism by which surface fluxes can influence the advection still need to be determined, but may be related to the projection of short (seasonal to interannual) time scale correlations between MLD and temperature onto the decadal time scale [see appendix A, Eq. (A9)].

The diagnostic framework developed here is eminently suitable for use with observations and multimodel ensembles. For observations, however, great care must be taken in analysis of errors as rates of change of both OHC and SST consist of a fine balance (i.e., a small residual) between large competing terms of opposite sign. In addition, decadal-scale observational analysis would require high-quality mixed layer depth observations that are still not available globally. Finally, we note that the new framework can be usefully applied to the CMIP6 model ensemble in order to establish the robustness of the results, and to reveal individual model deficiencies that could help usefully constrain climate change projections.

## Acknowledgments

This research was funded by the NERC ACSIS Programme (Grant NE/N018044/1), the NERC funded RAPID AMOC programme at 26°N, the DYNAMOC project (NE/M005097/1), the SMURPHS project (NE/N005686/1 and NE/N005767/1), and funding from the European Union Horizon 2020 research and innovation programme BLUE-ACTION (Grant 727852). The coupled climate model data are available from the Met Office Managed Archive Storage System (MASS). Analysis codes are available from the corresponding author on request.

### APPENDIX A

#### Derivation of SST Anomaly Equation

we first isolate the time derivative terms

then aggregate terms

where

decompose *h* and *T*_{a}, *X*, *T*_{−h}, *Q*_{NET}, and *Q*_{−h} into mean and anomaly components, denoted by an overbar and an asterisk respectively, in Eq. (A3)

take the mean of Eq. (A5)

consolidate terms, parameterize $Q\u2212h*=\lambda QNET*$, and divide by $h\xaf$

where

### APPENDIX B

#### Optimal Value for Diffusive Heat Flux Fraction *λ*

*λ*

As explained in section 4d we obtain an optimal value for *λ* by ensuring that the resulting mixed layer advection entrainment term $\mathbb{R}ML*$ has a maximum correlation with the full depth advection term *R*_{FD}. Figure B1a illustrates this correlation for the western (black) and eastern (red) SPNA for values of *λ* between 0.91 and 1.0. It is remarkable that such a maximum correlation with nonnegligible value exists, ~0.63 for the western and ~0.88 for the eastern SPNA. Corresponding lags are shown in Fig. B1b and indicate that the mixed layer term precedes the full-depth term by three years in the western SPNA and that the two terms are simultaneous in the eastern SPNA. For the purposes of this paper we choose a compromise value of *λ* = 0.99.

## REFERENCES

*Data Analysis Methods in Physical Oceanography*. Pergamon, 634 pp.

*Numerical Recipes: The Art of Scientific Computing*. Cambridge University Press, 818 pp.

*The Encyclopedia of Life Support Systems (EOLSS)*Tropical Meteorology, http://iprc.soest.hawaii.edu/users/xie/o-a.pdf.

## Footnotes

© 2019 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (http://www.ametsoc.org/PUBSReuseLicenses).