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

HAL
open science

# The Active and Passive Roles of the Ocean in Generating Basin-Scale Heat Content Variability

Dafydd Stephenson, Florian Sévellec

► To cite this version:

Dafydd Stephenson, Florian Sévellec. The Active and Passive Roles of the Ocean in Generating Basin-Scale Heat Content Variability. Geophysical Research Letters, 2021, 48 (19), 10.1029/2020GL091874. hal-03445784

HAL Id: hal-03445784

https://hal.science/hal-03445784

Submitted on 16 Jun 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Geophysical Research Letters®

## RESEARCH LETTER

10.1029/2020GL091874

**Key Points:**

* Surface ocean heat content variability is controlled by passive processes in all basins

* Full-depth ocean heat content variability is controlled by active ocean feedbacks in most basins

* Active ocean feedbacks act as a source of decadal predictability in the North Atlantic

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

D. Stephenson, D.Stephenson@noc.soton.ac.uk

**Citation:**

Stephenson, D., & Sévellec, F. (2021). The active and passive roles of the ocean in generating basin-scale heat content variability. *Geophysical Research Letters*, *48*, e2020GL091874. https://doi.org/10.1029/2020GL091874

Received 1 DEC 2020

Accepted 3 SEP 2021

© 2021. The Authors.
This is an open access article under
the terms of the Creative Commons
Attribution License, which permits use,
distribution and reproduction in any
medium, provided the original work is
properly cited.

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

## The Active and Passive Roles of the Ocean in Generating Basin-Scale Heat Content Variability

Dafydd Stephenson¹ ID and Florian Sévellec¹,²

¹Ocean and Earth Science, University of Southampton, Southampton, UK, ²Laboratoire d'Océanographie Physique et Spatiale, Université de Brest CNRS IRD Ifremer, Brest, France

**Abstract** The role of ocean circulation in transforming surface forcing into interannual-to-multidecadal oceanic variability is an area of ongoing debate. Here, a novel method, establishing exact causal links, is used to quantitatively determine the role of ocean active and passive processes in transforming stochastic surface forcing into heat content variability. To this end, we use a global ocean model in which the dynamical response to forcing can be switched on (fully active) or off (purely passive) and consider the resulting effect on heat content variance. While passive ocean processes mainly control the surface variance (over 92%) in all basins, most regions show the importance of active processes at depth. This role is particularly important for full-depth North Atlantic heat content, which we investigate further, highlighting signatures of the meridional overturning circulation in delaying the variance growth.

**Plain Language Summary** The ocean's role in climate is fundamental due to its ability to absorb significant amounts of heat relative to the other components of the Earth system. However, changes in heat can modify the ocean currents which transport that heat. The importance of this feedback effect remains uncertain, and so our study aims to determine how important this process is. We achieve this by alternately switching on and off the ability of simulated ocean currents to respond to changes in heat and salt driven by the atmosphere in a state-of-the art numerical simulation of the ocean. We then compare how variable the heat content of the ocean is in both "on" and "off" cases. We show that ocean circulation changes are unimportant near the surface, but in most regions they play a key role at depth. We look in detail at the North Atlantic, the region where circulation changes have the most important effect.

## 1. Introduction

It is well documented that the oceanic heat reservoir has a crucial role in climate; the ocean has absorbed over 90% of the excess energy associated with anthropogenic warming (Trenberth et al., 2014), for instance. However, this single number obscures the spatiotemporal heterogeneity of ocean heat content change, which is punctuated by hiatuses and surges (e.g., Meehl et al., 2011), geographically differential warming (e.g., Drijfhout et al., 2012), and varying impacts at different depths (e.g., Gleckler et al., 2016). The mechanisms underlying these variations remain challenging to disentangle due to the complexity of the climate system. This is particularly relevant on interannual-to-multidecadal timescales, where natural variability and external forcing have comparable amplitude (Meehl et al., 2009). Understanding these variations is thus crucial for modeling and predicting them.

The simplest explanation of heat content anomalies in the ocean is that they originate in the atmosphere, either via external forcing or natural, internal fluctuations, are fluxed into the mixed layer, and then passively circulated around the ocean interior along its preferred ventilation pathways. In this paradigm, the anomalous heat can be considered density compensated in that the ocean circulation does not change (e.g., Mauritzen et al., 2012). This approximation is often applied when modeling the long-term response to anthropogenic forcing (e.g., Marshall et al., 2015; Newsom et al., 2020; Zanna et al., 2019), with anomalous heat fluxes represented by a passive tracer. However, investigations of the validity of this approximation for heat uptake typically flag the North Atlantic as a region to which it is particularly ill-suited (Banks & Gregory, 2006; Garuba & Klinger, 2016, 2018; Xie & Vallis, 2012), due to the Atlantic Meridional Overturning Circulation (AMOC) and its link with heat storage in models (Kostov et al., 2014). The involvement of the AMOC in interannual-to-multidecadal ocean temperature variations remains a contentious issue, however.
Recent studies have argued that the predominant patterns of Atlantic Multidecadal Variability (AMV) in climate simulations featuring realistic ocean general circulation models (OGCMs) can be recreated by coupling a realistic atmosphere to a time-invariant “slab” ocean (Cane et al., 2017; Clement et al., 2015, 2016), suggesting that these patterns are purely passive. In this slab ocean case, commonalities with fully active ocean simulations can only be established statistically. On the other hand, the aforementioned passive tracer approach (by propagating a passive “temperature” tracer initially coincident with the active temperature field in a single simulation and considering their divergence) provides a more clear separation of the response into its passive and dynamical components. Nevertheless, statistical slab-OGCM comparisons remain the de facto standard for determining the role of the ocean in near-term regional low-frequency variability (Dommenget, 2010; Dommenget & Latif, 2002; Delworth et al., 2017; Wang & Dommenget, 2016; Zhang, 2017).

