Dynamical attribution of North Atlantic interdecadal predictability to 6 oceanic and atmospheric turbulence under 7 realistic and optimal stochastic forcing

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 timescales. In a global ocean model, we present i.) an optimisation framework to compute the most efficient noise patterns to generate uncertainty and ii.) a uniquely inexpensive, dynamical method for attributing sources of ocean uncertainty to internal (mesoscale eddy turbulence) and external (atmospheric) 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. We demonstrate that optimal noise patterns target features of the underlying circulation such as the North Atlantic Current and deep water formation regions. We then show that noise forcing in the actual climate system stimulates these patterns with various degrees of efficiency, ultimately leading to the growth of error. We reaffirm the established notion that higher frequency variations 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 which contribute the most to ocean error, accounting for up to 60% after 60 years of growth in the case of 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. 16 17 18 19 20 21 22 23 24

oceanic prediction error.

49
As the slow component of the climate system, the ocean is key to predicting variations on 50 timescales of seasons or longer. However, the ocean is now known to exhibit substantial variability which has solution (for initial condition zero) 138 u(t 0 , t 1 ) = where u is the surface temperature, t 0 and t 1 are the initial and final time, σ 2 is the variance of (2) is thus an Itô integral (Itô 1944). It may be noted that the response is 143 an Ornstein-Uhlenbeck process (Uhlenbeck and Ornstein 1930), such that variability generation follows the autocovariance function: 145 Cov(u(t 0 , t 1 ), u(t 0 , t 2 )) = σ 2 2λ e −λ|t 2 −t 1 | − e −λ(t 2 +t 1 −2t 0 ) .
This autocovariance function is weakly stationary in the limit t 0 → −∞ and so corresponds via where ω is the time frequency and S the PSD.

149
In this simple framework, the ocean therefore low-pass filters spectrally constant (white noise) 150 surface heat fluxes, producing a frequency spectrum which is constant (i.e., white noise) in the 151 limit of low frequency (ω λ) and follows an inverse square law (i.e., red noise) in the limit of 152 high frequency (ω λ). The transition frequency is determined by the ocean adjustment timescale 153 (i.e., λ). We will return to these classical results concerning Ornstein-Uhlenbeck processes in 154 Section 2c.

155
Although a useful first-order representation of the evolution of unpredictable surface temperature Noting that d |ū = N (|ū , t) dt and neglecting higher order terms leads to the complementary 175 equation of (5). In this context, A(t) is the Jacobian of the nonlinear system with respect to the 176 ocean state: The (zero initial condition) solution to (5) is given by where Ψ(t 1 , t 0 ) is the propagator matrix [the scalar Ψ(t 1 , t 0 ) = e −λ(t 1 −t 0 ) in the univariate case of Beginning from the last formula, we can diagnose the covariance between any two scalar-valued metrics of the ocean state which are linear. These metrics can be defined by the co-vectors |F 1,2 183 where the scalar product F 1,2 |u = u|F 1,2 are the Euclidean inner products of the co-vectors 184 and the ocean state vector anomaly. We have where s represents time. A multi-dimensional generalisation of Itô's isometry may be applied to 186 this expression (e.g., Section 3.6 of Duan and Wang 2014). In particular, the Itô integral terms may 187 be written as non-anticipatory (left) Riemann sums such that the right hand side of (9) becomes where i, j, k are discrete increment indices, and K is the total number of discrete increments.

190
Applying a transpose and Fubini's theorem: We note that, ∀ i = j, the increments of the Wiener processes do not overlap and so are independent 192 by definition, reducing the expression to a single sum in which the central outer product corresponds to a diagonal matrix, as the vectors are elementwise 194 independent. As Wiener increments are normally distributed as W (t k+1 − t k ) ∼ N (0, t k+1 − t k ), 195 in their infinitesimal limit the equation becomes 196 Cov( F 1 |u(t 0 , t 1 , F 2 |u(t 0 , t 1 ) = t 1 t 0 F 1 |Ψ(t 1 , t)ΣΨ † (t, t 1 )|F 2 dt.
Note that our solution generalizes the result heuristically derived by Sévellec et al. (2018). Similarly 197 to their approach, we remark that while it is standard to diagnose the variance evolution of a metric 198 by propagating many realisations of (8) as an ensemble and considering its spread, (13)    To determine the optimal Σ, we apply two constraints to the optimal variance source: its global 223 average has fixed amplitude, and any two points which are not independent have a correlation of ±1.

224
The former implies that the stochastic process has finite power (corresponding to band-limited white 225 noise), while the latter assumes that if two points covary, they must do so completely constructively 226 (as would be optimal). We begin by considering the general case, where the stochastic process is As outlined above, we partition the stochastic process into N regions such that points within the 233 regions are perfectly covarying, but are independent of points in other regions. Equivalently, we region), and define a binary projection B i ∈ R n×m i such that Following (8), the evolution of the state vector in response to stimulation in the i th region is The implication is that in the region, a single Wiener process is "shaped" 241 by a pattern of local amplitudes |L i .

