# Dynamical attribution of North Atlantic interdecadal predictability to oceanic and atmospheric turbulence under diagnosed and optimal stochastic forcing

Dafydd Stephenson, F. Sévellec

► To cite this version:

Dafydd Stephenson, F. Sévellec. Dynamical attribution of North Atlantic interdecadal predictability to oceanic and atmospheric turbulence under diagnosed and optimal stochastic forcing. Journal of Climate, 2021, pp.1-69. 10.1175/JCLI-D-21-0035.1. hal-03445776

HAL Id: hal-03445776

https://hal.science/hal-03445776

Submitted on 8 Aug 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Dynamical Attribution of North Atlantic Interdecadal Predictability to Oceanic and Atmospheric Turbulence under Diagnosed and Optimal Stochastic Forcing

DAFYDD STEPHENSONª AND F. SÉVELLECª,ᵇ

ª Ocean and Earth Science, University of Southampton, Southampton, United Kingdom

ᵇ Laboratoire d'Océanographie Physique et Spatiale, Univ. Brest CNRS IRD IFREMER, Brest, France

(Manuscript received 14 January 2021, in final form 23 July 2021)

**ABSTRACT:** Unpredictable variations in the ocean originate from both external atmospheric forcing and chaotic processes internal to the ocean itself, and are a crucial sink of predictability on interdecadal time scales. In a global ocean model, we present (i) an optimization framework to compute the most efficient noise patterns to generate uncertainty and (ii) a computationally inexpensive, dynamical method for attributing sources of ocean uncertainty to internal (mesoscale eddy-driven) and external (atmospherically driven) origins, sidestepping the more typical ensemble approach. These two methods are then applied to a range of metrics (heat content, volume transport, and heat transport) and time averages (monthly, yearly, and decadal) in the subtropical and subpolar North Atlantic. The optimal noise patterns create variability in integrated quantities of interest through features of the underlying circulation such as the North Atlantic Current and deep water formation regions. Meanwhile, noise forcing diagnosed from model representations of the actual climate system stimulates these theoretical patterns with various degrees of efficiency, ultimately leading to the growth of error. We reaffirm that higher-frequency variations in meridional transports are primarily wind driven, while surface buoyancy forcing is the ultimately dominant source of uncertainty at lower frequencies. For year-averaged quantities in the subtropics, it is mesoscale eddies that contribute the most to oceanic uncertainty, accounting for up to 60% after 60 years of growth for volume transport at 25°N. The impact of eddies is greatly reduced in the subpolar region, which we suggest may be explained by overall lower sensitivity to small-scale noise there.

**SIGNIFICANCE STATEMENT:** Climate does not change steadily; it naturally fluctuates around a general trend. The prediction of climate several decades to a century ahead depends mostly on the ability to anticipate future human activity, but for the coming years to a few decades ahead (when the future pathway of human activity is not yet fully apparent) natural fluctuations also have an important role. These fluctuations, however, cannot be perfectly predicted for long. The ability to predict them is limited, for example, by the build-up of unwelcome "noise" from erratic processes such as the weather. In this study, we look at the different sources of this noise, how important they are, and how they impact prediction accuracy of climatically important ocean quantities decades in the future. To achieve this, we use a unique computer simulation of the ocean, which works backward and describes how to most effectively create change. This uncovers the mechanisms by which the ocean most efficiently responds to noise, and also shows how this compares with the response to noise resembling that generated by the real ocean-atmosphere system. We demonstrate that in the climatically important region of the North Atlantic, unpredictable ocean circulation changes in the more tropical region are mostly due to oceanic mesoscale eddies (the oceanic equivalent of atmospheric storms). Farther north, in the subpolar region, however, it is the atmosphere that is primarily responsible for the development of oceanic prediction error.

**KEYWORDS:** North Atlantic Ocean; Meridional overturning circulation; Ocean circulation; Turbulence; Climate variability; Climate prediction; Interdecadal variability

# 1. Introduction

As the slow component of the climate system, the ocean is key to predicting variations on time scales of seasons or longer. However, the ocean is now known to exhibit substantial variability at all time scales. The predictability of these variations, and their attribution to different sources, is crucial to the understanding and prediction of climate, particularly on so-called

near-term time scales for which the anthropogenically forced signal is not yet dominant (Meehl et al. 2009).

Variations in the North Atlantic have long been hypothesized to be uniquely predictable due to interactions between its meridional overturning circulation (MOC) and anomalies in upper ocean heat content. In the late 1990s, an increase in computational resources allowed this hypothesis to be tested in state-of-the-art climate models using the prognostic technique of ensemble modeling [e.g., the review of Latif and Keenlyside (2011)]. In this framework, each member of a coupled climate model ensemble is initialized with a slightly perturbed atmospheric state. As the atmosphere has no predictability beyond a few weeks (Lorenz 1969), the atmospheric components of the ensemble rapidly diverge such that their differences in many ways resemble stochastic noise. The rate of divergence of the

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

*Corresponding author:* Dafydd Stephenson, d.stephenson@noc.soton.ac.uk
ocean components in response thus quantifies ocean predictability. Early studies using this methodology revealed enhanced predictability, often up to decades, in the North Atlantic sector against a background of strong MOC influence (Griffies and Bryan 1997; Grötzner et al. 1999; Collins and Sinha 2003; Msadek et al. 2010; Persechino et al. 2013). The implication that large-scale ocean dynamics slow error growth forced by the atmosphere is promising for near-term prediction in the region, but these studies collectively fail to account for oceanic mesoscale turbulence as an additional source of uncertainty. The relative importance of this source is becoming increasingly scrutinized as ocean components of cutting-edge climate models evolve toward eddying resolution (Kirtman et al. 2012; Small et al. 2014; Sein et al. 2018; Gutjahr et al. 2019; Held et al. 2019; Haarsma et al. 2020).

A new generation of studies is now addressing the question of attributing oceanic variability to internal (generated by chaotic oceanic processes) and external (atmospherically forced) origins using the prognostic ensemble approach in high-resolution ocean-only models (e.g., Sérazin et al. 2017; Leroux et al. 2018; Jamet et al. 2019). Each member has a common atmospheric forcing, but differing oceanic initial conditions. As such, the ensemble mean is taken to smooth out any intrinsic oceanic variability, such that its temporal variability is assumed to derive purely from fluctuations in the forcing. Contrarily, the ensemble spread, given their common atmospheric forcing, is assumed to come solely from intrinsic oceanic differences.

In this manner, Sérazin et al. (2017) conclude that ocean intrinsic variability is the dominant contributor to deep-ocean heat content fluctuations in the North Atlantic subtropical gyre and Gulf Stream regions, while Leroux et al. (2018) estimate that intrinsic MOC variability is 60% that of atmospheric at 26°N. In a regional model, Jamet et al. (2019) find that over half of the variability in the annually averaged Atlantic MOC at this latitude is intrinsic. Although oceanic variability forced at the domain boundaries will appear “external” in a regional model, this result agrees closely with the global model results of Grégorio et al. (2015). All of these studies show a shift in behavior at subpolar latitudes, where the atmospheric component dominates. It is unclear whether this shift arises from differing dynamics or the inability of eddy-permitting models to generate realistic turbulence in high-latitude regions (where deformation radii are greater), although this will likely become clear as fully eddy-resolving climate models become standard.

Despite the revolutionary advances in computing that now allow studies such as these to utilize ensembles containing as many as 50 members in a global, eddy-permitting ocean [as in Leroux et al. (2018)], such ensemble studies require national-scale computing infrastructure, inaccessible to many in the community. Furthermore, the ensemble approach does not allow a causal description of the translation of internal and external sources of unpredictable variability into expressed oceanic error growth or prediction uncertainty. An alternative framework, allowing dynamical attribution of the large-scale oceanic response to small perturbations (such as those from atmospheric fluxes or the mesoscale eddy field) is the adjoint method (Errico 1997). While the ensemble approach begins by applying small changes to the model state and then evaluates

their impact on oceanic metrics of interest, the adjoint method turns the problem around: it begins with an oceanic metric of interest and then describes its sensitivity to small changes in the model state.

This method has been applied to attributing Atlantic MOC fluctuations to different surface fluxes in the MITgcm by Pillar et al. (2016), and was used in the OPA model (the oceanic component of the model used herein) by Sévellec et al. (2018) to determine the relative impacts of unpredictable atmospheric forcing variability and initial condition uncertainty on the divergence of a theoretical ocean ensemble.

This study builds further on the theoretical ensemble approach of Sévellec et al. (2018). Here, we explore ocean error growth from two perspectives. In the first, we use an adjoint model to determine the most efficient patterns for stimulating ensemble divergence [the optimal stochastic perturbations (OSPs); Sévellec et al. 2007]. In this framework, the model is blind to sources of chaotic variability representative of the “real world” and instead describes how these sources should look in order to have the greatest effect on oceanic uncertainty. In this sense, the outcome describes, for different metrics, the sensitivity of their variance to different sources and locations, highlighting oceanic patterns of efficient error growth. The method shares much in common with the method of stochastic optimals (Farrell and Ioannou 1996), which has been previously applied in more idealized models of the ocean circulation such as box models (Lohmann and Schneider 1999; Tziperman and Ioannou 2002), a latitude-depth model (Sévellec et al. 2007), idealized coupled models (Chang et al. 2004; Zanna and Tziperman 2008), and a 3D planetary geostrophic model (Sévellec et al. 2009), typically in the context of stochastic excitation of internal modes. While the complexity of the problem has historically limited it to simple models, we are now able to determine both atmospheric and oceanic OSPs in a GCM under certain assumptions by separating the problem into smaller, regional OSP problems.

In the second perspective, we provide the model with stochastic representations of internal and external turbulent variability sources taken from model representations of the real climate system. This allows us to dynamically attribute ocean uncertainty to these different sources. The “realistic” sources are diagnosed from more complex models; the external, atmospheric component is calculated from a coupled non-eddying climate model, while the internal, mesoscale-eddy-driven component is calculated from an eddy-permitting ocean model. The attribution method is uniquely inexpensive—a single bidecadal simulation of a coupled climate model and of an eddy-permitting ocean model are used to compute the stochastic properties, while the highly efficient adjoint ocean model in a non-eddying (laminar) configuration can recreate a theoretically infinite ensemble with a single simulation (Sévellec and Sinha 2018).

The study proceeds as follows. In section 2, we outline the mathematical theory of stochastically forced ensembles that underlies our two approaches. This begins with a treatment of the classical, temporally uncorrelated (“white noise”) case, which provides the theoretical framework for deriving the OSPs. We then advance to time-correlated stochastic noise,
more appropriate for creating a more realistic representation of turbulence in the case of oceanic mesoscale eddies, although we highlight that our use of an eddy-permitting model, similarly to the global ensemble studies described above, means that a portion of the turbulent spectrum remains uncaptured by this representation. This choice was necessitated by computational cost, but the theory discussed here can be applied with a representation derived from any model. In section 3, we describe how this time-correlated representation is diagnosed, along with the three models used for the study and the configuration of our experiments. Our results are presented for both the optimal and diagnosed forcing cases in section 4 before being discussed along with our conclusions in section 5.

## 2. Theoretical framework: Variance of stochastically forced linear systems

While a more thorough exploration of the theory underpinning in this study follows in sections 2a–2c for interested readers, those who are interested in the key results will find them summarized in section 2d before being applied in section 4, and so can skip sections 2a–2c.

### a. Temporally uncorrelated forcing

A simple, classical model of low-frequency variability generation in the ocean is that of Hasselmann (1976). In it, mixed layer temperature changes are assumed to be a purely passive response to random, serially uncorrelated surface heat fluxes. These are absorbed and slowly “forgotten” by the ocean, which tends back toward its unperturbed state. The model is univariate and entirely determined by two parameters: the time scale on which this restoring occurs (parameterizing the ocean dynamics as a single memory term) and the volatility of the random fluxes (parameterizing the atmospheric forcing). It may be written as the stochastic differential equation

$$du = -\lambda u dt + \sigma dW, \quad (1)$$

which has the solution (for initial condition zero)

$$u(t_0, t_1) = \int_{t_0}^{t_1} e^{-\lambda(t_1-t)} \sigma dW(t), \quad (2)$$

where $u$ is the surface temperature, $t_0$ and $t_1$ are the initial and final time, $\sigma^2$ is the variance of temperature change induced by random surface atmospheric heat fluxes during a time increment $dt$, $\lambda^{-1}$ defines the $e$-folding time scale of the ocean dynamics (i.e., its memory), and $dW$ is an increment of a standard-normal Wiener process $W$ (akin to the distance of a random walk during the time increment $dt$). Equation (2) is thus an Itô integral (Itô 1944). It may be noted that the response, (2), is an Ornstein-Uhlenbeck process (Uhlenbeck and Ornstein 1930), a stochastic model that we will return to in section 2c.

Although a useful first-order representation of the evolution of unpredictable surface temperature variability (Frankignoul and Hasselmann 1977), the model is inherently limited by its treatment of a single forcing and response term, representing a spatial average of a single independent region of the ocean and

atmosphere (without accounting for any internal ocean processes, beyond a crude memory term). In a more realistic representation, atmospheric forcing may coherently influence multiple regions of the ocean, which may interact with each other through a range of variables and processes. If the dynamics of these interactions remain linear, (1) can be generalized to a nonautonomous linear system of stochastic differential equations:

$$d|\mathbf{u}\rangle = \mathbf{A}(t)|\mathbf{u}\rangle dt + \mathbf{L} d|\mathbf{W}(t)\rangle, \quad (3)$$

where $|\mathbf{u}\rangle$ is the ocean state vector anomaly, describing the response of each prognostic variable at each location, $|\mathbf{W}(t)\rangle$ is a vector of independent standard-normal Wiener processes, $\mathbf{A}(t)$ describes the linear interactions between all ocean variables and locations, and $\mathbf{L}$ is the lower-triangular matrix describing the stochastic atmospheric fluxes through the Cholesky decomposition $\Sigma = \mathbf{LL}^\dagger$ of their covariance matrix. In this decomposition, the dagger ($\dagger$) represents the adjoint defined by the Euclidean inner product. We have employed the “bra-ket” notation of Dirac (1939), in which vectors are represented as “kets” $|\mathbf{b}\rangle$ and covectors “bras” $\langle\mathbf{a}|$ such that closed bra-ket pairs become scalar through the Euclidean inner product $\langle\mathbf{a}|\mathbf{b}\rangle = c$, while open pairs $|\mathbf{a}\rangle\langle\mathbf{b}| = \mathbf{C}$ become matrices through the outer product. Although this is an improvement over the univariate Hasselmann model, there remain some simplifying assumptions about the nature of the noise forcing, which we now briefly reflect on. Foremost, it is taken that the noise forcing is Gaussian and additive, or state-independent. This is a common and effective assumption for large-scale-average quantities and climactic time scales (Sardeshmukh and Sura 2009) such as those considered here, but more complex processes beyond our study may be better represented by state-dependent or multiplicative noise, under which the mean state may “drift,” and Itô calculus (as used here) is an insufficient framework. An example is given by Sura and Sardeshmukh (2008), who note that high-frequency SST variations are often better modeled by a Hasselmann-type model in which the decay term ($\lambda$) itself is stochastic. For the interested reader, a thorough review of these types of considerations for stochastic climate modeling and their associated calculi is provided by Penland (2003).