In this study, we present an alternative approach to the question of regional heat content variability, using an adjoint model. Unlike a conventional model, which integrates anomalies forward in time, an adjoint model describes the sensitivity of a quantity of interest (here heat content) to past changes (here stochastic atmospheric forcing), establishing causes, rather than effects (Errico, 1997). As such, the method describes all possible perturbations leading to a particular ocean response, otherwise requiring a theoretically “infinite” ensemble of perturbed simulations. This has been leveraged to attribute the sources of temporal ocean variability in response to historical atmospheric forcing (Kostov et al., 2021; Pillar et al., 2016; Smith & Heimbach, 2019) and to establish the growth of oceanic variance in response to a representation of stochastic atmospheric forcing (Sévellec et al., 2018).

We use this approach to isolate the role of the ocean in modeled heat content variability, by projecting a diagnosed stochastic representation of atmospheric buoyancy and momentum fluxes onto passive and active surface adjoint sensitivity fields. In the passive case, buoyancy anomalies cannot change the circulation.

## 2. Method and Diagnostics

To characterize low-frequency ocean variability, Hasselmann (1976) and Frankignoul and Hasselmann (1977) developed an idealized, single-variable stochastic model of ocean surface temperature in response to random heat fluxes. These atmospheric fluxes can be seen as a continuous stream of small disturbances to ocean surface temperature, which accumulate and are slowly “forgotten”. This can be represented as

$$u(t) = \int_{0}^{t} e^{-\lambda(t-\tau)} L dW(\tau), \qquad (1)$$

where $u(t)$ is the ocean temperature anomaly at time $t$ (using $u(0) = 0$ without loss of generality), $\lambda$ is the inverse damping timescale representing the ocean dynamics, $W(\tau)$ is a standard-normal Wiener process with time index $\tau$ (such that for time $t$, the impact of $W$ originating at time $\tau$ is small if $t \gg \tau$), and $L^2$ describes the intensity of the stochastic fluxes (variance of their temperature impact per unit time).

A large body of work has generalized Equation 1 to multivariate systems featuring multiple interacting variables and locations, ranging from box models (Cessi, 1994; Chang et al., 2004; Griffies & Tziperman, 1995; Grötzner et al., 1999; Zanna & Tziperman, 2008) to GCMs (Hawkins & Sutton, 2009; Sévellec et al., 2018; Tziperman et al., 2008). In the higher order case, different variables and locations may all be represented by a single anomaly state vector: $|u\rangle$. Here we utilize the bra-ket notation of Dirac (1939), whereby vectors are written as “kets,” for example, $|b\rangle$, and co-vectors as “bras,” for example, $\langle al, a|$, such that bra-ket pairs become either scalar through the Euclidean inner product, for example $\langle al|b\rangle = c$ or matrices through the outer product, for example $|b\rangle\langle al| = C$. We may also represent more involved linear processes than simple exponential decay using (e.g., Farrell & Ioannou, 1996a, 1996b) the propagator, $\Psi$, of the ocean model which linearly maps anomalies $|u(0)\rangle$ to their later states $|u(t)\rangle = \Psi(t, 0)|u(0)\rangle$, with $\Psi(t, 0) = e^{-\lambda t}$ in the univariate case of Equation 1. This reads:

$$|u(t)\rangle = \int_{0}^{t} \Psi(t, \tau) L dW(\tau), \qquad (2)$$

where $|W\rangle$ is a vector of independent standard-normal Wiener processes and $\Sigma = LL^\dagger$ is a covariance matrix (describing the stochastic flux intensity similarly to $L^2$ in the univariate case of Equation 1 but also accounting for coherence between locations and variables). As before, $|u(0)\rangle$ is assumed zero-valued.
From this formula, one can obtain the outcome of a quantity of interest $\langle \mathbf{F}|\mathbf{u} \rangle$, such as heat content, in a fully active ($\Psi_A$) or purely passive ($\Psi_P$) ocean model. While the heat content variation in a fully active model is a classical problem of modern ocean physics, it is important to explicitly describe the routes by which its purely passive component can exhibit heat content anomalies. The first is the fluxing of heat content anomalies from the atmosphere which then propagate through the ocean by mean advection and diffusion. The second is the introduction of circulation anomalies by the wind. Although buoyancy anomalies cannot modify the circulation in the purely passive case, momentum fluxes may still create an anomalous circulation. This can then create heat content variations by redistributing the underlying mean ocean temperature field.

Given the quantity of interest $\langle \mathbf{F}|\mathbf{u} \rangle$, one can also compute its variance from Equation 2 as described by, for example, Sévellec et al. (2018). The covariance ($\sigma_{AP}$) between outcomes in the two configurations (fully active and purely passive) of the model can similarly be calculated to determine their common components. Using the defining property of the adjoint $\langle \mathbf{a}|\mathbf{X}|\mathbf{b} \rangle = \langle \mathbf{b}|\mathbf{X}^{\dagger}|\mathbf{a} \rangle$ (where $\mathbf{X}$ and $\mathbf{X}^{\dagger}$ are an operator and its adjoint) and following a multi-dimensional generalization of Itô's isometry (e.g., Section 3.6 of Duan & Wang, 2014), the covariance at time $t$ reads:

$$ 
\begin{aligned}
\sigma_{AP}(t) &= \operatorname{Cov}(\langle \mathbf{F}|\mathbf{u}_A(t) \rangle, \langle \mathbf{F}|\mathbf{u}_P(t) \rangle) \\
&= E\left[\langle \mathbf{F}|\int_0^t \Psi_A(t, \zeta) \mathbf{L} d\langle \mathbf{W}(\zeta) \rangle \langle \mathbf{F}|\int_0^t \Psi_P(t, \gamma) \mathbf{L} d\langle \mathbf{W}(\gamma) \rangle \right] \\
&= \int_0^t \langle \mathbf{F}|\Psi_A(t, \tau) \Sigma \Psi_P^{\dagger}(\tau, t) |\mathbf{F} \rangle d\tau,
\end{aligned}
\quad (3) $$