242
In order to determine the optimal shape of this pattern, we utilise the method of Lagrange where S i ∈ R m i ×m i is a (diagonal) volumetric weighting matrix. The corresponding Lagrange 247 function can be expressed as where the scalar γ i is the Lagrange multiplier. Maximizing the Lagrangian leads to effect. We note that left multiplication of (18) by L * i | S i results in 2) L

269
The above derivation applies to the case of N perfectly correlated independent regions, but we 270 may consider two specific cases of this in order to imitate conditions similar to the atmospherically 271 forced and eddy-driven variability felt by the ocean. In particular, we consider the two limiting is uncorrelated between all variables and locations, is taken as an idealized representation of small-276 scale noise in the ocean (i.e., noise induced by subgrid processes). These cases correspond to 277 solving a single eigenvalue problem vs. solving dim(|u ) (trivially scalar) eigenvalue problems.

278
In particular, for N = 1, the sole projection matrix is the identity matrix B 1 = I, while for 279 N = n = dim(|u ), the projection matrices become the standard basis vectors B i = |e i (i.e., e i 280 projects a scalar to the i th location of the full state vector).

281
In the former (everywhere perfectly covarying) case, (22) becomes where |L opt is the leading eigenvector of The latter (everywhere uncorrelated) case corresponds to the condition that every point is inde-284 pendent, and Σ opt is diagonal. The associated eigen"vector" problems are scalar, such that the 285 eigenspace is infinite. All terms in (22) are now scalars such that |L opt i can be seen to cancel, 286 while the matrices S i may be written as S i . Ultimately, where, solving (18) with B i = |e i , the eigenvalues γ opt i are trivially the diagonal elements of The sum of the eigenvalues is also the trace of this (scaled outer product) matrix, and is thus given 289 by the corresponding inner product. Therefore, from (24) the optimal stochastic covariance matrix 290 in the completely uncorrelated case is (where the diag[·] operator corresponds to the diagonal matrix with the same diagonal). We respec-292 tively use these two limiting cases to explore theoretical variance linked to idealized atmospheric 293 forcing (assuming perfect correlation everywhere over the surface and zero noise in the interior) and 294 ocean internal subgrid fluxes (assuming noise everywhere, with zero correlation between locations 295 and variables).