We now address how to apply this stochastic theory of linear systems to an ocean model. Realistic ocean models are not linear, but the above theory can still be applied under certain simplifying assumptions. In particular, if perturbations (such as those applied by turbulent fluctuations) are sufficiently small, then a linearization of the model, representable in the same way as the complementary equation of (3) can provide a first-order description of its anomalous behavior in response to them. Consider a nonlinear system such as that solved by a typical ocean general circulation model (GCM):

$$d|\mathbf{U}\rangle = \mathcal{N}(|\mathbf{U}\rangle, t) dt,$$

where $\mathcal{N}$ is a nonlinear operator, $t$ is time, and $|\mathbf{U}\rangle$ the full state vector. This equation can be considered an abstraction wherein the operator $\mathcal{N}$ represents all processes, equations and boundary conditions numerically solved by the GCM.
Expansion of the full state vector $|\mathbf{U}\rangle = |\bar{\mathbf{u}}\rangle + |\mathbf{u}\rangle$ (about a mean state $|\bar{\mathbf{u}}\rangle$) yields

$$d(|\bar{\mathbf{u}}\rangle + |\mathbf{u}\rangle) = [\mathcal{N}(|\bar{\mathbf{u}}\rangle, t) + \mathbf{A}(t)|\mathbf{u}\rangle + \mathcal{O}(|\mathbf{u}\rangle^2)]dt. \quad (4)$$

Noting that $d|\bar{\mathbf{u}}\rangle = \mathcal{N}(|\bar{\mathbf{u}}\rangle, t) dt$ and neglecting higher-order terms leads to the complementary equation of (3). In this context, $\mathbf{A}(t)$ is the Jacobian of the nonlinear system with respect to the ocean state:

$$\mathbf{A}(t) = \frac{\partial}{\partial |\bar{\mathbf{u}}\rangle} \mathcal{N}(|\bar{\mathbf{u}}\rangle, t). \quad (5)$$

The (zero initial condition) solution to (3) is given by

$$|\mathbf{u}(t_0, t_1)\rangle = \int_{t_0}^{t_1} \Psi(t_1, t) \mathbf{L} d|\mathbf{W}(t)\rangle, \quad (6)$$

where $\Psi(t_1, t_0)$ is the propagator matrix [the scalar $\Psi(t_1, t_0) = e^{-\lambda(t_1-t_0)}$ in the univariate case of (2)] which describes the linear response of the ocean at time $t_1$ to changes originating from time $t_0$.

Beginning from the last formula, we can diagnose the covariance between any two scalar-valued metrics of the ocean state that are linear. These metrics can be defined by the covectors $|\mathbf{F}_{1,2}\rangle$, where the scalar products $\langle \mathbf{F}_{1,2}|\mathbf{u}\rangle = \langle \mathbf{u}|\mathbf{F}_{1,2}\rangle$ are the Euclidean inner products of the covectors and the ocean state vector anomaly. We have

$$\mathrm{Cov}(\langle \mathbf{F}_1 | \mathbf{u}(t_0, t_1) \rangle, \langle \mathbf{F}_2 | \mathbf{u}(t_0, t_1) \rangle) = E \left[ \left\langle \mathbf{F}_1 \left| \int_{t_0}^{t_1} \Psi(t_1, t) \mathbf{L} d|\mathbf{W}(t)\rangle \right. \right. \right. \\ \left. \left. \times \left\langle \mathbf{F}_2 \left| \int_{t_0}^{t_1} \Psi(t_1, s) \mathbf{L} d|\mathbf{W}(s)\rangle \right. \right. \right. \\ \left. \left. \left. \right\rangle \right] \right], \quad (7)$$

where $s$ represents time and $E[\cdot]$ is the expectation operator. A multidimensional generalization of Itô's isometry may be applied to this expression (e.g., section 3.6 of Duan and Wang 2014), leading to

$$\mathrm{Cov}(\langle \mathbf{F}_1 | \mathbf{u}(t_0, t_1) \rangle, \langle \mathbf{F}_2 | \mathbf{u}(t_0, t_1) \rangle) = \int_{t_0}^{t_1} \langle \mathbf{F}_1 | \Psi(t_1, t) \mathbf{\Sigma} \Psi^{\dagger}(t, t_1) | \mathbf{F}_2 \rangle dt, \quad (8)$$

which generalizes the result of Sévellec et al. (2018). Similarly to their approach, we remark that while it is standard to diagnose the variance evolution of a metric by propagating many realizations of (6) as an ensemble and considering its spread, (8) does not require us to propagate any such realization. Instead, it describes the response of such an ensemble (in the theoretical limit of large ensemble size) using only the statistical properties ($\mathbf{\Sigma}$) of the noise. It further provides a dynamical link between the response of the target metrics $\langle \mathbf{F}_{1,2} \rangle$ and the stochastic source of variability represented by $\mathbf{\Sigma}$. Where this representation can be linearly partitioned into independent sources (e.g., internal and external, $\mathbf{\Sigma} = \mathbf{\Sigma}_I + \mathbf{\Sigma}_E$), the variance can be dynamically attributed to each. The only requirements of the method are that

1) our metrics of interest $\langle \mathbf{F}_{1,2} |$ are linear functions of the ocean state;

2) we have a linear model of ocean dynamics, $\Psi(t_1, t_0)$ [we take a linearized OGCM, which following (4) is valid for small variations about a trajectory, see section 3]; and

3) we have an adequate statistical description $\mathbf{\Sigma}$ of any stochastic sources of variability.

Regarding the latter point, two approaches may be taken: the properties of the stochastic processes may be diagnosed and prescribed [as in, e.g., Sévellec et al. (2018)] or they may be determined from the linear model itself [in the framework of an optimization problem, as in, e.g., Sévellec et al. (2007, 2009)]. We begin with the latter approach, which provides insight into the mechanisms by which sources of variability are translated into oceanic variance in a theoretical setting.

### b. Optimal stochastic perturbations

As $\mathbf{\Sigma}$ can be allowed to take any form in (8), the problem of variance estimation can be reformulated as an optimization question: what form should $\mathbf{\Sigma}$ take such that variance $\mathrm{Var}[\langle \mathbf{F} | \mathbf{u}(t_0, t_1) \rangle] = \int_{t_0}^{t_1} \langle \mathbf{F} | \Psi(t_1, t) \mathbf{\Sigma} \Psi^{\dagger}(t, t_1) | \mathbf{F} \rangle dt$ is maximal for a given metric $\langle \mathbf{F} |$? The solution to the problem, under certain conditions, can be determined dynamically from the linear model itself, allowing insight into the mechanisms behind oceanic uncertainty without explicitly prescribing sources of uncertainty.

To determine the optimal $\mathbf{\Sigma}$, we apply two constraints to the optimal variance source: its global average has fixed amplitude, and any two points that are not independent have a correlation of $\pm 1$. The former implies that the stochastic process has finite power (corresponding to band-limited white noise), while the latter assumes that if two points covary, they must do so completely constructively (as would be optimal). We begin by considering the general case, where the stochastic process is partitioned into $N$ such regions (where each point in the region is perfectly correlated), before considering the specific cases corresponding to the two limits of $N$: (i) $N = 1$ corresponding to a fully global correlation [as in Sévellec et al. (2007, 2009)] and (ii) $N = n$ (where $n$ is the dimension of the state vector $|\mathbf{u}\rangle$), corresponding to the absence of any correlation.

#### 1) GENERAL CASE

As outlined above, we partition the stochastic process into $N$ regions such that points within the regions are perfectly varying, but are independent of points in other regions. Equivalently, we separate $\mathbf{\Sigma} \in \mathbb{R}^{n \times n}$ into $N$ local matrices $\mathbf{\Sigma}_i \in \mathbb{R}^{m_i \times m_i}$ (where $m_i$ is the local dimension of the $i$th region), and define a binary projection $\mathbf{B}_i \in \mathbb{R}^{n \times m_i}$ such that

$$\mathbf{\Sigma} = \sum_{i=1}^{N} \mathbf{B}_i \mathbf{\Sigma}_i \mathbf{B}_i^{\dagger}. \quad (9)$$

Following (6), the evolution of the state vector in response to stimulation in the $i$th region is

$$|\mathbf{u}_i(t_0, t_1)\rangle = \int_{t_0}^{t_1} \Psi(t_1, t) \mathbf{B}_i \mathbf{L}_i d|\mathbf{W}_i(t)\rangle, \quad (10)$$

where $\mathbf{\Sigma}_i = \mathbf{L}_i \mathbf{L}_i^{\dagger}$ is the Cholesky decomposition of the local covariance matrix, equivalently to the global case. Fundamentally,
as the region is perfectly correlated, it may be written in terms of a single stochastic process. The vector $\mathbf{L}_i d|\mathbf{W}_i$ thus becomes $|\mathbf{L}_i\rangle dW_i$, such that $\mathbf{\Sigma}_i$ is the outer product $\mathbf{\Sigma}_i = |\mathbf{L}_i\rangle\langle\mathbf{L}_i|$. The implication is that in the region, a single Wiener process is “shaped” by a pattern of local amplitudes $|\mathbf{L}_i\rangle$.

To determine the optimal shape of this pattern, we utilize the method of Lagrange multipliers [consistently with Sévellec et al. (2007)]. In particular, we wish to maximize the local contribution to the variance

$$ \mathrm{Var}(\langle \mathbf{F} | \mathbf{u}_i \rangle) = \int_{t_0}^{t_1} \langle \mathbf{F} | \Psi(t_1, t) \mathbf{B}_i \mathbf{\Sigma}_i \mathbf{B}_i^\dagger \Psi(t, t_1) | \mathbf{F} \rangle dt \quad (11a) $$

under the constraint that the amplitude $\epsilon_i$ of $\mathbf{\Sigma}_i$ follows

$$ \mathrm{Tr}(\mathbf{S}_i \mathbf{\Sigma}_i) = \langle \mathbf{L}_i | \mathbf{S}_i | \mathbf{L}_i \rangle = \epsilon_i^2, \quad (11b) $$

where $\mathbf{S}_i \in \mathbb{R}^{m_i \times m_i}$ is a (diagonal) volumetric weighting matrix. The corresponding Lagrange function can be expressed as

$$ \mathcal{L}(\gamma_i, |\mathbf{L}_i\rangle, t_0, t_1) = \int_{t_0}^{t_1} \langle \mathbf{F} | \Psi(t_1, t) \mathbf{B}_i | \mathbf{L}_i \rangle^2 dt - \gamma_i (\langle \mathbf{L}_i | \mathbf{S}_i | \mathbf{L}_i \rangle - \epsilon_i^2), \quad (12) $$

where the scalar $\gamma_i$ is the Lagrange multiplier. Maximizing the Lagrangian leads to

$$ \left. \frac{\partial \mathcal{L}}{\partial |\mathbf{L}_i|} \right|_{\{\gamma_i^*, |\mathbf{L}_i^*\}\}} = 0, \\ \int_{t_0}^{t_1} (\mathbf{B}_i^\dagger \Psi^\dagger(t, t_1) |\mathbf{F}\rangle \langle \mathbf{F} | \Psi(t_1, t) \mathbf{B}_i) dt |\mathbf{L}_i^*\rangle - \gamma_i^* \mathbf{S}_i |\mathbf{L}_i^*\rangle = 0, \quad (13) $$

which holds when $\gamma_i^*$ and $|\mathbf{L}_i^*\rangle$ are an eigenvalue-eigenvector pair of

$$ \mathbf{S}_i^{-1} \int_{t_0}^{t_1} (\mathbf{B}_i^\dagger \Psi^\dagger(t, t_1) |\mathbf{F}\rangle \langle \mathbf{F} | \Psi(t_1, t) \mathbf{B}_i^\dagger) dt, \quad (14) $$

since $\mathbf{S}_i$ (as an operator representing a norm) is invertible. Any such eigenpair represents a particular solution to the optimization problem. All together these solutions are a set of “stochastic optimals” as described by Farrell and Ioannou (1996). While our framework is similar to that of Farrell and Ioannou (1996), we are optimizing a variance measure, rather than a norm, similarly to Tziperman and Ioannou (2002). Of this set of solutions, we seek the solution with the greatest effect. We note that left multiplication of (13) by $(\mathbf{L}_i^* | \mathbf{S}_i)$ results in

$$ \int_{t_0}^{t_1} \langle \mathbf{F} | \Psi(t_1, t) \mathbf{B}_i | \mathbf{L}_i^* \rangle^2 dt = \gamma_i^* \langle \mathbf{L}_i^* | \mathbf{S}_i | \mathbf{L}_i^* \rangle, $$

or, equivalently, $\mathrm{Var}[(\mathbf{F}|\mathbf{u}_i(t_0, t_1))] = \gamma_i^{*2}$, so that the Lagrangian multiplier $\gamma_i$ is essentially representing the variance that we wish to maximize. Hence the eigenvector $|\mathbf{L}_i^{\mathrm{opt}}\rangle$ corresponding to the universally optimal solution of (13) is that belonging to the leading eigenvalue $\gamma_i^{\mathrm{opt}}$. This vector describes the pattern of the local OSP, and it remains only to rescale it by a factor $\epsilon_i \langle \mathbf{L}_i^{\mathrm{opt}} | \mathbf{S}_i | \mathbf{L}_i^{\mathrm{opt}} \rangle^{-1/2}$ so that the associated optimal covariance matrix obeys the constraint (11b). This matrix is therefore

$$ \mathbf{\Sigma}_i^{\mathrm{opt}} = \epsilon_i^2 \frac{|\mathbf{L}_i^{\mathrm{opt}}\rangle \langle \mathbf{L}_i^{\mathrm{opt}}|}{\langle \mathbf{L}_i^{\mathrm{opt}} | \mathbf{S}_i | \mathbf{L}_i^{\mathrm{opt}} \rangle}. \quad (15) $$