where $\sigma_{AP}$ is the covariance between fully active and purely passive version of the model denoted by $\Psi_A$ and $\Psi_P$, respectively, and $E[\cdot]$ is the expectation of a stochastic Itô process. This leads to expressions for the variance of the fully active and purely passive component at time $t$:

$$ 
\begin{aligned}
\sigma_A^2(t) &= \operatorname{Var}(\langle \mathbf{F}|\mathbf{u}_A(t) \rangle) = \int_0^t \langle \mathbf{F}|\Psi_A(t, \tau) \Sigma \Psi_A^{\dagger}(\tau, t) |\mathbf{F} \rangle d\tau, \\
\sigma_P^2(t) &= \operatorname{Var}(\langle \mathbf{F}|\mathbf{u}_P(t) \rangle) = \int_0^t \langle \mathbf{F}|\Psi_P(t, \tau) \Sigma \Psi_P^{\dagger}(\tau, t) |\mathbf{F} \rangle d\tau.
\end{aligned}
\quad (4) $$

These equations describe the level of variance of the ocean heat content obtained after a time $t$ in response to stochastic forcing starting from rest in the fully active ($\sigma_A^2$) and purely passive ($\sigma_P^2$) cases. Of note is that the disappearance of the stochastic terms $d\langle \mathbf{W} \rangle$ means that it is not necessary to propagate any particular realization of random noise to obtain these properties—the variance of an “infinite” ensemble of such realizations can be determined using just the properties of the noise and a single run of the linear model. These will asymptotically tend towards their associated climatological heat content variance. The covariance describes how much of this variance is common to both, and can be normalised to give a Level of Agreement (LoA) between the purely passive and fully active cases, which we define as $\text{LoA}(t) = \frac{\sigma_{AP}(t)}{\sigma_A(t)\sigma_P(t)}$. If the LoA is unity at a given time, it is taken that anomalous heat content variation in the fully active ocean has been entirely controlled by purely passive processes.

These diagnostics have three requirements. Firstly, a linearized OGCM is needed to provide the propagator $\Psi_A$ and its adjoint $\Psi_A^{\dagger}$. Secondly, this propagator requires an isolated purely passive component $\Psi_P$ and its adjoint $\Psi_P^{\dagger}$. The model, its adjoint, and the purely passive configuration are described in Section 3.1. Lastly, we require a stochastic representation $\Sigma$ of surface fluxes. We diagnose this from a coupled climate model (also described in Section 3.1). In particular, we assume that buoyancy and momentum flux anomalies from the coupled simulation climatology follow a band-limited (therefore finite power), spatially covarying Gaussian white noise. At each location, the power spectral density (PSD) of the flux anomalies is therefore assumed constant up to a few days, and zero at higher frequency. The cutoff is determined by the e-folding decorrelation timescales of the fluxes (Figure 1, contours). We also have an implicit low-frequency limit imposed by the 20-year length of the coupled simulation. The elements of $\Sigma$ are then given by the (effectively constant) PSD averaged over this band.

It is important to remark on linearity and independence, which allow for further decomposition of the above diagnostics. As the model propagators are linear, we can consider the fully active model $\Psi_A$ to be the sum of the purely passive model $\Psi_P$ and a dynamical-only component $\Psi_D$ which encompasses just the feedback terms. Furthermore, the propagation of multiple quantities is equal to the propagation of their sum by linearity: $\Psi^{\dagger}(\langle \mathbf{F}_1 \rangle + \langle \mathbf{F}_2 \rangle) = \Psi^{\dagger} \langle \mathbf{F}_1 \rangle + \Psi^{\dagger} \langle \mathbf{F}_2 \rangle$.
<figure><img src="image_4.png" /></figure>

**Figure 1.** Local (co)variance (shading) and decorrelation time ($\lambda^{-1}$, contours) of surface fluxes in the coupled model. (a) Variance in rate of temperature change due to heat flux (HF). (b) Variance in rate of salinity change due to freshwater flux (FWF). (c) Covariance between rate of temperature and salinity change. (d and e) Variance in rate of zonal and meridional velocity change due to zonal and meridional momentum fluxes (ZMF and MMF), respectively. (f) Covariance between rate of zonal and meridional velocity change. Thick dashed, solid, and dotted black contours indicate decorrelation time ($\lambda^{-1}$) of 1, 2, and 3 days, respectively. Thin gray contours are intermediate values, separated by 0.5 days.

We additionally take surface buoyancy fluxes (described by $\Sigma_B$) and momentum fluxes (described by $\Sigma_M$) to be independent, and so the response to each can be determined separately, with $\Sigma = \Sigma_B + \Sigma_M$. We emphasize that the covariance between the buoyancy components (heat and freshwater fluxes) and between the momentum components (zonal and meridional fluxes) remain fully acknowledged. Using this to calculate $\sigma_{AP}$ separately in response to buoyancy only and momentum only allows the LoA to be partitioned accordingly, by modifying its numerator while retaining the denominator. Finally, although the diagnostics of $\sigma_A^2$, $\sigma_P^2$ and $\sigma_{AP}$ are scalar values, they can be computed elementwise without summation, such that the contribution of each variable at each location to the total can be isolated. Similarly, the time integral can be decomposed to obtain the contribution of any time interval. This permits us to see the surface distribution and timing of sources leading to the resulting (scalar) heat content variance.

# 3. Application to an OGCM

## 3.1. Model Description

Our stochastic representation $\Sigma$ in Equation 3 is constructed from thermal, haline, and zonal and meridional momentum fluxes diagnosed from a coupled climate model (Figure 1). Specifically, a 20 year simulation using the IPSL-CM5A-LR coupled model was run in its CMIP5 pre-industrial control configuration (Dufresne et al., 2013) with daily average output. Atmospheric variations at higher or lower frequencies than those captured by this simulation were taken not to influence heat content variance on the spatiotemporal scales considered here. The model was chosen as its ocean component is NEMO (v3.2) with its ORCA2 global configuration (2° nominal resolution with 31 vertical levels), similarly to our linearized ocean model (described below). The atmospheric component is the LMDZ5a model, with 3.75° × 1.9° horizontal resolution and 39 vertical levels (Hourdin et al., 2013).