296
A useful metric of the OSP is the ratio of the output variance to the input variance A * = 297 Var( F |u(t 0 , t 1 ) / 2 , which we term the amplification factor. Notably, for the globally perfect 298 covariance case, this is simply the associated eigenvalue For the globally decorrelated case, is the sum of the eigenvalues. Our considerations so far have involved stochastic forcing with varying levels of spatial coherence, 303 but which is serially decorrelated (therefore band-limited white noise). While this allows an 304 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 306 the diagnosed case. Indeed, while white noise is typically considered an acceptable representation 307 of atmospheric variability (which decays on timescales much shorter than those of the oceanic 308 large scale; Hasselmann 1976), the ocean mesoscale eddy field evolves much more slowly (e.g.,

309
Chelton et al. 2007). To realistically represent this using diagnosed fluxes, we therefore extend our 310 framework to include temporally correlated stochastic forcing. We consider again the Ornstein-

311
Uhlenbeck case, which is a simple example of a temporally correlated stochastic process.

312
We begin by modifying (5) where, as before, |u defines the state vector anomaly, A defines the system's linear dynamics

317
[for instance via the Jacobian of a corresponding nonlinear system, as in (7)], and, in contrast 318 to the previous cases, |X is the forcing from continuous, spatiotemporally correlated stochastic 319 processes. The zero-initial-condition solution is given by where the complementary equation and therefore the propagator matrix, Ψ(t 1 , t 0 ), are notably 321 identical to (8). As in (9), we seek the covariance between two metrics of the state vector, given by 322 Cov( F 1 |u(t 0 , t 1 , F 2 |u(t 0 , t 1 ) = Ornstein-Uhlenbeck processes [such as (2) with t 0 → −∞], a multivariate generalization of (3) 326 gives where λ is a diagonal matrix of reciprocal e-folding times of the anomalous fluxes at each location, 328 and LL † = Σ is their spatial covariance matrix. As these quantities can be diagnosed from an 329 appropriate dataset, we can use this formulation to diagnose the variance growth.

330
In the proceeding section we diagnose (from realistic models) λ and Σ for the cases of external 331 (atmospheric; λ E , Σ E ) and internal (oceanic turbulent mesoscale eddy driven; λ I , Σ I ) turbulent 332 fluxes, assessing the appropriateness of the Ornstein-Uhlenbeck representation. We then proceed 333 to attribute the variance of different metrics in response to these sources, which, following (29) 334 and assuming independence between the internal and external components is given by The variance may be broken down further still, by writing the covariance matrices as the sum of their the λ E , Σ E terms), and, in the latter case, the separate contributions of the covarying zonal and 339 meridional momentum fluxes. The final term of 31 can accordingly be split into: ) are the external noise properties for the buoyancy, and zonal and merid-341 ional momentum fluxes, respectively, Σ u,v E is for the zonal and meridional covariance term.

342
Finally, in addition to separating the variance into contributions from different variables, we note 343 that we can also isolate contributions from different regions of space. The inner products of (31) 344 represent spatial integrals of local contributions to the total variance (integrated over volume in the 345 internal case and over area in the external case). An alternative formulation of (31) is therefore where V I and V E are continuous functions representing the respective internal variance contribution

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

401
To evaluate the goodness of fit of the Ornstein-Uhlenbeck process representation, we compare advection equation (at high Péclet number, such that diffusive processes can be neglected) reads where T and V are the scalar and tridirectional vector fields of temperature and of velocity, 418 respectively, ∇· is the tridirectional gradient operator, ·|· is the inner product, · is a tridirectional 419 spatial averaging operator, and · is its associated spatial fluctuation. This separation is such that 420 the lower-resolution model (LRM) is able to resolve temperatures at the scale of the spatial average 421 (e.g., T ), while the higher-resolution model (HRM) resolves the sum of the spatial average and its 422 fluctuation (e.g., T = T + T ). We are interested in the mean effect of the small scale on the large 423 scale following application of the spatial averaging operator. Applying this operator, the equation As before, we consider the large-scale climatological cycle to be common to the HRM and the 426 LRM, so we separate (35) into a trajectory (·) and a temporal fluctuation (· ): where ∂ t T + V |∇ T = − V |∇ T is the trajectory component common to both LRM and HRM.

428
The unresolved component in the LRM is therefore where smaller, second order time-fluctuating terms ( V |∇ T ) have been neglected. The flux 430 terms on the left-hand side represent large-scale interactions between the temporal mean and its 431 fluctuations, while the right hand side describes the temporal fluctuation of small-scale fluxes.

432
We note that the latter term also contains interactions between the climatology and fluctuations, 433 as can be seen by separating its interior velocity and temperature components into their own where similar considerations may be made for the internal eddy-driven salt flux (F I S ). 439 We apply this approach to determine the turbulent eddy heat and salt fluxes unresolved in our Atlantic. In all cases, monthly, annually, and decadally averaged quantities are considered. We now consider (using the limiting cases of Section 2b) the spatially correlated external and 487 decorrelated internal OSPs of the metrics of Section 3d in our linearized ocean model (Section 3a).

488
The sensitivity of the metric to different potential sources of variability is indicated by the amplifi-489 cation factor (Table 1)

515
The optimal OHC perturbation in the uncorrelated case ( Fig. 5a and b) shows sensitivity to in the real climate system, as derived in Section 2c. Following (31) and (32), application of 535 each component of the stochastic forcing separately allows the resultant variance evolution to be 536 partitioned accordingly (Fig. 6). There is a substantial difference between the nature of month-and 537 decade-averaged transport metrics, both in the variance amplitude and in the impacts of different 538 sources, as in the OSP case (shown by the amplification factors of and MHT, where it is the external buoyancy contribution which is the slowest to develop.

559
Following (33), we consider the spatial distributions of these contributions to the 60-yr variance 560 for the annually averaged case, within the transition between the two discussed (month-and decade-561 average) cases ( Fig. 7; where the zonal and meridional momentum flux covariance contribution 562 is not shown). There is generally a high level of agreement between the patterns shown in the 563 optimal case (i.e., what the ocean "wants"; Figs 3, 4, and 5) and the realistic case (Fig. 7). This 564 is linked to the overall relative constant shape of the realistic forcing (i.e., what the ocean "gets"; 565 Fig. 1). Although we remind of the contrast between the two frameworks (i.e., white vs. temporally 566 correlated noise) when making any such comparisons.

567
In particular, volume and heat transport variability are primarily driven by ocean internal buoy- Applying the considerations of Section 4a1 to the subpolar region, differences emerge in the 589 amplitude of the response to the optimal stochastic forcing ( Table 2). For subpolar MVT, the 590 correlated surface OSP is much more effective at generating variability than in the subtropics, Under prescribed, realistic sources of variability, the subpolar region is dominated by external 623 forcing (Fig. 11), which accounts for up to 94% of the total variance after 60 years in the case 624 of month-averaged heat transport. As in the subtropics, the meridional transport metrics exhibit 625 a regime shift when moving from month-averaged quantities (up to 86% momentum-driven) to This has allowed us to compare the commonalities between the optimal and actual cases (albeit 671 in a limited way, given their differences in spatiotemporal correlation). We have further been 672 able, in the diagnosed case, to dynamically attribute variability to its origins. The latter ability 673 notably forgoes the more typical ensemble attribution approach, which generally necessitates many  avg. T 2. As in Table 1, but for subpolar OSPs, (whose spatial distributions are shown in Figures 8,9, and 10).   Fig. 7, but for subpolar ocean metrics. We note again that the differing contributions 1011 (as shown in Figure 11) lead to large differences in the color scales between panels. . . . . 63

LIST OF FIGURES
Var (      F . 12. As in Fig. 7, but for subpolar ocean metrics. We note again that the differing contributions (as shown in Figure 11) lead to large differences in the color scales between panels.