Our local magnitude $\epsilon_i$ may be chosen arbitrarily, and so, although the $N$ regions correspond to $N$ independent problems, we seek an optimal scaling $\epsilon_i$ that maximizes their individual contribution to the overall variance, while constraining the total magnitude $\sum_{i=1}^N \epsilon_i^2 = \epsilon^2$. In particular, we note that the total variance $\mathrm{Var}[(\mathbf{F}|\mathbf{u}(t_0, t_1))] = \sum_{i=1}^N \epsilon_i^2 \gamma_i^{\mathrm{opt}}$ following the above. This may be alternatively rewritten as an inner product $\mathrm{Var}[(\mathbf{F}|\mathbf{u}(t_0, t_1))] = \langle \mathbf{E} | \gamma \rangle$, where $|\mathbf{E}\rangle$ and $|\gamma\rangle$ are vectors of dimension $N$ concatenating all the amplitudes $(\epsilon_i^2)$ and optimal variances $(\gamma_i^{\mathrm{opt}})$, respectively, of the local optimal shape $(|\mathbf{L}_i^{\mathrm{opt}}\rangle)$ for the $N$ regions. As the inner product is maximal for parallel vectors (i.e., $|\mathbf{E}\rangle$ parallel to $|\gamma\rangle$), it follows after some algebra that

$$ \epsilon_i^2 = \frac{\epsilon_i^2 \gamma_i^{\mathrm{opt}}}{\sum_{i=1}^{N} \gamma_i^{\mathrm{opt}}}. \quad (16) $$

Hence, for these choices of $\epsilon_i$, we obtain the global covariance matrix describing the overall solution to the OSP problem in terms of local OSPs:

$$ \mathbf{\Sigma}_i^{\mathrm{opt}} = \frac{\epsilon_i^2}{\sum_{i=1}^{N} \gamma_i^{\mathrm{opt}}} \sum_{i=1}^{N} \gamma_i^{\mathrm{opt}} \mathbf{B}_i \frac{|\mathbf{L}_i^{\mathrm{opt}}\rangle \langle \mathbf{L}_i^{\mathrm{opt}}|}{\langle \mathbf{L}_i^{\mathrm{opt}} | \mathbf{S}_i | \mathbf{L}_i^{\mathrm{opt}} \rangle} \mathbf{B}_i^\dagger, \quad (17) $$

where, as described above, $|\mathbf{L}_i^{\mathrm{opt}}\rangle$ and $\gamma_i^{\mathrm{opt}}$ are the leading eigenpair of

$$ \mathbf{S}_i^{-1} \int_{t_0}^{t_1} \mathbf{B}_i \Psi^\dagger(t, t_1) |\mathbf{F}\rangle \langle \mathbf{F} | \Psi(t_1, t) \mathbf{B}_i^\dagger dt. $$

## 2) LIMITING CASES

The above derivation applies to the case of $N$ perfectly correlated independent regions, but we may consider two specific cases of this in order to imitate conditions similar to the atmospherically forced and eddy-driven variability felt by the ocean. In particular, we consider the two limiting cases: $N = 1$ and $N = \mathrm{dim}(|\mathbf{u}\rangle)$. The former case, where the forcing is everywhere perfectly correlated, can be applied to the surface layer as an idealized representation of the large-scale coherent patterns of the atmosphere (Sévellec et al. 2007, 2009). The latter case, where the forcing is uncorrelated between all variables and locations, is taken as an idealized representation of small-scale noise in the ocean (i.e., noise induced by subgrid processes). These cases correspond to solving a single eigenvalue problem versus solving $\mathrm{dim}(|\mathbf{u}\rangle)$ (trivially scalar) eigenvalue problems. In particular, for $N = 1$, the sole projection matrix is the identity matrix $\mathbf{B}_1 = \mathbf{I}$, while for $N = n = \mathrm{dim}(|\mathbf{u}\rangle)$, the projection matrices become the standard basis vectors $\mathbf{B}_i = |e_i\rangle$ (i.e., $|e_i\rangle$ projects a scalar to the $i$th location of the full state vector).

In the former (everywhere perfectly covarying) case, (17) becomes

$$ \mathbf{\Sigma}_i^{\mathrm{opt}} = \epsilon_i^2 \frac{|\mathbf{L}_i^{\mathrm{opt}}\rangle \langle \mathbf{L}_i^{\mathrm{opt}}|}{\langle \mathbf{L}_i^{\mathrm{opt}} | \mathbf{S}_i | \mathbf{L}_i^{\mathrm{opt}} \rangle}, \quad (18) $$
where $|\mathbf{L}^{\text{opt}}\rangle$ is the leading eigenvector of

$$ \mathbf{S}^{-1} \int_{t_0}^{t_1} \Psi^{\dagger}(t, t_1) |\mathbf{F}\rangle \langle \mathbf{F}| \Psi(t_1, t) dt. $$

The latter (everywhere uncorrelated) case corresponds to the condition that every point is independent, and $\mathbf{\Sigma}^{\text{opt}}$ is diagonal. The associated eigen“vector” problems are scalar, such that the eigenspace is infinite. All terms in (17) are now scalars such that $|\mathbf{L}_i^{\text{opt}}\rangle$ can be seen to cancel, while the matrices $\mathbf{S}_i$ may be written as $S_i$. Ultimately,

$$ \mathbf{\Sigma}^{\text{opt}} = \frac{\epsilon^2}{\sum_{i=1}^{N} \gamma_i^{\text{opt}}} \sum_{i=1}^{N} |\mathbf{e}_i\rangle \frac{\gamma_i^{\text{opt}}}{S_i} \langle \mathbf{e}_i|, \quad (19) $$

where, solving (13) with $\mathbf{B}_i = |\mathbf{e}_i\rangle$, the eigenvalues $\gamma_i^{\text{opt}}$ are trivially the diagonal elements of

$$ \mathbf{S}^{-1} \int_{t_0}^{t_1} \Psi^{\dagger}(t, t_1) |\mathbf{F}\rangle \langle \mathbf{F}| \Psi(t_1, t) dt. $$

The sum of the eigenvalues is also the trace of this (scaled outer product) matrix, and is thus given by the corresponding inner product. Therefore, from (19) the optimal stochastic covariance matrix in the completely uncorrelated case is

$$ \mathbf{\Sigma}^{\text{opt}} = \frac{\epsilon^2}{\int_{t_0}^{t_1} \langle \mathbf{F}|\Psi(t_1, t)\mathbf{S}^{-1}\Psi^{\dagger}(t, t_1)|\mathbf{F}\rangle dt} \times \operatorname{diag} \left[ \mathbf{S}^{-1} \int_{t_0}^{t_1} \Psi^{\dagger}(t, t_1) |\mathbf{F}\rangle \langle \mathbf{F}|\Psi(t_1, t) dt \right] \mathbf{S}^{-1} \quad (20) $$

(where the diag[·] operator corresponds to the diagonal matrix with the same diagonal). We respectively use these two limiting cases to explore theoretical variance linked to idealized atmospheric forcing (assuming perfect correlation everywhere over the surface and zero noise in the interior) and ocean internal subgrid fluxes (assuming noise everywhere, with zero correlation between locations and variables).

A useful metric of the OSP is the ratio of the output variance to the input variance $A_* = \text{Var}[(\mathbf{F}|\mathbf{u}(t_0, t_1))/\epsilon^2]$, which we term the response ratio. Notably, for the globally perfect covariance case, this is simply the associated eigenvalue

$$ A_* = \gamma^{\text{opt}}. \quad (21a) $$

For the globally decorrelated case,

$$ A_* = \int_{t_0}^{t_1} \langle \mathbf{F} | \Psi(t_1, t) \mathbf{S}^{-1} \Psi^{\dagger}(t, t_1) | \mathbf{F} \rangle dt \quad (21b) $$

is the sum of the eigenvalues.

### c. Temporally correlated forcing

Our considerations so far have involved stochastic forcing with varying levels of spatial coherence, but which is serially decorrelated (therefore band-limited white noise). While this allows an idealized, theoretical exploration of variance generation mechanisms in the optimal case, it is inadequate for

realistically representing turbulent fluxes in the climate system, as we wish to in the diagnosed case. Indeed, while white noise is typically considered an acceptable representation of atmospheric variability (which decays on time scales much shorter than those of the oceanic large scale; Hasselmann 1976), the ocean mesoscale eddy field evolves much more slowly (e.g., Chelton et al. 2007). To represent this using diagnosed fluxes, we therefore extend our framework to include temporally correlated stochastic forcing. We consider again the Ornstein–Uhlenbeck case, which is a simple example of a temporally correlated stochastic process.

We begin by modifying (3) such that anomalous fluxes are now modeled by a continuous, time-integrable stochastic process (contrary to the former, white noise case, where they were everywhere discontinuous and representable only in the framework of distribution theory). The equation becomes

$$ d|\mathbf{u}\rangle = [\mathbf{A}(t)|\mathbf{u}\rangle + |\mathbf{X}(t)\rangle] dt, \quad (22) $$

where, as before, $|\mathbf{u}\rangle$ defines the state vector anomaly, $\mathbf{A}$ defines the system's linear dynamics [for instance via the Jacobian of a corresponding nonlinear system, as in (5)], and, in contrast to the previous cases, $|\mathbf{X}\rangle$ is the forcing from continuous, spatiotemporally correlated stochastic processes. The zero-initial-condition solution is given by

$$ |\mathbf{u}(t_0, t_1)\rangle = \int_{t_0}^{t_1} \Psi(t_1, t) |\mathbf{X}(t)\rangle dt, \quad (23) $$

where the complementary equation and therefore the propagator matrix $\Psi(t_1, t_0)$ are notably identical to (6). As in (7), we seek the covariance between two metrics of the state vector, given by

$$ \operatorname{Cov}[\langle \mathbf{F}_1 | \mathbf{u}(t_0, t_1) \rangle, \langle \mathbf{F}_2 | \mathbf{u}(t_0, t_1) \rangle] = \int_{t_0}^{t_1} \int_{t_0}^{t_1} \langle \mathbf{F}_1 | \Psi(t_1, t) E[|\mathbf{X}(t)\rangle] \times \langle \mathbf{X}(s) | \Psi^{\dagger}(s, t_1) | \mathbf{F}_2 \rangle dt ds, \quad (24) $$

where the term $E[|\mathbf{X}(t)\rangle\langle\mathbf{X}(s)|]$ gives the spatiotemporal covariance matrix of the forcing. In the white noise case, the autocorrelation conceptually corresponds to the Dirac delta function, leading to $E[|\mathbf{X}(t)\rangle\langle\mathbf{X}(s)|] = \delta(t-s)\mathbf{LL}^{\dagger}$, consistently with (8). For a vector $|\mathbf{X}\rangle$ of saturated Ornstein–Uhlenbeck processes, the autocorrelation function is given by (see, e.g., section 1.2 of Lindner 2009)

$$ E[|\mathbf{X}(t)\rangle\langle\mathbf{X}(s)|] = e^{-\lambda t} \mathbf{LL}^{\dagger} e^{-\lambda s}, \quad (25) $$

where $\lambda$ is a diagonal matrix of reciprocal $e$-folding times of the anomalous fluxes at each location, and $\mathbf{LL}^{\dagger} = \mathbf{\Sigma}$ is their spatial covariance matrix. As these quantities can be diagnosed from an appropriate dataset, we can use this formulation to diagnose the variance growth.

In the proceeding section we diagnose (from more realistic models that are assumed to capture a substantial portion of atmospheric and oceanic turbulence) $\lambda$ and $\mathbf{\Sigma}$ for the cases of external (atmospheric; $\lambda_E$, $\mathbf{\Sigma}_E$) and internal (oceanic mesoscale eddy driven; $\lambda_I$, $\mathbf{\Sigma}_I$) turbulent fluxes, assessing the appropriateness of the Ornstein–Uhlenbeck representation. We then proceed to attribute the variance of different metrics in response to these
sources, which, following (24) and assuming independence between the internal and external components is given by

$$ \mathrm{Var}[\langle \mathbf{F} | \mathbf{u}(t_0, t_1) \rangle] = \int_{t_0}^{t_1} \int_{t_0}^{t_1} \langle \mathbf{F} | \Psi(t_1, t) e^{-\lambda_t^i \Sigma_i} e^{-\lambda_i^j s} \Psi^\dagger(s, t_1) | \mathbf{F} \rangle \, dt \, ds \\ + \int_{t_0}^{t_1} \int_{t_0}^{t_1} \langle \mathbf{F} | \Psi(t_1, t) e^{-\lambda_E^i \Sigma_E} e^{-\lambda_E^j s} \Psi^\dagger(s, t_1) | \mathbf{F} \rangle \, dt \, ds. \quad (26) $$

The variance may be broken down further still, by writing the covariance matrices as the sum of their different components. For example, we are interested in the independent contributions of buoyancy and momentum fluxes to the externally forced component Var_E of the variance (corresponding to the $\lambda_E$, $\Sigma_E$ terms), and, in the latter case, the separate contributions of the covarying zonal and meridional momentum fluxes. The final term of (26) can accordingly be split into

$$ \mathrm{Var}_E[\langle \mathbf{F} | \mathbf{u}(t_0, t_1) \rangle] = \int_{t_0}^{t_1} \int_{t_0}^{t_1} \langle \mathbf{F} | \Psi(t_1, t) \exp(-\lambda_E^b t) \Sigma_E^b \exp(-\lambda_E^{b\dagger} s) \Psi^\dagger(s, t_1) | \mathbf{F} \rangle \, dt \, ds, \\ + \int_{t_0}^{t_1} \int_{t_0}^{t_1} \langle \mathbf{F} | \Psi(t_1, t) \exp(-\lambda_E^u t) \Sigma_E^u \exp(-\lambda_E^{u\dagger} s) \Psi^\dagger(s, t_1) | \mathbf{F} \rangle \, dt \, ds, \\ + \int_{t_0}^{t_1} \int_{t_0}^{t_1} \langle \mathbf{F} | \Psi(t_1, t) \exp(-\lambda_E^v t) \Sigma_E^v \exp(-\lambda_E^{v\dagger} s) \Psi^\dagger(s, t_1) | \mathbf{F} \rangle \, dt \, ds, \\ + \int_{t_0}^{t_1} \int_{t_0}^{t_1} \langle \mathbf{F} | \Psi(t_1, t) \exp(-\lambda_E^u t) \Sigma_E^{u\dagger} \exp(-\lambda_E^{u\dagger} s) \Psi^\dagger(s, t_1) | \mathbf{F} \rangle \, dt \, ds, \quad (27) $$