The linear ocean model which provides the propagator $\Psi_A$ of Equation 3 is NEMOTAM (Vidard et al., 2015), which is derived from NEMO v3.4 (Madec, 2012) and is used in its ORCA2-LIM configuration. The model configuration is similar to that detailed in Stephenson et al. (2020), which also discusses the implementation of the purely passive configuration in detail. The nonlinear model, which provides the simulation about which NEMOTAM is linearized, is forced by a single representative year (CORE normal year forcing, phase 1, v2.0; Large & Yeager, 2004, 2009). This allows the ocean response described by $\Psi$ in Equation 2 to
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Level of Agreement (LoA) between purely passive and fully active simulations in generating the accumulated final (60-year) heat content variance due to buoyancy (red) and momentum (green) surface stochastic fluxes, determined by calculating $\sigma_{AP}$ in response to each. LoA is shown for the three cases (surface layer—corresponding to SST, upper 1,500 m, and full-depth heat content). Largest bar plot shows the case for the total global ocean heat content variance, smaller inner plots show regional values. Thinner dashed black lines signify a LoA of 50%. Black solid lines on the map mark the boundaries of the regions in our definitions.</figcaption></figure>

be that of a “typical” year, cleanly separated from interannual surface variations, which were applied separately through the matrix $\Sigma$ in Equation 3.

## 3.2. Results