where $(\lambda_E^{b,u,v}, \Sigma_E^{b,u,v})$ are the external noise properties for the buoyancy fluxes and for the zonal and meridional momentum fluxes, respectively; $\Sigma_E^{u,v}$ is for the zonal and meridional covariance term.

Finally, in addition to separating the variance into contributions from different variables, we note that we can also isolate contributions from different regions of space. The inner products of (26) represent spatial integrals of local contributions to the total variance (integrated over volume in the internal case and over area in the external case). An alternative formulation of (26) is therefore

$$ \mathrm{Var}[\langle \mathbf{F} | \mathbf{u}(t_0, t_1) \rangle] = \int_{\Omega} \mathcal{V}_I(x, y, z, t_0, t_1) \, dV \\ + \int_{\Omega_0} \mathcal{V}_E(x, y, t_0, t_1) \, dA, \quad (28) $$

where $\mathcal{V}_I$ and $\mathcal{V}_E$ are continuous functions representing the respective internal variance contribution per unit volume and external variance contribution per unit area, $\Omega$ and $dV$ represent the ocean interior and a volume increment, respectively, $\Omega_0$ and $dA$ represent the ocean surface and an area increment, respectively, and $x, y, z$ are the zonal, meridional, and vertical coordinates. The corresponding integrands are thus spatial distributions of variance contributions. This can be applied to both (26) and (27) without loss of generality.

### d. Summary of theoretical framework

We have shown that, by representing the ocean as a linear system forced by temporally uncorrelated white Gaussian noise (WGN) in (6), certain properties may be derived analytically given only the properties of the noise (its spatial coherence and intensity, represented by the covariance matrix $\Sigma$). In particular, we calculated the covariance of any two scalar-valued linear metrics of the ocean state (and, as a corollary, the variance of a single such metric, for instance ocean heat content) using just

this linear representation and noise covariance matrix in (8). In section 2b, we considered one option for calculating a noise covariance matrix, seeking the noise distribution that most efficiently stimulates variance in a metric of interest, using a formulation of the Lagrange multiplier method for calculating stochastic optimals. As the covariance matrix describes both the magnitude and coherence of the noise, we were able to solve this problem for a given noise amplitude in the case of both a globally coherent noise (with every location perfectly correlated or anticorrelated) in (19) and an everywhere decorrelated noise in (20), with each case transforming the input variance of the noise into an output variance in the target metric via a linear rescaling (“response ratio”) given in (21a) and (21b), respectively. Finally, in addition to the consideration of spatial correlation in sections 2a and 2b, we introduced temporally correlated noise in section 2c, motivated by the slower evolution of the ocean’s mesoscale eddy field and its poor agreement with WGN. We found that by introducing a second matrix $\lambda$ corresponding to the decorrelation time of random fluxes in the system, we could diagnose the variance of a metric of choice similarly to section 2a, and that this could be decomposed into contributions from oceanic (internal) and atmospheric (external) noise in (26).

## 3. Model configurations, methods, and experimental design

### a. Linear ocean model configuration

As outlined in section 2, we use a linear ocean model to provide the propagator matrix $\Psi$, which is used both to derive our OSPs [following (17)] and to evolve our prescribed, diagnosed stochastic processes [following (24)]. The model is v3.4 of the NEMO GCM (Madec et al. 2012) whose routines are linearized in the tangent-linear and adjoint model (TAM) package NEMOTAM (Vidard et al. 2015). The model is run in the nominal 2° ORCA2 configuration with 31 vertical levels in
partial-step z coordinates, subject to repeated CORE normal year forcing (Large and Yeager 2004). More details can be found in Stephenson et al. (2020). We note that the same ocean model is common to our linear propagator $\Psi$, the coupled climate model (section 3b) used to diagnose our stochastic external flux representation ($\lambda_E$, $\Sigma_E$) and the high-resolution ocean model (section 3c) used to diagnose our stochastic internal flux representation ($\lambda_I$, $\Sigma_I$). We therefore consider seasonal variations of the oceanic large scale to be common to all three cases, which are explicitly captured in the trajectory $|\bar{u}|$ of (4). In this sense, our stochastic representations are of anomalies $[|\bar{u}|]$ in (4)] from this shared climatology, which are unresolved in the low-resolution, ocean-only model. We note that this intermodel approach is subject to the same approximation as other studies which diagnose forcings from higher-fidelity models (or indeed observations): the ocean circulation and vertical structure is assumed to be the same in all cases. Projecting from one model to another will thus naturally lead to misrepresentations around technically challenging features (such as the Gulf Stream separation, for example), or areas where isopycnals do not align between models. The impact of such mismatches, similarly to that of nonstationary eddy variations such as seasonal cycles, is taken to be small (Fig. 1 shows the climatologies of the models from which these anomalies are taken).

## b. Diagnosis of stochastic atmospheric fluxes

To represent the effects of anomalous surface fluxes by an Ornstein–Uhlenbeck process we diagnose the parameters $\lambda_E$, $\Sigma_E$ from the outputs of a coupled climate model. In particular, we use the IPSL-CM5A-LR coupled model, which was run for 20 years in its CMIP5 preindustrial control configuration (cf. Dufresne et al. 2013). The ocean component of the model is NEMO-ORCA2 (v3.2) which has the same (ORCA2) configuration as our linearized ocean model. To cleanly separate atmospherically forced variability from internally forced turbulent ocean variability (which is diagnosed separately; section 3c), the ocean component of the chosen climate model is laminar, such that there is no internal turbulent variability (Grégorio et al. 2015). The atmospheric component is the LMDZ5a model, with a horizontal resolution of 3.75° × 1.9° and 39 levels in the vertical (Hourdin et al. 2013).

To isolate the impact of external forcing, the 20-yr time series of daily averaged surface zonal ($\tau_x$) and meridional ($\tau_y$) wind stress, heat ($Q$) and freshwater fluxes (FWF) produced by the coupled model were considered. As described in section 3a, the climatologies of these fluxes were taken to be present in the trajectory of the linear model (via its repeated annual forcing) and so were removed. The remaining anomalies were then linearly mapped to a corresponding external-flux-induced rate of change in ocean surface zonal and meridional velocity [$F_u^E = \tau_x/(\rho_0\Delta z)$ and $F_v^E = \tau_y/(\rho_0\Delta z)$, respectively, where $\Delta z$ is the uppermost model layer thickness and $\rho_0$ is a reference density], sea surface temperature (SST; $F_u^E = Q/(\rho_0 c_p \Delta z)$, where $c_p$ is specific heat capacity), and sea surface salinity (SSS; $F_s^E = S_0 \text{FWF}/(\rho_0 \Delta z)$, where $S_0$ is a reference salinity). The covariance and e-folding decorrelation time of these time series (Figs. 2a–e) were then used to construct the stochastic representation following (25).

The variance of the heat flux term (Fig. 2a) is broadly distributed away from the tropics with regions of intense focus such as western boundary currents, while the freshwater flux variance term is conversely highest in the tropics (Fig. 2b). Their covariance (Fig. 2c) reflects this difference such that increasing $F_T^E$ corresponds to salinification in these regions of highest variance in $F_T^E$ and freshening in regions of highest $F_s^E$ variance. Both temperature and salinity changes are most persistent at low latitudes (Figs. 2a,b, contours). For wind stress-induced surface velocity changes, zonal and meridional variances show broadly similar spatial patterns, focused at high latitudes (Figs. 2d and 2e, respectively). The zonal component is notably more intense and more persistent (Figs. 2d,e, contours).

The matrix $\Sigma_E$ was populated using the covariances of these time series with the corresponding time series of each dependent variable at every other location (Fig. 2 shows the lead diagonal of $\Sigma_E$). The term $\lambda_E$ is a diagonal matrix of local e-folding times calculated from the lag-autocorrelation of the time series (shown by contours in Fig. 2). Buoyancy and momentum fluxes were assumed independent of each other, but their components (temperature and salinity for the former, meridional and zonal momentum for the latter) are allowed to spatially covary.

To evaluate the goodness of fit of the Ornstein–Uhlenbeck process representation, we compare the power spectral density (PSD) of a theoretically perfect process with matching parameters at each location with the PSD produced by the time series. To evenly weight all frequencies, we use the root-mean-square logarithmic error (RMSLE) metric, normalized by the mean of the logarithm of the PSD. This effectively corresponds to taking a normalized root-mean-square error, but in logarithmic space, such that all frequencies contribute uniformly. For comparison we also evaluate the error in the same way when the more traditional Gaussian white noise representation (i.e., constant PSD) is used to fit the model outputs (Figs. 3b,c,e,f). This reveals that while a white noise representation fits the spectrum of surface fluxes very well, the Ornstein–Uhlenbeck model is still almost everywhere an improvement (Figs. 3h,j,k,l). For surface temperature and zonal momentum fluxes (Figs. 3b,c), the improvement is universal, but more pronounced in the tropics, where decay times are longer (Figs. 2a,d). This tropical improvement is also evident in surface freshwater and meridional momentum flux (Figs. 3e,f), but with a very slight increase in error in the latter in a small region west of the Agulhas and in the former in the Southern Ocean and subpolar North Atlantic. We note that this error increase is localized to regions where fluxes are weak (Figs. 2b,e), but the error metric is not weighted by the intensity of the fluxes, and is based purely on the shape of the power spectrum. The error increase is furthermore far smaller in magnitude than the corresponding error decrease elsewhere.

## c. Diagnosis of ocean mesoscale eddy fluxes

In addition to the variability driven by turbulent atmospheric processes, processes creating variability exist within the ocean interior that are also unresolved by our laminar ocean-only model, due to the coarseness of its spatial discretization. To represent this, we utilize spatiotemporal Reynolds averaging,
<figure><img src="image_2.png" /><figcaption>FIG. 1. Climatological properties of the models used to diagnose our stochastic forcing. (left) Mean (a) heat and (b) freshwater fluxes from the two-decade coupled model run used to diagnose the external stochastic forcing. Arrows show the strength of the corresponding mean wind stress. (middle) Mean (c) sea surface temperature and (d) salinity from the LRM, which provides the linearized model propagator, averaged over the 60-yr trajectory run. Arrows show the strength of the surface velocity field. (right) Mean (e) sea surface temperature and (f) salinity from the HRM from which the internal flux representation is calculated. Arrows again show the strength of the surface velocity field.</figcaption></figure>

in which large-scale temperature variations are potentially impacted by small-scale anomalies in a purely advective transport framework. For the temperature, the associated advection equation (at high Péclet number, such that diffusive processes can be neglected) reads

$$ \partial_t T + \langle \mathbf{V} | \nabla T \rangle = 0, \\ \partial_t (\hat{T} + \bar{T}) + \langle \hat{\mathbf{V}} + \bar{\mathbf{V}} | \nabla (\hat{T} + \bar{T}) \rangle = 0, \quad (29) $$

where T and V are the scalar and tridirectional vector fields of temperature and of velocity, respectively, $\nabla(\cdot)$ is the tridirectional

gradient operator, $\langle \cdot | \cdot \rangle$ is the inner product, $\hat{\cdot}$ is a tridirectional spatial averaging operator, and $\bar{\cdot}$ is its associated spatial fluctuation. This separation is such that the lower-resolution model (LRM) is able to resolve temperatures at the scale of the spatial average (e.g., $\hat{T}$), while the higher-resolution model (HRM) resolves the sum of the spatial average and its fluctuation (e.g., $T = \hat{T} + \bar{T}$). We are interested in the mean effect of the small scale on the large scale following application of the spatial averaging operator. Applying this operator, the equation reduces to

$$ \partial_t \hat{T} + \langle \hat{\mathbf{V}} | \nabla \hat{T} \rangle = - \langle \hat{\mathbf{V}} | \nabla \bar{T} \rangle. \quad (30) $$
<figure><img src="image_3.png" /><figcaption>FIG. 2. Leading diagonal of flux covariance matrices (shading) and flux decorrelation times (contours) for (a)–(f) external (atmospheric) and (g)–(i) internal (ocean mesoscale eddy) turbulent fluxes. Contours are separated by half a day and increase in darkness, with thicker, solid contours at 0.5 (lightest), 1.5, and 2.5 (darkest) days. In the latter case, quantities are depth averaged, and contours are separated by 10 days with thicker contours at 10 (lightest), 30, and 50 (darkest) days.</figcaption></figure>

As before, we consider the large-scale climatological cycle to be common to the HRM and the LRM, so we separate (30) into a trajectory (·) and a temporal fluctuation (·'):

$$ \partial_t (\overline{T} + \hat{T}') + \langle \overline{V} + \hat{V}' | \nabla (\overline{T} + \hat{T}') \rangle = -[\overline{\langle V | \nabla T \rangle} + \overline{\langle V | \nabla T' \rangle}], \quad (31) $$

where $\partial_t \overline{T} + \langle \overline{V} | \nabla \overline{T} \rangle = - \overline{\langle V | \nabla T \rangle}$ is the trajectory component common to both LRM and HRM. The unresolved component in the LRM is therefore

$$ \partial_t \hat{T}' + \langle \overline{V} | \nabla \hat{T}' \rangle + \langle \hat{V}' | \nabla \overline{T} \rangle = - \overline{\langle V | \nabla T \rangle}', \quad (32) $$

where smaller, second-order time-fluctuating terms ($\langle \hat{V}' | \nabla \hat{T} \rangle$) have been neglected. The flux terms on the left-hand side represent large-scale interactions between the temporal mean and its fluctuations, while the right-hand side describes the temporal fluctuation of small-scale fluxes. We note that the latter term also contains interactions between the climatology and fluctuations, as can be seen by separating its interior velocity and temperature components into their own time-mean and fluctuating terms. Of those various flux terms, we are interested in the stationary eddy-driven component (i.e., the transport of small-scale buoyancy fluctuations by small-scale current fluctuations). Thus, neglecting large-scale

and seasonal to subseasonal interactions, we are left with the (eddy-driven) internal-flux-induced rate of change in temperature:

$$ F_T' = - \overline{\langle V' | \nabla T' \rangle}', \quad (33) $$

where similar considerations may be made for the internal eddy-driven salt flux ($F_s^l$).

We apply this approach to determine the turbulent eddy heat and salt fluxes unresolved in our LRM (the linear model of section 3a) using an eddy-permitting ocean model (the HRM), which is assumed to resolve a significant portion of these fluxes, despite its own limitations in resolution. In particular, NEMO (v3.5) was run for 20 years in its 1/4°, 75-level ORCA025 configuration, with climatological forcing. The configuration is based on that of Grégorio et al. (2015), who produce the forcing by creating a mean year from the Drakkar Forcing Set (Brodeau et al. 2010). A smoothly forced ocean-only model was chosen to minimize the impact of turbulent atmospheric fluxes (which were determined separately; section 3b). The spatial averaging of (29) was undertaken by averaging all grid points in the HRM which fall within a single grid cell of the LRM.

As in the external case, the time series of resolved internal turbulent fluxes [$F_T^l$ and $F_s^l$, following (33)] were used to determine (via the lag-autocorrelation e-folding time) $\lambda_I$ and (via
<figure><img src="image_4.png" /><figcaption>FIG. 3. Map of error in fitting power spectra of resolved (top) internal (depth-averaged values shown) and (middle), (bottom) external turbulent fluxes to the theoretical power spectrum (Lindner 2009) of (a)–(f) a white Gaussian noise and (g)–(l) an Ornstein–Uhlenbeck process: internal heat flux in (a) and (g), external heat flux in (b) and (h), external momentum flux (zonal component) in (c) and (i), internal salt flux in (d) and (j), and external momentum flux (meridional component) in (f) and (l). The error is given as the ratio of root-mean-square logarithmic error (RMSLE) to the log-mean value of the spectrum (equivalent to the RMSE:mean ratio in log space), such that the whole spectrum is weighted evenly. The error is based solely on the shape of the spectrum, and not the intensity of the fluxes as shown in Fig. 2.</figcaption></figure>
the covariance with other locations) $\Sigma_I$. Owing to the much greater number of elements in $\Sigma_I$ due to the vertical dimension, technical constraints prohibit a fully global treatment of spatial covariance. Instead, we assume spatial covariance to occur only locally: within a $(3 \times 2^\circ)^2 = (6^\circ)^2$ area (i.e., in a nine-point horizontal neighborhood of each location), and throughout the corresponding vertical. Features larger than this would be resolved by the LRM (e.g., Griffies and Treguier 2013). Testing this demonstrated a substantial difference in our results compared with the case with no local covariance, but only a minor difference relative to a larger, 16-point neighborhood. This assumption allows us to use a sparse matrix representation of $\Sigma_I$, reducing computational demand to the same order as that of $\Sigma_E$.

The temperature ($F_I^T$: Fig. 2g) and salinity ($F_S^I$: Fig. 2h) components of the subgrid fluxes can be seen to exhibit generally similar variance distributions, with almost indistinguishable decorrelation time scales (Fig. 2, contours). Common to both components is the strong imprint of the Gulf Stream, Agulhas, Zapiola Gyre, Kuroshio, and Antarctic Circumpolar Current. Their covariance (Fig. 2i) emphasizes these common regions and is effectively everywhere positive, while salinity flux variability uniquely shows strong signatures in the Amazon and Niger outflow regions. There is some latitudinal dependence of decay time (as may be expected from the changing deformation radius; e.g., Chelton et al. 1998) but decay times $\lambda_I$ largely reflect the variance itself, $\Sigma_I$. For example, the shortest times (on the order of days), at the equator, may also be found at much higher latitudes in turbulent regions such as the Gulf Stream. Meanwhile, the gyre interiors show greater persistence, up to many months in the Pacific, and these are the regions where the fluxes are also weakest. These quiescent, persistent regions are understandably where a constant-spectrum approximation (with instantaneous decay) fits most poorly. Consistently this is where the greatest improvements are seen when moving from a Gaussian white noise representation of internal fluxes to an Ornstein–Uhlenbeck process representation (section 2b; Fig. 3), as fitting the error is based on the shape of the spectrum of fluxes, and is not weighted by their intensity.

### d. Experiment design

As described in section 2, we can use our linear model configuration and stochastic approach to analyze the variance evolution of any linear, scalar-valued function of the ocean state, in both a theoretical (optimized stochastic representation) and realistic (diagnosed stochastic representation) context. We choose to focus on a range of climatically relevant metrics: the meridional volume transport (MVT, integrated from the surface to the depth of maximum overturning), full-depth meridional heat transport (MHT), and ocean heat content (OHC; over the present depth range of the majority of the Argo fleet, 0–2000 m). These metrics are calculated for the subtropical [at 25°N for MVT (0–870 m) and MHT (full depth), from 15° to 40°N for OHC (0–2000 m)] and subpolar [at 55°N for MVT (0–1200 m) and MHT (full depth), from 40° to 65°N for OHC (0–2000 m)] North Atlantic. In all

cases, monthly, annually, and decadally averaged quantities are considered.

## 4. Results

### a. Subtropical North Atlantic

#### 1) OPTIMAL STOCHASTIC PERTURBATIONS

We now consider (using the limiting cases of section 2b) the spatially correlated external and decorrelated internal OSPs of the metrics of section 3d in our linearized ocean model (section 3a). The sensitivity of the metric to different potential sources of variability is indicated by the response ratio (Table 1), following (21). For instance, the correlated surface zonal momentum flux OSP of yearly MHT has a response ratio of $0.28 \times 10^{-3}$ PW$^2$ per $(\text{cm s}^{-1})^2 \text{s}^{-1}$ (or $0.28 \times 10^{-3}$ PW$^2 \text{cm}^2 \text{s}^{-3}$). This implies that a stochastic surface zonal momentum flux following the correlated OSP which has a magnitude of $1 (\text{cm s}^{-1})^2 \text{s}^{-1}$ will induce a response in annually averaged MHT with a variance of $0.28 \times 10^{-3}$ PW$^2$ across a large ensemble. This noise amplitude is purely illustrative, as the response ratio is a linear rescaling that is independent of the strength of the noise. The response ratios for MVT and MHT suggest a change in regime when averaging times are increased. For these metrics, sensitivity to large-scale spatially correlated buoyancy fluxes at the surface remains relatively constant at all time scales, producing a response of similar amplitude. Conversely, sensitivity to internal, spatially decorrelated buoyancy fluxes falls sharply with increasing average time, particularly for MVT. Surface momentum flux sensitivity also sees a sharp decline from monthly to annual time scales for both MVT and MHT. OHC variability exhibits no apparent regime shift of this nature, with a steady sensitivity to changes in all variables across all time scales.

So as to understand the mechanisms of variability generation in the model, we now consider the spatial distribution of the perturbations (for year-averaged quantities) in more detail (Figs. 4, 5, and 6). Although the OSP framework does not permit the isolation of individual mechanisms (as it combines every possible mechanism which may impact a desired metric, across all time scales), we may in some cases qualitatively compare clear patterns with prior studies which have highlighted them, particularly studies of adjoint sensitivity fields similar to those which underlie the complex structure of the OSP. The optimal perturbations for MVT and MHT (Figs. 4 and 5; shading) are broadly similar. In the uncorrelated, internal case (panels a and b) the perturbation can have no large-scale structure and simply reflects the sensitivity of the metric in question to each individual point acting alone. These are greatest in the Gulf Stream, and along the evaluation line of the metrics. The large-scale patterns of the correlated external buoyancy forcing, however, reflect strongly the model mean state. In particular, subtropical meridional transport variability displays a strong sensitivity to subpolar surface buoyancy fluxes, reflected as a large-scale gradient across the northern boundary of the subtropical gyre. Wind sensitivity displays very consistent patterns indicating stimulation of Ekman transport (in the case of zonal wind) and western boundary
<table><caption>TABLE 1. Normalized response ratios of various OSPs [following (21)] for the North Atlantic subtropical region. The unit of the response ratio is that of output variance (given in the row header) divided by that of the normalized noise amplitude (given in the column header). A stochastic forcing with the spatial distribution of the OSP (Figs. 4, 5, and 6) and with unit amplitude will thus stimulate a response in the target metric with the given response variance. Note that, in line with (26), the response ratio is independent of the amplitude of the noise in the linear framework. The amplitudes shown in the table have unit amplitude for illustrative purposes. Left-hand columns correspond to perfectly correlated surface OSPs. Right-hand columns correspond to totally uncorrelated internal OSPs. Note that the amplitude units differ between the correlated and uncorrelated cases.</caption><thead><tr><th rowspan="2"></th><th rowspan="2">Avg time</th><th colspan="4">Surface (correlated)</th><th colspan="2">Full-depth (uncorrelated)</th></tr><tr><th>T (×10<sup>-3</sup> K<sup>2</sup> s<sup>-1</sup>)</th><th>S (×10<sup>-4</sup> psu<sup>2</sup> s<sup>-1</sup>)</th><th>u [(cm s<sup>-1</sup>)<sup>2</sup> s<sup>-1</sup>]</th><th>v [(cm s<sup>-1</sup>)<sup>2</sup> s<sup>-1</sup>]</th><th>T (×10<sup>-3</sup> K<sup>2</sup> day<sup>-1</sup>)</th><th>S (×10<sup>-4</sup> psu<sup>2</sup> day<sup>-1</sup>)</th></tr></thead><tbody><tr><td rowspan="3">MVT (Sv<sup>2</sup>)</td><td>30 days</td><td>0.23</td><td>0.60</td><td>0.14</td><td>0.08</td><td>21.74</td><td>53.86</td></tr><tr><td>1 yr</td><td>0.22</td><td>0.56</td><td>0.05</td><td>0.04</td><td>12.02</td><td>30.43</td></tr><tr><td>10 yr</td><td>0.19</td><td>0.50</td><td>0.04</td><td>0.03</td><td>2.94</td><td>7.81</td></tr><tr><td rowspan="3">MHT (×10<sup>-3</sup> PW<sup>2</sup>)</td><td>30 days</td><td>1.16</td><td>2.73</td><td>1.88</td><td>0.52</td><td>82.88</td><td>201.71</td></tr><tr><td>1 yr</td><td>1.10</td><td>2.62</td><td>0.28</td><td>0.23</td><td>53.95</td><td>134.33</td></tr><tr><td>10 yr</td><td>0.89</td><td>2.16</td><td>0.20</td><td>0.15</td><td>13.65</td><td>37.13</td></tr><tr><td rowspan="3">OHC (×10<sup>-3</sup> K<sup>2</sup>)</td><td>30 days</td><td>2.74</td><td>6.60</td><td>0.19</td><td>0.23</td><td>12.89</td><td>42.24</td></tr><tr><td>1 yr</td><td>2.77</td><td>6.57</td><td>0.19</td><td>0.22</td><td>12.68</td><td>41.63</td></tr><tr><td>10 yr</td><td>2.55</td><td>6.22</td><td>0.18</td><td>0.21</td><td>10.38</td><td>35.56</td></tr></tbody></table>

transport change combined with eastern boundary up- or downwelling (in the case of meridional wind). Upwelling directly impacts the volume transport locally through geostrophy (Hirschi et al. 2007; Kanzow et al. 2010; Polo et al. 2014), but this pattern has also been observed in other sensitivity studies to trigger pressure anomalies which reach great distances along the eastern boundary (Pillar et al. 2016; Jones et al. 2018).

The optimal OHC perturbation in the uncorrelated case (Figs. 6a,b) shows sensitivity to buoyancy fluxes throughout the region, but particularly at the subpolar-subtropical gyre interface, which has been highlighted as a key region for variability generation in the Atlantic (Buckley and Marshall 2016). Also clear, but less pronounced, are local peaks around the Agulhas retroflection and the Zapiola Gyre. The correlated surface OSPs are notably different in the cases of temperature and salinity due to the ability of surface temperature fluxes to impact heat content variability both through passive and active mechanisms. These mechanisms, respectively inducing heat content anomalies through anomalous heat fluxes (passive) and through dynamical density changes that redistribute the existing heat reservoir (active), can sometimes conflict (Stephenson and Sévellec 2021). The active mechanisms are made clear by the correlated salinity OSP, which cannot stimulate OHC passively and which instead stimulates clear positive/negative density gradients across the northern boundary of the North Atlantic and South Atlantic subtropical gyres, as well as a local peak in the deep water formation region of the model (Stephenson et al. 2020). The temperature perturbation echoes this, but with a distribution which is almost everywhere equally signed, so as to passively stimulate heat content. The momentum flux perturbations (Figs. 6e,f) are generally more complex but can still be seen to broadly coincide with predominantly zonal streamlines and coastal regions in the zonal and meridional cases, respectively. There is a notable focus along the subpolar-subtropical gyre interface for the zonal momentum flux.

## 2) DYNAMICAL ATTRIBUTION OF SUBTROPICAL VARIANCE

Having explored the theoretical patterns and mechanisms by which oceanic variability can be optimally stimulated in our model, we turn our attention to the ways in which it is stimulated in a model representation of the real climate system, as derived in section 2c. Following (26) and (27), application of each component of the stochastic forcing separately allows the resultant variance evolution to be partitioned accordingly (Fig. 7). There is a substantial difference between the nature of month- and decade-averaged transport metrics, both in the variance amplitude and in the impacts of different sources, as in the OSP case (shown by the response ratios of Table 1). External momentum fluxes are responsible for 52% of month-averaged MVT and for 63% of month-averaged MHT by the end of the 60-yr simulation, but just 9% and 10%, respectively, for decade-averaged MHT. Similarly, the external buoyancy component contributes just 4% to month-averaged MVT variance at 60 years, but over 50% in the decade-averaged case. For year-averaged MVT and MHT, the ocean internal component is the dominant contributor to the final variance, at
<figure><img src="image_5.png" /><figcaption>FIG. 4. Optimal stochastic perturbation for year-averaged subtropical meridional volume transport in (a),(b) the fully spatially uncorrelated cases (depth averages shown) and (c)–(f) the perfectly spatially correlated surface-only cases. Streamlines show time-averaged volume transport over the upper 2 km in the trajectory. Dashed lines show 25°N, the latitude at which the meridional volume transport is evaluated.</figcaption></figure>

60% for MVT and 58% for MHT. In addition to differences between monthly and decadal metrics in the final (60 yr) variability, a difference in the evolution of this variance is also apparent. Contributions from all sources are fairly steady in time for MVT and MHT for the quickly-saturating month-averaged case. For 10-yr average MVT and MHT, there is a more notable shift. Following initialization, external momentum and internal buoyancy fluxes are the main causes of error growth. However, the contribution of wind peaks abruptly, while the eddy component grows for around 6 years, peaking at nearly 80% of the total uncertainty. On longer time scales, the mesoscale-eddy-driven component falls to slightly less than half of the total contribution over the remainder of the simulation. During this stage, it is the more slowly acting external buoyancy component that develops and contributes the remaining variance.

Notably, as in the OSP perspective, the components of the OHC variance after 60 years are consistent across different time averages, with an almost equal contribution (around 45% each) from external and internal buoyancy fluxes. This follows the slow growth of the internal component, which, at its lowest, contributes only around 25% of the total uncertainty. This is in contrast with the MVT and MHT, where it is the external buoyancy contribution that is the slowest to develop.

Following (28), we consider the spatial distributions of these contributions to the 60-yr variance for the annually averaged case, within the transition between the two discussed (month- and decade-average) cases (Fig. 8; where the zonal and meridional momentum flux covariance contribution is not shown). There is generally a high level of agreement between the patterns shown in the optimal case
<figure><img src="image_6.png" /><figcaption>FIG. 5. As in Fig. 4, but for subtropical meridional heat transport (evaluated at 25°N, denoted by the dashed line).</figcaption></figure>

(i.e., what the ocean “wants”; Figs. 4, 5, and 6) and the realistic case (Fig. 8). This is linked to the overall relative constant shape of the diagnosed, prescribed forcing (i.e., what the ocean “gets”; Fig. 2). Although we remind of the contrast between the two frameworks (i.e., white versus temporally correlated noise) when making any such comparisons.

In particular, volume and heat transport variability are primarily driven by ocean internal buoyancy fluxes local to the western boundary, and by remote external buoyancy fluxes in the subpolar region. Zonal surface momentum fluxes, consistently with the OSP, almost exclusively stimulate a zonal band along the evaluation line (Figs. 8g,h), while in the meridional case a combination of western boundary current and eastern along-shelf stimulation pervade. The agreement between the prescribed (temporally correlated) and optimal (white noise) forcing is less apparent in the case of OHC. Internal buoyancy fluxes affecting heat content variability can be predominantly traced in the prescribed case to highly focused sources in the

noisiest regions of the Atlantic (Fig. 8c vs Figs. 2g,h), while the optimal white noise perturbation is more evenly distributed throughout the Atlantic with a local peak in the subtropical-subpolar “transition zone” (Buckley and Marshall 2016). The distribution in the prescribed case also exhibits a selection of locations which make a negative contribution, particularly north of the North Atlantic Current. These arise from the covariance of neighboring points with an otherwise strong contribution gradient, and act as a compensatory “source” of predictability relative to that which would stem from a spatially decorrelated representation. External buoyancy fluxes contribute over a broader area than the internal case, with the most concentrated contributions in the remote subpolar region. The contribution from zonal wind is almost exclusively along the evaluation region’s boundaries, whereas in the meridional case (as also seen in the OSP) the western coasts of Europe and South America have the clearest impact.
<figure><img src="image_7.png" /><figcaption>FIG. 6. As in Fig. 4, but for subtropical ocean heat content (evaluated between 15° and 40°N, denoted by the two dashed lines).</figcaption></figure>

b. *Subpolar North Atlantic*

1) OPTIMAL STOCHASTIC PERTURBATIONS

Applying the considerations of section 4a(1) to the sub-
polar region, differences emerge in the amplitude of the
response to the optimal stochastic forcing (Table 2). For
subpolar MVT, the correlated surface OSP is much more
effective at generating variability than in the subtropics,

particularly on annual time scales (for which the response
ratio is around 4 times as large in the subpolar region). For
MHT, the values are similar in both regions. The opposite is
apparent in the spatially uncorrelated case, where for ex-
ample the response of monthly MHT to its uncorrelated
optimal noise perturbation is over 6 times as large in the
subtropics as in the subpolar region. OHC again shows
consistent behavior across all time averages, but is much
<figure><img src="image_8.png" /></figure>

FIG. 7. Attribution of uncertainty following initialization for the subtropical ocean metrics [(a)–(c) MVT, (d)–(f) MHT, and (g)–(i) OHC] over different averaging times [(top) month, (middle), year, and (bottom) decade], following (26) and (27). Green and red shading indicate variance due to external (atmospheric) momentum and buoyancy fluxes, respectively. Blue shading indicates variance due to internal buoyancy fluxes (due to oceanic mesoscale eddy forcing). Dashed white contours show percentages (inset text) of the total variance. Shaded boxes show the averaging window over which the metric is evaluated. Variance due to surface momentum fluxes is partitioned into zonal (dark green) and meridional (light green) components, where shading between them indicates covariance.

more sensitive to external momentum and internal buoyancy changes than in the subtropics.

The OSP for meridional volume transport (Fig. 9) shows a much more concentrated spatial distribution than its subtropical equivalent. In the uncorrelated ocean interior case, almost all of the weight is focused at the core of the subpolar gyre (Figs. 9a,b). For the perfectly correlated surface case, this hotspot, coincident with the surface outcrop of the model North Atlantic Deep Water (Stephenson et al. 2020), is complemented by a dipole pattern crossing the North Atlantic Current (Figs. 9c,d). This dipole resembles the surface sensitivity of the least damped interdecadal mode of variability (corresponding to a large-scale thermal Rossby wave) present in an earlier version of the model (Sévellec and Fedorov 2013). As with the subtropical metric, the optimal momentum flux patterns are an east-west band in the zonal case and a predominantly eastern-boundary-following pattern in the meridional case.

While having many common features with that of MVT, the optimal pattern for MHT (Fig. 10) is much less focused, neglecting the hotspot of the northwest Atlantic for a more spread out distribution. The optimal internal perturbation consists of buoyancy fluxes throughout the subpolar gyre, as

well as in the subtropical-subpolar intergyre region. In the correlated case, the dipole feature between gyres (also visible for the subtropical case; Fig. 5) is more heavily emphasized. In addition to the familiar features of the velocity OSPs, fainter bands encircle the subpolar gyre.

The OSPs of subpolar OHC variance (Fig. 11) exhibit many similar behaviors to those described for other metrics. The uncorrelated interior noise favors the subtropical-subpolar gyre boundary, while the correlated surface heat flux pattern targets oppositely the deep water outcrop regions and the wider North Atlantic, with a particular focus on the North Atlantic Current. Similarly to the correlated OSP of subtropical heat content, the correlated subpolar zonal velocity OSP displays a complex arrangement of alternating bands which broadly coincide with strongly zonal currents in the trajectory, while the meridional pattern predominantly targets coastal upwelling and downwelling (i.e., alongshore velocity/momentum fluxes) in these same regions.

## 2) DYNAMICAL ATTRIBUTION OF SUBPOLAR VARIANCE

Under prescribed sources of variability, the subpolar region is dominated by external forcing (Fig. 12), which accounts for
<figure><img src="image_9.png" /><figcaption>FIG. 8. Spatial distribution of sources of accumulated variance for subtropical ocean metrics [(top) MVT, (middle) MHT, and (bottom) OHC] after 60 years of simulation following (26), (27), and (28). Variance per unit volume due to internal buoyancy fluxes is depth integrated to give the (a)–(c) water column total contribution per unit area while (d)–(l) the variance due to external momentum [zonal component in (g)–(i), meridional component in (j)–(l)], and buoyancy fluxes in (d)–(f) are surface distributions of contribution per unit area. Dashed lines show the latitude (MVT and MHT metrics) or region (OHC metric) where the metric is evaluated. Note the differing (sometimes by orders of magnitude) color scales, reflecting the differing contributions shown in Fig. 7.</figcaption></figure>

up to 94% of the total variance after 60 years in the case of month-averaged heat transport. As in the subtropics, the meridional transport metrics exhibit a regime shift when moving from month-averaged quantities (up to 86% momentum-driven) to decade-averaged quantities (where over 60% of the final variance can be attributed to surface buoyancy fluxes). For all time averages, momentum fluxes contribute most of the early-stage error growth of MVT and MHT following initialization, but the buoyancy component becomes more established over

the first decade. MVT and MHT are much less variable overall than in the subtropics, while heat content variance is slightly higher, again showing consistent behavior across all considered time averages. Also notable is that, despite full convergence not being reached after the 60 years, heat content seemingly shows a higher degree of saturation in the subpolar region than in the subtropical region.

The spatial patterns of subpolar variance origins in response to prescribed fluxes (Fig. 13; where the zonal and meridional
<table><caption>TABLE 2. As in Table 1, but for subpolar OSPs, (whose spatial distributions are shown in Figs. 9, 10, and 11).</caption><thead><tr><th rowspan="2"></th><th rowspan="2">Avg time</th><th colspan="4">Surface (correlated)</th><th colspan="2">Full-depth (uncorrelated)</th></tr><tr><th>T (×10<sup>-3</sup> K<sup>2</sup> s<sup>-1</sup>)</th><th>S (×10<sup>-4</sup> psu<sup>2</sup> s<sup>-1</sup>)</th><th>u [(cm s<sup>-1</sup>)<sup>2</sup> s<sup>-1</sup>]</th><th>v [(cm s<sup>-1</sup>)<sup>2</sup> s<sup>-1</sup>]</th><th>T (×10<sup>-3</sup> K<sup>2</sup> day<sup>-1</sup>)</th><th>S (×10<sup>-4</sup> psu<sup>2</sup> day<sup>-1</sup>)</th></tr></thead><tbody><tr><td rowspan="3">MVT (Sv<sup>2</sup>)</td><td>30 days</td><td>0.39</td><td>1.36</td><td>0.26</td><td>0.18</td><td>7.88</td><td>27.60</td></tr><tr><td>1 yr</td><td>0.85</td><td>2.09</td><td>0.53</td><td>0.20</td><td>3.79</td><td>12.72</td></tr><tr><td>10 yr</td><td>0.30</td><td>1.07</td><td>0.18</td><td>0.12</td><td>0.92</td><td>2.96</td></tr><tr><td rowspan="3">MHT (×10<sup>-3</sup> PW<sup>2</sup>)</td><td>30 days</td><td>1.71</td><td>4.27</td><td>0.47</td><td>0.38</td><td>13.74</td><td>36.14</td></tr><tr><td>1 yr</td><td>1.54</td><td>4.03</td><td>0.41</td><td>0.34</td><td>11.24</td><td>29.84</td></tr><tr><td>10 yr</td><td>0.83</td><td>2.25</td><td>0.19</td><td>0.15</td><td>3.69</td><td>10.45</td></tr><tr><td rowspan="3">OHC (×10<sup>-3</sup> K<sup>2</sup>)</td><td>30 days</td><td>4.76</td><td>12.35</td><td>1.99</td><td>1.46</td><td>48.31</td><td>147.22</td></tr><tr><td>1 yr</td><td>4.71</td><td>12.30</td><td>1.98</td><td>1.46</td><td>47.85</td><td>145.97</td></tr><tr><td>10 yr</td><td>4.33</td><td>11.29</td><td>1.87</td><td>1.39</td><td>38.63</td><td>121.61</td></tr></tbody></table>

momentum flux covariance contribution is not shown) are generally less similar to the corresponding optimal perturbations [section 4b(1)] than in the subtropics (section 4a), although we again treat comparisons between the two frameworks with caution. We further note that the prescribed forcing may project more strongly onto nonleading eigenvectors, but that we restrict our comparisons to the leading eigenvector, as it corresponds to the optimal forcing. The differences are particularly clear for internal buoyancy fluxes, which for all metrics share a common maximum at around 40°N, far south of the corresponding peaks in the uncorrelated OSPs. For MVT there is a large contribution on the evaluation line west of Scotland, apparently coincident with a local peak in the uncorrelated OSP, but the most sensitive region in the central subpolar gyre is only weakly stimulated. As in the subtropical region, negative contributions flank the Gulf Stream and its extension, acting as a compensatory “source” of predictability offsetting its covarying sinks. Variance due to (temporarily correlated) prescribed external buoyancy fluxes more closely agrees with the (white noise) spatially correlated OSP. In particular, the northern portion of the optimal dipole shape is discernible for MHT, while the deep water outcrop hotspot can be faintly recognized, along with the west-European shelf in the case of MVT. Heat content variability due to external buoyancy fluxes largely coincides with the most concentrated region of the correlated OSP, in the North Atlantic Current, but shows little agreement elsewhere. The external momentum flux components are qualitatively similar for all three metrics, again stimulating transport across constant latitude lines in the zonal case [where the noise input (Fig. 2) constructively stimulates the most sensitive regions (Figs. 9–11)], while highlighting the coasts for the meridional case. Both zonal (in the case of MVT and OHC) and meridional (in the case of MHT) momentum flux contributions are offset by a weakly negative compensation bordering the regions of strongest positive variance stimulation.

## 5. Discussion and conclusions

The climate system contains a number of sinks of predictability or, equivalently, sources of uncertainty, from which unpredictable noise can grow and eventually overwhelm predictable signal (such as that provided to an initialized forecast). In this study, we have considered the sources (and compensatory sinks) of uncertainty in metrics of the North Atlantic from two perspectives. In the first perspective, two distinct optimal stochastic perturbations were calculated, encapsulating the patterns which generate maximum variance in the metric. These are a representation of the sensitivity of the metric to random forcing. The perturbations differ by their spatial coherence: one being fully spatially uncorrelated, with the other fully correlated over the surface layer. These are the extrema of possible spatial correlation, and respectively mimic, in an idealized sense, the behavior of stochastic fluxes due to (mesoscale) oceanic turbulence and (synoptic scale) atmospheric turbulence. In the second perspective, the optimal stochastic forcing was instead replaced with a prescribed, model-diagnosed stochastic representation of
<figure><img src="image_10.png" /><figcaption>FIG. 9. As in Fig. 4, but for subpolar meridional volume transport (evaluated at 55°N denoted by the dashed line).</figcaption></figure>

these sources, including spatiotemporal covariance. The
properties of the representations were diagnosed from more
complex (fully coupled and eddy-permitting) models. This
has allowed us to compare the commonalities between the
optimal and representative cases (albeit in a limited way,

given their differences in spatiotemporal correlation). We
have further been able, in the diagnosed case, to dynami-
cally attribute variability to its origins. The latter ability
notably forgoes the more typical ensemble attribution ap-
proach, which generally necessitates many simulations in a

<figure><img src="image_11.png" /><figcaption>FIG. 10. As in Fig. 4, but for subpolar meridional heat transport (evaluated at 55°N denoted by the dashed line).</figcaption></figure>
OHC (SUBPOLAR)

<figure><img src="image_12.png" /><figcaption>FIG. 11. As in Fig. 4, but for subpolar heat content (evaluated between 40°N and 65°N denoted by the two dashed lines).</figcaption></figure>

high-complexity model, and cannot ensure causality. The sources determined by these two perspectives can qualitatively be seen as what the ocean “wants” (in order to maximize variability) and what the ocean “gets” (in the real world). Regions where the ocean “gets what it wants” offer particularly poor prospects for prediction, as both the sources of uncertainty and their mechanisms of transformation into large-scale uncertainty play a role. In such locations, any signal provided to an initialized model is more likely to degrade due

to the coaction of noise and the ocean's dynamical response to it. This corresponds in our results to faster growth rates and greater uncertainties.

Variations on the OSP technique have been utilized in the context of optimal excitation of MOC variability in a number of studies [a thorough review is provided by Monahan et al. (2008)]. However, due to the complexity of the problem these studies are typically undertaken in an idealized context, utilizing either box models (e.g., Tziperman and Ioannou 2002;
<figure><img src="image_13.png" /><figcaption>FIG. 12. As in Fig. 7, but for subpolar ocean metrics</figcaption></figure>

Zanna and Tziperman 2008) or idealized ocean models (e.g., Sévellec et al. 2007, 2009). We have adapted the framework to a global OGCM by reducing the covariance matrix to block diagonal form and considering its limiting cases. We note (e.g., Farrell and Ioannou 1996) the close relationship between optimal stochastic forcings and optimal initial perturbations: the former is in a sense a linear combination of the latter such that the coefficients are determined by the OSP approach. As the linear optimal perturbation of a linear ocean metric is simply a rescaling of the adjoint sensitivity field (Sévellec et al. 2007), we may consider the sources highlighted by the OSP in the context of past adjoint sensitivity studies, where they appear robust across differing models, metrics, and time scales. Recurring mechanisms evident in our study include, for instance, the along-shelf stimulation by meridional wind and subsequent triggering of coastal pressure anomalies, particularly along the west coast of Africa. This pattern has been stressed by Jones et al. (2018) in an adjoint sensitivity study of Labrador Sea heat content, Loose et al. (2020) regarding heat transport across the Greenland–Scotland Ridge, and Pillar et al. (2016) in the context of meridional overturning in the subtropics. The latter study additionally analyzes fainter alternating bands of wind stress sensitivity as also seen here, concluding that these communicate pressure anomalies via topographically steered Rossby waves.

Common to the surface thermohaline OSPs of all metrics considered here is a large-scale buoyancy gradient pivoting on the North Atlantic Current, which has in dynamical studies been seen to stimulate subtropical (Pillar et al. 2016; Kostov et al. 2019) and subpolar (Sévellec et al. 2017) volume transport, as well as basinwide (Sévellec and Fedorov 2017) and Labrador Sea (Jones et al. 2018) heat content. This is joined by a “hotspot” common to the heat content and subpolar volume transport OSPs in both the correlated and uncorrelated cases which is associated with the passive transport of buoyancy anomalies via deep water pathways (Sévellec and Fedorov 2015; Stephenson et al. 2020).

To estimate the extent to which these intrinsic ocean sensitivities are exploited by actual sources of stochastic variability, and to quantify the respective contribution of these sources to oceanic uncertainty, we then considered the metrics from the second, prescriptive, perspective. A number of studies have dynamically attributed oceanic changes to prescribed external surface forcings using adjoint methods (Pillar et al. 2016; Sévellec et al. 2018; Smith and Heimbach 2019) but the relative quantification of internal oceanic mesoscale eddy contributions has thus far been restricted to a resource-intensive ensemble framework (e.g., Bessières et al. 2017). The associated sensitivity to initial conditions suggests that these contributions may present a key sink of predictive skill in initialized high-resolution climate models however, and so are of increasing
<figure><img src="image_14.png" /><figcaption>FIG. 13. As in Fig. 8, but for subpolar ocean metrics. We note again that the differing contributions (as shown in Fig. 12) lead to large differences in the color scales between panels.</figcaption></figure>

importance. By incorporating temporal correlation, we have presented a realistic stochastic representation (an Ornstein-Uhlenbeck process) of the slowly evolving ocean mesoscale that can also be projected onto the adjoint sensitivity fields. This stochastic representation fits the power spectrum of modeled eddy buoyancy fluxes much more closely than Gaussian white noise, which is the more commonly employed framework when considering atmospherically driven low-frequency variability [e.g., the review of Farneti (2017), and references therein]. This has allowed us to bypass the ensemble approach in exchange for the more numerically efficient dynamical method for both oceanic (internal) and atmospheric (external) sources of error growth.

The diagnosed stochastic forcing approach reveals a regime change in meridional transport variability for longer time averages. In particular, we have shown that surface momentum fluxes dominate for month-averaged transport metrics while surface buoyancy fluxes take over for decade averages. This regime shift is well documented (Dong and Sutton 2001; Hirschi et al. 2007; Polo et al. 2014) but we find that in the early stages of the error growth, and for annual averages, it is ocean internal buoyancy fluxes, due to mesoscale eddies, that form the greatest contribution in the subtropics. As early-stage growth is when the signal-to-noise ratio diminishes most rapidly, it may be internal sources that present the greatest barrier to subtropical predictability. Our results indicate that these

sources ultimately account for up to 60% of annually averaged volume transport variability at 25°N. This quantification broadly agrees with the varying estimates of ensemble studies (albeit at the higher end; e.g., Grégorio et al. 2015; Jamet et al. 2019), which typically place a local peak in internal oceanic contributions to MVT variability near 25°N (our subtropical metric latitude) with a corresponding trough near 55°N (our subpolar metric latitude) consistently with the decrease we show here. This coincides with an overall lower maximum amplitude of meridional transport variability at 55°N, which is reached much more slowly, suggesting that observations near this latitude (e.g., from the OSNAP observational array; Lozier et al. 2017) will remain predictable for longer than those near 25°N (e.g., from the RAPID observational array; Cunningham et al. 2007), particularly for lower frequencies. We did not find any such regime shift in the case of ocean heat content, whose variability for all time averages is dominated by external forcing, particularly in the more quiescent subpolar region [consistent with the ensemble study of Sérazin et al. (2017)]. Conversely to MVT and MHT, the slower growth in OHC uncertainty in the subtropics suggests it is more predictable there than in the subpolar region, where greater depths are frequently exposed to atmospheric forcing.

When comparing the theoretically deduced (white noise) OSPs with the sources of variability in response to diagnosed (temporally correlated) stochastic forcing, a general overlap
was observed in the subtropical region. This suggests efficient stimulation of the preferred mechanisms of the ocean, despite the differing temporal correlation of the two frameworks. This was less true of the subpolar region, which may go some way to explaining the smaller diagnosed variance there relative to the subtropics, despite its higher sensitivity to surface forcing (quantified via response ratios) in the optimal framework. Regarding the subsurface component, it is commonly discussed that the smaller deformation radius at higher latitudes necessitates an ocean model with a fully eddy-resolving resolution in order to faithfully represent the internal contribution. As such, this contribution is likely underrepresented in eddy-permitting ensemble studies, which typically portray it as very minor (e.g., Grégorio et al. 2015; Leroux et al. 2018). This also impacts our own approach of diagnosing mesoscale eddy fluxes in an eddy-permitting model. However, we reinforce that even without prescribed forcing, the theoretical OSP framework has allowed us to quantify the subtropical sensitivity of meridional transports to spatially uncorrelated noise as being many times as large as those in the subpolar region. It is thus apparent that these large-scale oceanic metrics are simply less affected by small-scale noise in this region, potentially offering increased benefit from targeted monitoring systems, which can be used to initialize forecasts in such a way that the initialized signal is better protected from degradation by chaotic internal noise. Notably, the opposite was found to be true for subpolar heat content, which also displayed faster growing uncertainty in the diagnosed case.

Previous studies investigating interactions between the oceanic mesoscale and the low-frequency large scale (such as the interactions considered here) present conflicting behavior. While some studies show constructive stimulation of low-frequency variability (e.g., Berloff et al. 2007; Arbic et al. 2014), others show its destruction by small-scale noise (e.g., LaCasce and Pedlosky 2004; Hochet et al. 2020; Sévellec et al. 2021). Other approaches still suggest that ocean eddies may increase predictability by facilitating surface–depth communication, increasing persistence (e.g., Kravtsov 2020; Zhang et al. 2021). The framework of our study describes variability from a linear, ensemble perspective in which any divergence in phase space constitutes an irreversible accumulation of error (a source of uncertainty). This framework is not well suited to isolating such destructive feedbacks, but we have seen that some contributors to the net positive error growth are weakly negative. This slows this growth and restores some predictability. This is particularly apparent along the boundaries of noisy regions such as the North Atlantic Current.

We finally comment on some other limitations of the approach. While computationally efficient, we have used a linearized model under the assumption of small deviations from a trajectory, alongside a stationary, band-limited stochastic representation of dynamical processes which, in reality, are highly intricate. For example, our internal turbulent buoyancy flux representation cannot encompass coherent interbasin exchanges, which have been speculated to be an important mechanism of Atlantic MOC variability (e.g., Biastoch et al. 2008). While a coupled climate model was used to determine

the surface fluxes, the modeled ocean response is unable to interact with these, precluding the existence of any coupled feedbacks and associated modes, which may have a pronounced impact on interdecadal variability (e.g., Liu 2012). Despite these drawbacks, the analytical foundation of the framework leads to an efficient and thorough method for investigating the sources of oceanic variance and associated impacts on predictability. The result is a causal attribution of oceanic uncertainty (otherwise requiring a theoretically infinite ensemble) that can be cleanly partitioned into its sources and locations.

*Acknowledgments.* This research was supported by the U.K. Natural Environmental Research Council (SMURPHS, NE/N005767/1 and MESO-CLIP, NE/K005928/1 and SPITFIRE, NE/L002531/1) and by the DECLIC and Meso-Var-Clim projects funded through the French CNRS/INSU/LEFE program. The authors kindly thank Victor Estella Perez and Juliette Mignot for performing the coupled model simulation and Thierry Penduff and Jean-Marc Molines for supplying source code and forcing files for the eddying model simulation, as well as Simon Müller for his assistance with the adjoint model configuration.

*Data availability statement.* Data used to produce this work are available as follows. Output used from the IPSL-CM5A model is located at https://doi.org/10.5281/zenodo.4300471. Source code, configuration files and climatological forcing files for NEMO v3.5 in the ORCA025 configuration are available at https://doi.org/10.5281/zenodo.4473198. The source code for NEMO/NEMOTAM v3.4 in the ORCA2 configuration is available from the NEMO team at https://forge.ipsl.jussieu.fr/nemo/svn/NEMO/releases/release-3.4, with configuration and normal year forcing files for ORCA2-LIM located at https://doi.org/10.5281/zenodo.1471702. Modifications to model source code specific to our experiments, scripts to run the experiments themselves, and diagnostic scripts used here can be found at https://doi.org/10.5281/zenodo.5014706.

## REFERENCES

Arbic, B. K., M. Müller, J. G. Richman, J. F. Shriver, A. J. Morten, R. B. Scott, G. Sérazin, and T. Penduff, 2014: Geostrophic turbulence in the frequency–wavenumber domain: Eddy-driven low-frequency variability. *J. Phys. Oceanogr.*, **44**, 2050–2069, https://doi.org/10.1175/JPO-D-13-054.1.

Berloff, P., A. M. C. Hogg, and W. Dewar, 2007: The turbulent oscillator: A mechanism of low-frequency variability of the wind-driven ocean gyres. *J. Phys. Oceanogr.*, **37**, 2363–2386, https://doi.org/10.1175/JPO3118.1.

Bessières, L., and Coauthors, 2017: Development of a probabilistic ocean modelling system based on NEMO 3.5: Application at eddying resolution. *Geosci. Model Dev.*, **10**, 1091–1106, https://doi.org/10.5194/gmd-10-1091-2017.

Biastoch, A., C. W. Böning, and J. Lutjeharms, 2008: Agulhas leakage dynamics affects decadal variability in Atlantic overturning circulation. *Nature*, **456**, 489–492, https://doi.org/10.1038/nature07426.

Brodeau, L., B. Barnier, A.-M. Treguier, T. Penduff, and S. Gulev, 2010: An ERA40-based atmospheric forcing for global ocean
circulation models. Ocean Modell., 31, 88-104, https://doi.org/
10.1016/j.ocemod.2009.10.005.

Buckley, M. W., and J. Marshall, 2016: Observations, inferences, and mechanisms of the Atlantic meridional overturning circulation: A review. *Rev. Geophys.*, **54**, 5-63, https://doi.org/10.1002/2015RG000493.

Chang, P., R. Saravanan, T. DelSole, and F. Wang, 2004: Predictability of linear coupled systems. Part I: Theoretical analyses. *J. Climate*, **17**, 1474–1486, https://doi.org/10.1175/1520-0442(2004)017<1474:POLCSP>2.0.CO;2.

Chelton, D. B., R. A. deSzoeke, M. G. Schlax, K. El Naggar, and N. Siwertz, 1998: Geographical variability of the first baroclinic Rossby radius of deformation. *J. Phys. Oceanogr.*, **28**, 433–460, https://doi.org/10.1175/1520-0485(1998)028<0433:GVOTFB>2.0.CO;2.

—, M. G. Schlax, R. M. Samelson, and R. A. de Szoeke, 2007: Global observations of large oceanic eddies. *Geophys. Res. Lett.*, **34**, L15606, https://doi.org/10.1029/2007GL030812.

Collins, M., and B. Sinha, 2003: Predictability of decadal variations in the thermohaline circulation and climate. *Geophys. Res. Lett.*, **30**, 1306, https://doi.org/10.1029/2002GL016504.

Cunningham, S. A., and Coauthors, 2007: Temporal variability of the Atlantic meridional overturning circulation at 26.5°N. *Science*, **317**, 935–938, https://doi.org/10.1126/science.1141304.

Dirac, P. A. M., 1939: A new notation for quantum mechanics. *Mathematical Proceedings of the Cambridge Philosophical Society*, Cambridge University Press, Vol. 35, 416–418.

Dong, B.-W., and R. Sutton, 2001: The dominant mechanisms of variability in Atlantic ocean heat transport in a coupled ocean-atmosphere GCM. *Geophys. Res. Lett.*, **28**, 2445–2448, https://doi.org/10.1029/2000GL012531.

Duan, J., and W. Wang, 2014: Stochastic calculus in Hilbert space. *Effective Dynamics of Stochastic Partial Differential Equations*, Elsevier, 21–45.

Dufresne, J.-L., and Coauthors, 2013: Climate change projections using the IPSL-CM5 Earth System Model: From CMIP3 to CMIP5. *Climate Dyn.*, **40**, 2123–2165, https://doi.org/10.1007/s00382-012-1636-1.

Errico, R. M., 1997: What is an adjoint model? *Bull. Amer. Meteor. Soc.*, **78**, 2577–2592, https://doi.org/10.1175/1520-0477(1997)078<2577:WIAAM>2.0.CO;2.

Farneti, R., 2017: Modelling interdecadal climate variability and the role of the ocean. *Wiley Interdiscip. Rev.: Climate Change*, **8**, e441, https://doi.org/10.1002/wcc.441.

Farrell, B. F., and P. J. Ioannou, 1996: Generalized stability theory. Part I: Autonomous operators. *J. Atmos. Sci.*, **53**, 2025–2040, https://doi.org/10.1175/1520-0469(1996)053<2025:GSTPIA>2.0.CO;2.

Frankignoul, C., and K. Hasselmann, 1977: Stochastic climate models. Part II: Application to sea-surface temperature anomalies and thermocline variability. *Tellus*, **29**, 289–305, https://doi.org/10.3402/tellusa.v29i4.11362.

Grégorio, S., T. Penduff, G. Sérazin, J.-M. Molines, B. Barnier, and J. Hirschi, 2015: Intrinsic variability of the Atlantic meridional overturning circulation at interannual-to-multidecadal time scales. *J. Phys. Oceanogr.*, **45**, 1929–1946, https://doi.org/10.1175/JPO-D-14-0163.1.

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

—, and A. M. Treguier, 2013: *Ocean Circulation Models and Modeling*. International Geophysics Series, Vol. 103, Elsevier, 521–551.

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

Gutjahr, O., D. Putrasahan, K. Lohmann, J. H. Jungclaus, J.-S. Storch, N. Brüggemann, H. Haak, and A. Stössel, 2019: Max Planck Institute Earth System Model (MPI-ESM1.2) for the High-Resolution Model Intercomparison Project (HighResMIP). *Geosci. Model Dev.*, **12**, 3241–3281, https://doi.org/10.5194/gmd-12-3241-2019.

Haarsma, R., and Coauthors, 2020: HighResMIP versions of EC-Earth: EC-Earth3P and EC-Earth3P-HR—Description, model computational performance and basic validation. *Geosci. Model Dev.*, **13**, 3507–3527, https://doi.org/10.5194/gmd-13-3507-2020.

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

Held, I., and Coauthors, 2019: Structure and performance of GFDL’s CM4.0 climate model. *J. Adv. Model. Earth Syst.*, **11**, 3691–3727, https://doi.org/10.1029/2019MS001829.

Hirschi, J. J., P. D. Killworth, and J. R. Blundell, 2007: Subannual, seasonal, and interannual variability of the North Atlantic meridional overturning circulation. *J. Phys. Oceanogr.*, **37**, 1246–1265, https://doi.org/10.1175/JPO3049.1.

Hochet, A., T. Huck, O. Arzel, F. Sévellec, A. Colin de Verdière, M. Mazloff, and B. Cornuelle, 2020: Direct temporal cascade of temperature variance in eddy-permitting simulations of multidecadal variability. *J. Climate*, **33**, 9409–9425, https://doi.org/10.1175/JCLI-D-19-0921.1.

Hourdin, F., and Coauthors, 2013: Impact of the LMDZ atmospheric grid configuration on the climate and sensitivity of the IPSL-CM5A coupled model. *Climate Dyn.*, **40**, 2167–2192, https://doi.org/10.1007/s00382-012-1411-3.

Itó, K., 1944: Stochastic integral. *Proc. Imperial Academy*, **20**, 519–524, https://doi.org/10.3792/pia/1195572786.

Jamet, Q., W. Dewar, N. Wienders, and B. Deremble, 2019: Spatiotemporal patterns of chaos in the Atlantic overturning circulation. *Geophys. Res. Lett.*, **46**, 7509–7517, https://doi.org/10.1029/2019GL082552.

Jones, D. C., G. Forget, B. Sinha, S. A. Josey, E. J. Boland, A. J. Meijers, and E. Shuckburgh, 2018: Local and remote influences on the heat content of the Labrador Sea: An adjoint sensitivity study. *J. Geophys. Res. Oceans*, **123**, 2646–2667, https://doi.org/10.1002/2018JC013774.

Kanzow, T., and Coauthors, 2010: Seasonal variability of the Atlantic meridional overturning circulation at 26.5°N. *J. Climate*, **23**, 5678–5698, https://doi.org/10.1175/2010JCLI3389.1.

Kirtman, B. P., and Coauthors, 2012: Impact of ocean model resolution on CCSM climate simulations. *Climate Dyn.*, **39**, 1303–1328, https://doi.org/10.1007/s00382-012-1500-3.

Kostov, Y., H. L. Johnson, and D. P. Marshall, 2019: AMOC sensitivity to surface buoyancy fluxes: The role of air-sea feedback mechanisms. *Climate Dyn.*, **53**, 4521–4537, https://doi.org/10.1007/s00382-019-04802-4.

Kravtsov, S., 2020: Dynamics and predictability of hemispheric-scale multidecadal climate variability in an observationally constrained mechanistic model. *J. Climate*, **33**, 4599–4620, https://doi.org/10.1175/JCLI-D-19-0778.1.

LaCasce, J., and J. Pedlosky, 2004: The instability of Rossby basin modes and the oceanic eddy field. *J. Phys. Oceanogr.*, **34**, 2027–2041, https://doi.org/10.1175/1520-0485(2004)034<2027:TIORBM>2.0.CO;2.
Large, W. G., and S. G. Yeager, 2004: Diurnal to decadal global forcing for ocean and sea-ice models: The data sets and flux climatologies. *NCAR Tech. Note NCAR/TN-460+STR*, 105 pp., https://doi.org/10.5065/D6KK98Q6.

Latif, M., and N. S. Keenlyside, 2011: A perspective on decadal climate variability and predictability. *Deep-Sea Res. II*, **58**, 1880–1894, https://doi.org/10.1016/j.dsr2.2010.10.066.

Leroux, S., T. Penduff, L. Bessières, J.-M. Molines, J.-M. Brankart, G. Sérazin, B. Barnier, and L. Terray, 2018: Intrinsic and atmospherically forced variability of the AMOC: Insights from a large-ensemble ocean hindcast. *J. Climate*, **31**, 1183–1203, https://doi.org/10.1175/JCLI-D-17-0168.1.

Lindner, B., 2009: A brief introduction to some simple stochastic processes. *Stochastic Methods in Neuroscience*, C. Laing, and G. J. Lord, Eds., Oxford University Press, 1–28.

Liu, Z., 2012: Dynamics of interdecadal climate variability: A historical perspective. *J. Climate*, **25**, 1963–1995, https://doi.org/10.1175/2011JCLI3980.1.

Lohmann, G., and J. Schneider, 1999: Dynamics and predictability of Stommel's box model. A phase-space perspective with implications for decadal climate variability. *Tellus*, **51A**, 326–336, https://doi.org/10.3402/tellusa.v51i2.12314.

Loose, N., P. Heimbach, H. Pillar, and K. Nisancioglu, 2020: Quantifying dynamical proxy potential through shared adjustment physics in the North Atlantic. *J. Geophys. Res. Oceans*, **125**, e2020JC016112, https://doi.org/10.1029/2020JC016112.

Lorenz, E. N., 1969: The predictability of a flow which possesses many scales of motion. *Tellus*, **21**, 289–307, https://doi.org/10.3402/tellusa.v21i3.10086.

Lozier, S., and Coauthors, 2017: Overturning in the Subpolar North Atlantic Program: A new international ocean observing system. *Bull. Amer. Meteor. Soc.*, **98**, 737–752, https://doi.org/10.1175/BAMS-D-16-0057.1.

Madec, G., and Coauthors, 2012: NEMO ocean engine: Version 3.4. Note du Pôle de Modélisation, Institut Pierre-Simon Laplace (IPSL), 356 pp.

Meehl, G. A., and Coauthors, 2009: Decadal prediction: Can it be skillful? *Bull. Amer. Meteor. Soc.*, **90**, 1467–1486, https://doi.org/10.1175/2009BAMS2778.1.

Monahan, A. H., J. Alexander, and A. J. Weaver, 2008: Stochastic models of the meridional overturning circulation: Time scales and patterns of variability. *Philos. Trans. Roy. Soc.*, **A366**, 2525–2542, https://doi.org/10.1098/rsta.2008.0045.

Msadek, R., K. Dixon, T. Delworth, and W. Hurlin, 2010: Assessing the predictability of the Atlantic meridional overturning circulation and associated fingerprints. *Geophys. Res. Lett.*, **37**, L19608, https://doi.org/10.1029/2010GL044517.

Penland, C., 2003: A stochastic approach to nonlinear dynamics. *Bull. Amer. Meteor. Soc.*, **84**, ES43–ES52, https://doi.org/10.1175/BAMS-84-7-Penland.

Persechino, A., J. Mignot, D. Swingedouw, S. Labetoulle, and E. Guilyardi, 2013: Decadal predictability of the Atlantic meridional overturning circulation and climate in the IPSL-CM5A-LR model. *Climate Dyn.*, **40**, 2359–2380, https://doi.org/10.1007/s00382-012-1466-1.

Pillar, H. R., P. Heimbach, H. L. Johnson, and D. P. Marshall, 2016: Dynamical attribution of recent variability in Atlantic overturning. *J. Climate*, **29**, 3339–3352, https://doi.org/10.1175/JCLI-D-15-0727.1.

Polo, I., J. Robson, R. Sutton, and M. A. Balmaseda, 2014: The importance of wind and buoyancy forcing for the boundary density variations and the geostrophic component of the

AMOC at 26°N. *J. Phys. Oceanogr.*, **44**, 2387–2408, https://doi.org/10.1175/JPO-D-13-0264.1.

Sardeshmukh, P. D., and P. Sura, 2009: Reconciling non-Gaussian climate statistics with linear dynamics. *J. Climate*, **22**, 1193–1207, https://doi.org/10.1175/2008JCLI2358.1.

Sein, D. V., and Coauthors, 2018: The relative influence of atmospheric and oceanic model resolution on the circulation of the North Atlantic Ocean in a coupled climate model. *J. Adv. Model. Earth Syst.*, **10**, 2026–2041, https://doi.org/10.1029/2018MS001327.

Sérazin, G., and Coauthors, 2017: A global probabilistic study of the ocean heat content low-frequency variability: Atmospheric forcing versus oceanic chaos. *Geophys. Res. Lett.*, **44**, 5580–5589, https://doi.org/10.1002/2017GL073026.

Sévellec, F., and A. V. Fedorov, 2013: The leading, interdecadal eigenmode of the Atlantic meridional overturning circulation in a realistic ocean model. *J. Climate*, **26**, 2160–2183, https://doi.org/10.1175/JCLI-D-11-00023.1.

—, and —, 2015: Optimal excitation of AMOC decadal variability: Links to the subpolar ocean. *Prog. Oceanogr.*, **132**, 287–304, https://doi.org/10.1016/j-pocean.2014.02.006.

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

—, and B. Sinha, 2018: Predictability of decadal Atlantic meridional overturning circulation variations. *Oxford Research Encyclopedia of Climate Science*, Oxford University Press, https://doi.org/10.1093/acrefore/9780190228620.013.81.

—, M. Ben Jelloul, and T. Huck, 2007: Optimal surface salinity perturbations influencing the thermohaline circulation. *J. Phys. Oceanogr.*, **37**, 2789–2808, https://doi.org/10.1175/2007JPO3680.1.

—, T. Huck, M. Ben Jelloul, and J. Vialard, 2009: Nonnormal multidecadal response of the thermohaline circulation induced by optimal surface salinity perturbations. *J. Phys. Oceanogr.*, **39**, 852–872, https://doi.org/10.1175/2008JPO3998.1.

—, A. V. Fedorov, and W. Liu, 2017: Arctic sea-ice decline weakens the Atlantic meridional overturning circulation. *Nat. Climate Change*, **7**, 604–610, https://doi.org/10.1038/nclimate3353.

—, H. A. Dijkstra, S. S. Drijfhout, and A. Germe, 2018: Dynamical attribution of oceanic prediction uncertainty in the North Atlantic: Application to the design of optimal monitoring systems. *Climate Dyn.*, **51**, 1517–1535, https://doi.org/10.1007/s00382-017-3969-2.

—, A. C. Naveira Garabato, and T. Huck, 2021: Damping of climate-scale oceanic variability by mesoscale eddy turbulence. *J. Phys. Oceanogr.*, **51**, 491–503, https://doi.org/10.1175/JPO-D-20-0141.1.

Small, R. J., and Coauthors, 2014: A new synoptic scale resolving global climate simulation using the Community Earth System Model. *J. Adv. Model. Earth Syst.*, **6**, 1065–1094, https://doi.org/10.1002/2014MS000363.

Smith, T., and P. Heimbach, 2019: Atmospheric origins of variability in the South Atlantic meridional overturning circulation. *J. Climate*, **32**, 1483–1500, https://doi.org/10.1175/JCLI-D-18-0311.1.

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

—, S. A. Müller, and F. Sévellec, 2020: Tracking water masses using passive-tracer transport in NEMO v3.4 with NEMOTAM: Application to North Atlantic Deep Water and North Atlantic
Subtropical Mode Water. *Geosci. Model Dev.*, **13**, 2031–2050, https://doi.org/10.5194/gmd-13-2031-2020.

Sura, P., and P. D. Sardeshmukh, 2008: A global view of non-Gaussian SST variability. *J. Phys. Oceanogr.*, **38**, 639–647, https://doi.org/10.1175/2007JPO3761.1.

Tziperman, E., and P. J. Ioannou, 2002: Transient growth and optimal excitation of thermohaline variability. *J. Phys. Oceanogr.*, **32**, 3427–3435, https://doi.org/10.1175/1520-0485(2002)032<3427: TGAOEO>2.0.CO;2.

Uhlenbeck, G. E., and L. S. Ornstein, 1930: On the theory of the Brownian motion. *Phys. Rev.*, **36**, 823–841, https://doi.org/10.1103/PhysRev.36.823.

Vidard, A., P. A. Bouttier, and F. Vigilant, 2015: NEMOTAM: Tangent and adjoint models for the ocean modelling platform NEMO. *Geosci. Model Dev.*, **8**, 1245–1257, https://doi.org/10.5194/gmd-8-1245-2015.

Zanna, L., and E. Tziperman, 2008: Optimal surface excitation of the thermohaline circulation. *J. Phys. Oceanogr.*, **38**, 1820–1830, https://doi.org/10.1175/2008JPO3752.1.

Zhang, W., B. Kirtman, L. Siqueira, A. Clement, and J. Xia, 2021: Understanding the signal-to-noise paradox in decadal climate predictability from CMIP5 and an eddying global coupled model. *Climate Dyn.*, **56**, 2895–2913, https://doi.org/10.1007/s00382-020-05621-8.