We now apply the derivations of Section 2 to attribute the generation of variance in heat content (($F$ | in Equations 3 and 4) in the fully active simulation to its different sources. We evaluate heat content over three depth ranges (10 m, 1,500 m, and full-depth), which effectively correspond to sea surface temperature, heat content in the upper ocean, and the total heat content, respectively. We also consider both the global ocean and a seven-region partition of it (Figure 2, black lines). These regions are the Arctic Ocean (>70°N), the North (≥35°N, ≤70°N) and intertropical (≥-35°N, <35°N) Atlantic and Pacific, the Indian Ocean (≥-35°N), and the Southern Ocean (<-35°N).

Our analysis reveals that the LoA between purely passive and fully active heat content variance after 60 years varies significantly depending on the depth extent and geographical region (Figure 2, bars). The LoA is extremely high for sea surface temperature variance in all regions. This ranges from 92% in the intertropical Pacific to over 99% in the Southern Ocean, with a majority stimulated by buoyancy forcing. This implies that the purely passive uptake of heat controls temperature variability at the surface. There is a dramatic reduction in agreement when heat content is computed over a thicker layer. For the upper-1,500 m heat content, variance common to both the purely passive and fully active simulations accounts for as little as 31% in the case of the Indian Ocean, and just over half (52%) globally. The nature of stimulation of the purely passive component also changes over this depth range, shifting to a primarily wind-driven regime for all regions except the Arctic Ocean.

When heat content is defined over the full depth, it generally follows similar patterns to upper-1,500 m heat content, with notable exceptions in the North Atlantic and Arctic oceans. For those basins, another dramatic reduction in correspondence between the purely passive and fully active simulations occurs, with the LoA reducing to 27% and 25%, respectively. More subtle reductions can be seen elsewhere, and only in the North Pacific and Southern Ocean does the purely passive component still dominate the fully active simulation at full depth. It is worth noting the substantial impact (>50%) of purely passive wind effects in these regions.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Evolution of full-depth heat content variance in response to stochastic surface forcing in the purely passive (dotted lines) and fully active (solid lines) simulations. The difference between these lines is linked to the dynamical-only component, which may act destructively (passive > active) or constructively (active > passive). Thinner lines show separately the buoyancy-forced (red) and wind-driven (green) components. Stars mark the point at which 50% of the final (60-year) variance is reached.</figcaption></figure>

While the LoA provides a useful quantification of the ultimate role of the purely passive component of the ocean, it does not describe in detail the differences between the purely passive and fully active simulations (e.g., the timing of the variance growth or its source location). To tackle this question, we consider the time-evolving variance growth for each, along with its components (Figure 3). We focus on the full-depth case, where differences between these components are greatest. Similar decompositions have been considered for surface (Figure S1) and upper ocean (Figure S2) cases, and exhibit similar (but less significant) behavior.

The temporal evolution of the variability reveals that the purely passive and fully active simulations differ in both magnitude and timing. As discussed in Section 2, linearity permits the decomposition of the fully active model into the sum of the purely passive component and a remaining dynamical-only component. The difference between evolving variance in the fully active model and in the purely passive model (Figure 3, solid and dotted lines, respectively) can thus be attributed to internal ocean feedbacks belonging to this dynamical-only component, which are not always constructive. Indeed, the variance in the fully active simulation is often weaker than that of its purely passive counterpart. This suggests that certain behavior is possible only in the purely passive case, and is canceled out by the dynamical-only term in the fully active simulation. This is particularly visible for heat content variance in the Indian and Southern Oceans (dominated by wind stress). There, after two decades, most of the variance growth of the purely passive
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Surface sources of (a, c, e, g) buoyancy- and (b, d, f, h) wind-stimulated full-depth heat content variance in the North Atlantic, integrated over years 0–10 (upper panels) and 10–60 (lower panels) in the (a, b, e, f) fully active and (b, d, g, h) purely passive simulations, as well as in the dynamical-only diagnosed component (contours in c, d, g, h). The global integrals of the fully active and purely passive fields for the upper panels produce the values shown in Figure 3 at 10 years. The further addition of the global integral of the fields from the lower panels provides the values shown in Figure 3 at 60 years. The dynamical-only component is defined as the difference between the fully active and purely passive simulations. Solid and dashed contours indicate positive and negative values, respectively, with contour intervals of 0.05 (EJ)² km⁻² for buoyancy, and of 0.2 (EJ)² km⁻² for momentum.</figcaption></figure>

component stimulated by wind stress is canceled by the dynamical-only component. An example of similar behavior is provided by Cronin and Tozuka (2016), who demonstrate that Ekman transport is determined not only by wind stress and latitude (as in the classical analysis of Ekman, 1905), but also local geostrophic shear. In this perspective, Ekman transport has both a purely passive and dynamical-only component, which can act against each other.

As a measure of the rate at which the climatological variance is approached, we consider the time taken for the full-depth variance in each simulation to reach half of its final (60 year) value (Figure 3, stars). Following on from the previous discussion, the dynamical-only momentum component in the Indian Ocean acts to reduce the final variance achieved in the active case, thereby reaching half of this value in only 7 years, as opposed to 19 years in the purely passive simulation. At the opposite extreme, for the Arctic and North Atlantic, the dynamical-only contribution slows the variance evolution substantially. This effect is ultimately dominated by the buoyancy component, and is made up of similar contributions from temperature and salinity (not shown). Resultantly, in the North Atlantic, half of the final value is reached in only 3 years in the purely passive simulation, compared with 21 years in the fully active case.

The source of this continued growth in the active North Atlantic, even after the purely passive component appears to have saturated, corresponds to a change in the nature of the ocean response to buoyancy stimulation in the fully active simulation after 10 years. To determine the origin of this, we consider separately the surface distribution of the variance accumulated during the first 10 years (Figures 4a–4d) and from 10 to 60 years (Figures 4e–4h). This is determined from the elementwise computation of the variance, prior to summation, as outlined in Section 2.

In the first decade, the passive and active simulations maintain a high LoA (above 75%) and their spatial patterns are similar. Focusing on buoyancy forcing, the relatively focused region reflects the model's deep water formation site, as described in the passive tracer study of Stephenson et al. (2020). The difference between the fully active and purely passive distributions is the dynamical-only contribution (Figure 4c, contours). This corresponds to a large-scale dipole. The negative peak of the dipole overlies the positive contribution
by the purely passive component, having a slight compensating effect (Figure 4c). On decadal timescales, positive contributions to variance growth in both the purely passive and dynamical-only components coincide in location, and so the two components act constructively (shading and contours in Figure 4g).

The primary difference between stimulation by wind for the fully active and purely passive components in the first decade is the intensity of the induced variance (Figures 4b and 4d). Both components are dominated by Ekman transport across a zonal band defining the region's boundary (35°N), with an alternating pattern at the grid scale stimulating a mechanism apparently unique to the passive ocean. The addition of the negative dynamical-only component reduces the intensity of this pattern, however. Also notable in the fully active case is a seemingly persistent (Figures 4b and 4f) stimulation of variance at the subtropical-subpolar inter-gyre interface, as well as stimulation (both positive and negative) in coastal regions of the eastern North Atlantic and Greenland Sea.

## 4. Discussion and Conclusions

We have considered the stimulation of variance in ocean heat content by surface atmospheric noise. We evaluated heat content over a range of different regions and depths in a linearized global ocean model, comparing purely passive and fully active realizations of the ocean model. In the purely passive framework, temperature anomalies either arise due to random surface heat fluxes (and can be passively transported by the mean flow), or due to random surface momentum fluxes (which redistribute existing heat). However, these resulting temperature anomalies are unable to modify the ocean circulation.

In contrast to the established techniques of using a passive tracer (e.g., Banks & Gregory, 2006; Garuba & Klinger, 2016, 2018; Marshall et al., 2015; Xie & Vallis, 2012) or a slab ocean model (e.g., Clement et al., 2015; Dommenget, 2010; Dommenget & Latif, 2002; Wang & Dommenget, 2016) to investigate the role of the ocean, we have utilized a novel adjoint-based approach (Sévellec et al., 2018). The use of an adjoint model has uniquely allowed us to causally attribute heat content variance to different variables, times, and locations at the surface, by projecting onto surface sensitivity fields a representative stochastic representation of atmospheric fluxes diagnosed from a coupled climate model.

Our findings for the surface ocean (i.e., sea surface temperature) are that at least 92% of the variance in the fully active simulation is in agreement with its purely passive component. This is consistent with studies which suggest that oceanic dynamics are not needed to generate surface decadal variability (e.g., Cane et al., 2017; Clement et al., 2015, 2016). While our passive and active variance patterns may in some cases express a high (normalised) LoA, a purely passive model could still misestimate the amplitude of heat content variance, as the purely passive component can be offset or reinforced by the corresponding dynamical-only component in a fully active ocean. In our experiments, this discrepancy occurs over greater depth ranges (Figures 3 and S2), but the amplitudes of surface ocean responses in the passive and active cases are roughly equivalent (Figure S1), as suggested by Clement et al. (2015) in the North Atlantic, for instance.

The dynamical redistribution of existing heat by currents arising from buoyancy anomalies has been shown in past studies to substantially impact heat uptake (e.g., Banks & Gregory, 2006; Xie & Vallis, 2012), particularly in the North Atlantic. However, we have shown that the passive redistribution of the existing heat reservoir by wind anomalies is often more important in the context of heat content variability. This leads to a driving role for the passive component over several regions and depths, in particular the deep Southern and North Pacific oceans. Nevertheless, the full-depth North Atlantic also stands out here as having an important role for ocean feedbacks, with the dynamical-only component acting to slow the growth of heat content variance. We considered the time taken to reach 50% of the variance at the end of the simulation (as a conservative estimate of the time taken to reach half the asymptotic variance), and found that the fully active model takes 7 times as long (21 years) to reach this point as the purely passive simulation (3 years) in this region. This has potential consequences for climate predictability, as the variance growth can also be seen as the accumulation of error following model initialization (Sévellec et al., 2018). The time taken to reach half of the climatological variability is often taken as a measure of the upper limit of predictability, beyond which noise dominates the predictable signal (e.g., Griffies & Bryan, 1997; Grötzner et al., 1999). The reason for this delay in the fully active North Atlantic is a shift in the response to buoyancy forcing. On sub-decadal timescales, the dynamical-only component slows variance growth, before sustaining it on
timescales greater than 10 years, resulting in an “S”-shaped growth curve. In exploring the spatial distribu-
tion of the components of the fully active simulation, we have observed a basin-scale dipole pattern in the
North Atlantic. These patterns echo earlier sensitivity studies of the region in predecessors of our model
(e.g., Sévellec & Fedorov, 2017). These studies relate North Atlantic heat content sensitivity to an ocean-only
mode of variability in which heat content and AMOC anomalies feed back on each other via basin-scale
thermal Rossby wave propagation (Sévellec & Fedorov, 2013) consistently with observations of the AMV
(Frankcombe et al., 2009).

There are a number of considerations which are not accounted for in our approach. First, our conclusions
are likely oversimplified by our use of atmospheric variability sources alone in a linear, laminar model. In
a recent ensemble study at eddy-permitting resolution, Sérazin et al. (2017) suggested that a substantial
portion of ocean heat content variability is intrinsic, generated by chaotic, nonlinear processes within the
ocean. This suggests that we underestimate the role of the dynamical-only component by restricting it to
large-scale, laminar feedbacks. This has been addressed in a separate study (Stephenson & Sévellec, 2021).
In addition, the role of coupling in the stimulation of interdecadal variability is an entire field of research on
its own (e.g., the review of Liu, 2012). Here, our model uses an uncoupled ocean and a stochastic representa-
tion of the atmosphere. This limits the conclusions of our work, in particular for sea surface temperature
(where the surface boundary conditions have more impact). Furthermore, our stochastic representation is
of limited bandwidth, effectively averaging the power spectrum of a two-decade coupled simulation. The
result is a stationary (although globally coherent) white noise representation of daily-to-bidecadal atmos-
pheric variability. We emphasize, however, that these simplifications have allowed us to use an adjoint
ocean model to causally attribute the surface sources of heat content variability exactly, and with limited
computational expense, an approach which offers several unique advantages of its own.

Data Availability Statement

These model outputs can be found at https://doi.org/10.5281/zenodo.4300471 while code used for our ad-
joint simulations and analysis may be found at https://doi.org/10.5281/zenodo.4721438.

Acknowledgments

This research was supported by the
Natural and Environmental Re-
search Council UK (SMURPHS, NE/
N005767/1 and MESO-CLIP, NE/
K005928/1 and the SPITFIRE DTP) and
by the DECLIC and Meso-Var-Clim
projects funded through the French
CNRS/INSU/LEFE program. The au-
thors kindly thank Juliette Mignot and
Victor Estella Perez for providing the
coupled model outputs used to create
the stochastic representation, and
Simon Müller for his assistance with
the adjoint configuration.

References

Banks, H. T., & Gregory, J. M. (2006). Mechanisms of ocean heat uptake in a coupled climate model and the implications for tracer based predictions of ocean heat uptake. *Geophysical Research Letters*, *33*(7), L07608. https://doi.org/10.1029/2005gl025352

Cane, M. A., Clement, A. C., Murphy, L. N., & Bellomo, K. (2017). Low-pass filtering, heat flux, and Atlantic multidecadal variability. *Journal of Climate*, *30*(18), 7529–7553. https://doi.org/10.1175/jcli-d-16-0810.1

Cessi, P. (1994). A simple box model of stochastically forced thermohaline flow. *Journal of Physical Oceanography*, *24*(9), 1911–1920.
https://doi.org/10.1175/1520-0485(1994)024<1911:asbmos>2.0.co;2

Chang, P., Saravanan, R., DelSole, T., & Wang, F. (2004). Predictability of linear coupled systems. Part i: Theoretical analyses. *Journal of Climate*, *17*(7), 1474–1486. https://doi.org/10.1175/1520-0442(2004)017<1474:polcsp>2.0.co;2

Clement, A., Bellomo, K., Murphy, L. N., Cane, M. A., Mauritsen, T., Rädel, G., & Stevens, B. (2015). The Atlantic Multidecadal Oscillation without a role for ocean circulation. *Science*, *350*(6258), 320–324. https://doi.org/10.1126/science.aab3980

Clement, A., Cane, M. A., Murphy, L. N., Bellomo, K., Mauritsen, T., & Stevens, B. (2016). Response to comment on “the atlantic multidecadal oscillation without a role for ocean circulation”. *Science*, *352*(6293), 1527. https://doi.org/10.1126/science.aaf2575

Cronin, M. F., & Tozuka, T. (2016). Steady state ocean response to wind forcing in extratropical frontal regions. *Scientific Reports*, *6*, 28842.
https://doi.org/10.1038/srep28842

Delworth, T. L., Zeng, F., Zhang, L., Zhang, R., Vecchi, G. A., & Yang, X. (2017). The central role of ocean dynamics in connecting the North Atlantic Oscillation to the extratropical component of the Atlantic Multidecadal Oscillation. *Journal of Climate*, *30*(10), 3789–3805.
https://doi.org/10.1175/jcli-d-16-0358.1

Dirac, P. A. M. (1939). A new notation for quantum mechanics. In *Mathematical proceedings of the Cambridge philosophical society* (Vol. 35, pp. 416–418). https://doi.org/10.1017/s0305004100021162

Dommenget, D. (2010). The slab ocean el niño. *Geophysical Research Letters*, *37*(20), L20701. https://doi.org/10.1029/2010gl044888

Dommenget, D., & Latif, M. (2002). Analysis of observed and simulated sst spectra in the midlatitudes. *Climate Dynamics*, *19*(3–4), 277–288. https://doi.org/10.1007/s00382-002-0229-9

Drijfhout, S., Van Oldenborgh, G. J., & Cimatoribus, A. (2012). Is a decline of amoc causing the warming hole above the North Atlantic in observed and modeled warming patterns? *Journal of Climate*, *25*(24), 8373–8379. https://doi.org/10.1175/jcli-d-12-00490.1

Duan, J., & Wang, W. (2014). Stochastic calculus in hilbert space. In *Effective dynamics of stochastic partial differential equations* (pp. 21–47). Elsevier. https://doi.org/10.1016/b978-0-12-800882-9.00003-2

Dufresne, J.-L., Foujols, M.-A., Denvil, S., Caubel, A., Marti, O., Aumont, O., et al. (2013). Climate change projections using the IPSL-CM5 Earth System Model: From CMIP3 to CMIP5. *Climate Dynamics*, *40*(9–10), 2123–2165. https://doi.org/10.1007/s00382-012-1636-1

Ekman, V. W. (1905). On the influence of the Earth's rotation on ocean currents. *Ark. Mat. Astron. Fys.*, 2, 1–52.

Errico, R. M. (1997). What is an adjoint model? *Bulletin of the American Meteorological Society*, *78*(11), 2577–2591. https://doi.org/10.1175/1520-0477(1997)078<2577:wiaam>2.0.co;2
Farrell, B. F., & Ioannou, P. J. (1996a). Generalized stability theory. Part i: Autonomous operators. *Journal of the Atmospheric Sciences*, 53(14), 2025–2040. https://doi.org/10.1175/1520-0469(1996)053<2025:gstpia>2.0.co;2

Farrell, B. F., & Ioannou, P. J. (1996b). Generalized stability theory. Part ii: Nona autonomous operators. *Journal of the Atmospheric Sciences*, 53(14), 2041–2053. https://doi.org/10.1175/1520-0469(1996)053<2041:gstpin>2.0.co;2

Frankcombe, L. M., Dijkstra, H. A., & Von der Heydt, A. (2009). Noise-induced multidecadal variability in the North Atlantic: Excitation of normal modes. *Journal of Physical Oceanography*, 39(1), 220–233. https://doi.org/10.1175/2008jpo3951.1

Frankignoul, C., & Hasselmann, K. (1977). Stochastic climate models, part ii application to sea-surface temperature anomalies and thermocline variability. *Tellus*, 29(4), 289–305. https://doi.org/10.3402/tellusa.v29i4.11362

Garuba, O. A., & Klinger, B. A. (2016). Ocean heat uptake and interbasin transport of the passive and redistributive components of surface heating. *Journal of Climate*, 29(20), 7507–7527. https://doi.org/10.1175/jcli-d-16-0138.1

Garuba, O. A., & Klinger, B. A. (2018). The role of individual surface flux components in the passive and active ocean heat uptake. *Journal of Climate*, 31(15), 6157–6173. https://doi.org/10.1175/jcli-d-17-0452.1

Gleckler, P. J., Durack, P. J., Stouffer, R. J., Johnson, G. C., & Forest, C. E. (2016). Industrial-era global ocean heat uptake doubles in recent decades. *Nature Climate Change*, 6(4), 394–398. https://doi.org/10.1038/nclimate2915

Griffies, S. M., & Bryan, K. (1997). Predictability of North Atlantic multidecadal climate variability. *Science*, 275(5297), 181–184. https://doi.org/10.1126/science.275.5297.181

Griffies, S. M., & Tziperman, E. (1995). A linear thermohaline oscillator driven by stochastic atmospheric forcing. *Journal of Climate*, 8(10), 2440–2453. https://doi.org/10.1175/1520-0442(1995)008<2440:altodb>2.0.co;2

Grötzner, A., Latif, M., Timmermann, A., & Voss, R. (1999). Interannual to decadal predictability in a coupled ocean–atmosphere general circulation model. *Journal of Climate*, 12(8), 2607–2624. https://doi.org/10.1175/1520-0442(1999)012<2607:itdpia>2.0.co;2

Hasselmann, K. (1976). Stochastic climate models part I. Theory. *Tellus*, 28(6), 473–485. https://doi.org/10.3402/tellusa.v28i6.11316

Hawkins, E., & Sutton, R. (2009). The potential to narrow uncertainty in regional climate predictions. *Bulletin of the American Meteorological Society*, 90(8), 1095–1108. https://doi.org/10.1175/2009bams2607.1

Hourdin, F., Foujols, M.-A., Codron, F., Guemas, V., Dufresne, J.-L., Bony, S., et al. (2013). Impact of the LMDZ atmospheric grid configuration on the climate and sensitivity of the IPSL-CM5A coupled model. *Climate Dynamics*, 40(9–10), 2167–2192. https://doi.org/10.1007/s00382-012-1411-3

Kostov, Y., Armour, K. C., & Marshall, J. (2014). Impact of the Atlantic meridional overturning circulation on ocean heat storage and transient climate change. *Geophysical Research Letters*, 41(6), 2108–2116. https://doi.org/10.1002/2013gl058998

Kostov, Y., Johnson, H. L., Marshall, D. P., Heimbach, P., Forget, G., Holliday, N. P., et al. (2021). Distinct sources of interannual subtropical and subpolar Atlantic overturning variability. *Nature Geoscience*, 14, 491–495. https://doi.org/10.1038/s41561-021-00759-4

Large, W. G., & Yeager, S. (2009). The global climatology of an interannually varying air–sea flux data set. *Climate Dynamics*, 33(2–3), 341–364. https://doi.org/10.1007/s00382-008-0441-3

Large, W. G., & Yeager, S. G. (2004). *Diurnal to decadal global forcing for ocean and sea-ice models: The data sets and flux climatologies*. Tech. Rep. Boulder, CO, USA: National Center for Atmospheric Research Boulder.

Liu, Z. (2012). Dynamics of interdecadal climate variability: A historical perspective. *Journal of Climate*, 25(6), 1963–1995. https://doi.org/10.1175/2011jcli3980.1

Madec, G. (2012). *The NEMO team: Nemo ocean engine–Version 3.4. Note du Pôle de modélisation*. France: Institut Pierre-Simon Laplace (IPSL).

Marshall, J., Scott, J. R., Armour, K. C., Campin, J.-M., Kelley, M., & Romanou, A. (2015). The ocean’s role in the transient response of climate to abrupt greenhouse gas forcing. *Climate Dynamics*, 44(7–8), 2287–2299. https://doi.org/10.1007/s00382-014-2308-0

Mauritzen, C., Melsom, A., & Sutton, R. T. (2012). Importance of density-compensated temperature change for deep North Atlantic Ocean heat uptake. *Nature Geoscience*, 5(12), 905–910. https://doi.org/10.1038/ngeo1639

Meehl, G. A., Arblaster, J. M., Fasullo, J. T., Hu, A., & Trenberth, K. E. (2011). Model-based evidence of deep-ocean heat uptake during surface-temperature hiatus periods. *Nature Climate Change*, 1(7), 360–364. https://doi.org/10.1038/nclimate1229

Meehl, G. A., Goddard, L., Murphy, J., Stouffer, R. J., Boer, G., Danabasoglu, G., et al. (2009). Decadal prediction: Can it be skillful? *Bulletin of the American Meteorological Society*, 90(10), 1467–1486. https://doi.org/10.1175/2009bams2778.1

Newsom, E., Zanna, L., Khatiwala, S., & Gregory, J. M. (2020). The influence of warming patterns on passive ocean heat uptake. *Geophysical Research Letters*, 47, e2020GL088429. https://doi.org/10.1029/2020gl088429

Pillar, H. R., Heimbach, P., Johnson, H. L., & Marshall, D. P. (2016). Dynamical attribution of recent variability in Atlantic overturning. *Journal of Climate*, 29(9), 3339–3352. https://doi.org/10.1175/jcli-d-15-0727.1

Sérazin, G., Jaymond, A., Leroux, S., Penduff, T., Bessières, L., Llovel, W., et al. (2017). A global probabilistic study of the ocean heat content low-frequency variability: Atmospheric forcing versus oceanic chaos. *Geophysical Research Letters*, 44(11), 5580–5589. https://doi.org/10.1002/2017gl073026

Sévellec, F., Dijkstra, H. A., Drijfhout, S. S., & Germe, A. (2018). Dynamical attribution of oceanic prediction uncertainty in the North Atlantic: Application to the design of optimal monitoring systems. *Climate Dynamics*, 51(4), 1517–1535. https://doi.org/10.1007/s00382-017-3969-2

Sévellec, F., & Fedorov, A. V. (2013). The leading, interdecadal eigenmode of the Atlantic meridional overturning circulation in a realistic ocean model. *Journal of Climate*, 26(7), 2160–2183. https://doi.org/10.1175/jcli-d-11-00023.1

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

Smith, T., & Heimbach, P. (2019). Atmospheric origins of variability in the South Atlantic meridional overturning circulation. *Journal of Climate*, 32(5), 1483–1500. https://doi.org/10.1175/jcli-d-18-0311.1

Stephenson, D., Müller, S. A., & Sévellec, F. (2020). Tracking water masses using passive-tracer transport in NEMO v3.4 with NEMOT-AM: Application to North Atlantic Deep Water and North Atlantic Subtropical Mode Water. *Geoscientific Model Development*, 13(4), 2031–2050. https://doi.org/10.5194/gmd-13-2031-2020

Stephenson, D., & Sévellec, F. (2021). Dynamical attribution of North Atlantic interdecadal predictability to oceanic and atmospheric turbulence under diagnosed and optimal stochastic forcing. *Journal of Climate*, 1–69.

Trenberth, K. E., Fasullo, J. T., & Balmaseda, M. A. (2014). Earth’s energy imbalance. *Journal of Climate*, 27(9), 3129–3144. https://doi.org/10.1175/jcli-d-13-00294.1

Tziperman, E., Zanna, L., & Penland, C. (2008). Nonnormal thermohaline circulation dynamics in a coupled ocean–atmosphere gcm. *Journal of Physical Oceanography*, 38(3), 588–604. https://doi.org/10.1175/2007jpo3769.1
Vidard, A., Bouttier, P. A., & Vigilant, F. (2015). NEMOTAM: Tangent and adjoint models for the ocean modelling platform NEMO. *Geophysical Model Development*, 8(4), 1245–1257. https://doi.org/10.5194/gmd-8-1245-2015

Wang, G., & Dommenget, D. (2016). The leading modes of decadal sst variability in the southern ocean in cmip5 simulations. *Climate Dynamics*, 47(5–6), 1775–1792. https://doi.org/10.1007/s00382-015-2932-3

Xie, P., & Vallis, G. K. (2012). The passive and active nature of ocean heat uptake in idealized climate change experiments. *Climate Dynamics*, 38(3–4), 667–684. https://doi.org/10.1007/s00382-011-1063-8

Zanna, L., Khatiwala, S., Gregory, J. M., Ison, J., & Heimbach, P. (2019). Global reconstruction of historical ocean heat storage and transport. *Proceedings of the National Academy of Sciences of the United States of America*, 116(4), 1126–1131. https://doi.org/10.1073/pnas.1808838115

Zanna, L., & Tziperman, E. (2008). Optimal surface excitation of the thermohaline circulation. *Journal of Physical Oceanography*, 38(8), 1820–1830. https://doi.org/10.1175/2008jpo3752.1

Zhang, R. (2017). On the persistence and coherence of subpolar sea surface temperature and salinity anomalies associated with the Atlantic multidecadal variability. *Geophysical Research Letters*, 44(15), 7865–7875. https://doi.org/10.1002/2017gl074342