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

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


of 11
Recent studies have argued that the predominant patterns of Atlantic Multidecadal Variability (AMV) in climate simulations featuring realistic ocean general circulation models (OGCMs) can be recreated by coupling a realistic atmosphere to a time-invariant "slab" ocean (Cane et al., 2017;Clement et al., 2015Clement et al., , 2016, suggesting that these patterns are purely passive. In this slab ocean case, commonalities with fully active ocean simulations can only be established statistically. On the other hand, the aforementioned passive tracer approach (by propagating a passive "temperature" tracer initially coincident with the active temperature field in a single simulation and considering their divergence) provides a more clear separation of the response into its passive and dynamical components. Nevertheless, statistical slab-OGCM comparisons remain the de facto standard for determining the role of the ocean in near-term regional low-frequency variability (Dommenget, 2010;Dommenget & Latif, 2002;Delworth et al., 2017;Wang & Dommenget, 2016;Zhang, 2017).
In this study, we present an alternative approach to the question of regional heat content variability, using an adjoint model. Unlike a conventional model, which integrates anomalies forward in time, an adjoint model describes the sensitivity of a quantity of interest (here heat content) to past changes (here stochastic atmospheric forcing), establishing causes, rather than effects (Errico, 1997). As such, the method describes all possible perturbations leading to a particular ocean response, otherwise requiring a theoretically "infinite" ensemble of perturbed simulations. This has been leveraged to attribute the sources of temporal ocean variability in response to historical atmospheric forcing (Kostov et al., 2021;Pillar et al., 2016;Smith & Heimbach, 2019) and to establish the growth of oceanic variance in response to a representation of stochastic atmospheric forcing (Sévellec et al., 2018).
We use this approach to isolate the role of the ocean in modeled heat content variability, by projecting a diagnosed stochastic representation of atmospheric buoyancy and momentum fluxes onto passive and active surface adjoint sensitivity fields. In the passive case, buoyancy anomalies cannot change the circulation.

Method and Diagnostics
To characterize low-frequency ocean variability, Hasselmann (1976) and Frankignoul and Hasselmann (1977) developed an idealized, single-variable stochastic model of ocean surface temperature in response to random heat fluxes. These atmospheric fluxes can be seen as a continuous stream of small disturbances to ocean surface temperature, which accumulate and are slowly "forgotten". This can be represented as A large body of work has generalized Equation 1 to mutivariate systems featuring multiple interacting variables and locations, ranging from box models (Cessi, 1994;Chang et al., 2004;Griffies & Tziperman, 1995;Grötzner et al., 1999; to GCMs (Hawkins & Sutton, 2009;Sévellec et al., 2018;Tziperman et al., 2008). In the higher order case, different variables and locations may all be represented by a single anomaly state vector: |u . Here we utilize the bra-ket notation of Dirac (1939), whereby vectors are written as "kets," for example, |b , and co-vectors as "bras," for example, a| , such that bra-ket pairs become either scalar through the Euclidean inner product, for example    a b | c or matrices through the outer product, for example | | b a   C . We may also represent more involved linear processes than simple exponential decay using (e.g., Farrell & Ioannou, 1996a, 1996b) the propagator,  , of the ocean model which linearly maps anomalies | ( ) u 0  to their later states | ( ) t e in the univariate case of Equation 1. This reads: where |W  is a vector of independent standard-normal Wiener processes and  † E Σ LL is a covariance matrix (describing the stochastic flux intensity similarly to 2 E L in the univariate case of Equation 1 but also accounting for coherence between locations and variables). As before, | ( ) u 0  is assumed zero-valued.
where  AP E is the covariance between fully active and purely passive version of the model denoted by  E A and  E P , respectively, and  is the expectation of a stochastic Itô process. This leads to expressions for the variance of the fully active and purely passive component at time t : These equations describe the level of variance of the ocean heat content obtained after a time t in response to stochastic forcing starting from rest in the fully active (  2 A E ) and purely passive (  2 P E ) cases. Of note is that the disappearance of the stochastic terms d|W  means that it is not necessary to propagate any particular realization of random noise to obtain these properties-the variance of an "infinite" ensemble of such realizations can be determined using just the properties of the noise and a single run of the linear model. These will asymptotically tend towards their associated climatological heat content variance. The covariance describes how much of this variance is common to both, and can be normalised to give a Level of Agreement (LoA) between the purely passive and fully active cases, which we define as LoA is unity at a given time, it is taken that anomalous heat content variation in the fully active ocean has been entirely controlled by purely passive processes.
These diagnostics have three requirements. Firstly, a linearized OGCM is needed to provide the propagator  E A and its adjoint  † E A . Secondly, this propagator requires an isolated purely passive component  E P and its adjoint  † E P . The model, its adjoint, and the purely passive configuration are described in Section 3.1. Lastly, we require a stochastic representation  E of surface fluxes. We diagnose this from a coupled climate model (also described in Section 3.1). In particular, we assume that buoyancy and momentum flux anomalies from the coupled simulation climatology follow a band-limited (therefore finite power), spatially covarying Gaussian white noise. At each location, the power spectral density (PSD) of the flux anomalies is therefore assumed constant up to a few days, and zero at higher frequency. The cutoff is determined by the e-folding decorrelation timescales of the fluxes (Figure 1, contours). We also have an implicit low-frequency limit imposed by the 20-year length of the coupled simulation. The elements of  E are then given by the (effectively constant) PSD averaged over this band.
It is important to remark on linearity and independence, which allow for further decomposition of the above diagnostics. As the model propagators are linear, we can consider the fully active model  E A to be the sum of the purely passive model  E P and a dynamical-only component  E D which encompasses just the feedback terms. Furthermore, the propagation of multiple quantities is equal to the propagation of their sum by linearity:

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

Model Description
Our stochastic representation  E in Equation 3 is constructed from thermal, haline, and zonal and meridional momentum fluxes diagnosed from a coupled climate model ( Figure 1). Specifically, a 20 year simulation using the IPSL-CM5A-LR coupled model was run in its CMIP5 pre-industrial control configuration (Dufresne et al., 2013) with daily average output. Atmospheric variations at higher or lower frequencies than those captured by this simulation were taken not to influence heat content variance on the spatiotemporal scales considered here. The model was chosen as its ocean component is NEMO (v3.2) with its ORCA2 global configuration (  2 E nominal resolution with 31 vertical levels), similarly to our linearized ocean model (described below). The atmospheric component is the LMDZ5a model, with 3.    75 1.9 E horizontal resolution and 39 vertical levels (Hourdin et al., 2013). The linear ocean model which provides the propagator  E A of Equation 3 is NEMOTAM (Vidard et al., 2015), which is derived from NEMO v3.4 (Madec, 2012) and is used in its ORCA2-LIM configuration. The model configuration is similar to that detailed in Stephenson et al. (2020), which also discusses the implementation of the purely passive configuration in detail. The nonlinear model, which provides the simulation about which NEMOTAM is linearized, is forced by a single representative year (CORE normal year forcing, phase 1, v2.0; Large & Yeager, 2004). This allows the ocean response described by  E in Equation 2 to be that of a "typical" year, cleanly separated from interannual surface variations, which were applied separately through the matrix  E in Equation 3.

Results
We now apply the derivations of Section 2 to attribute the generation of variance in heat content (  | E F in Equations 3 and 4) in the fully active simulation to its different sources. We evaluate heat content over three depth ranges (10 m, 1,500 m, and full-depth), which effectively correspond to sea surface temperature, heat content in the upper ocean, and the total heat content, respectively. We also consider both the global ocean and a seven-region partition of it (Figure 2,  Our analysis reveals that the LoA between purely passive and fully active heat content variance after 60 years varies significantly depending on the depth extent and geographical region (Figure 2, bars). The LoA is extremely high for sea surface temperature variance in all regions. This ranges from 92% in the intertropical Pacific to over 99% in the Southern Ocean, with a majority stimulated by buoyancy forcing. This implies that the purely passive uptake of heat controls temperature variability at the surface. There is a dramatic reduction in agreement when heat content is computed over a thicker layer. For the upper-1,500 m heat content, variance common to both the purely passive and fully active simulations accounts for as little as 31% in the case of the Indian Ocean, and just over half (52%) globally. The nature of stimulation of the purely passive component also changes over this depth range, shifting to a primarily wind-driven regime for all regions except the Arctic Ocean.
When heat content is defined over the full depth, it generally follows similar patterns to upper-1,500 m heat content, with notable exceptions in the North Atlantic and Arctic oceans. For those basins, another dramatic reduction in correspondence between the purely passive and fully active simulations occurs, with the LoA reducing to 27% and 25%, respectively. More subtle reductions can be seen elsewhere, and only in the North Pacific and Southern Ocean does the purely passive component still dominate the fully active simulation at full depth. It is worth noting the substantial impact (>50%) of purely passive wind effects in these regions. in response to each. LoA is shown for the three cases (surface layer-corresponding to SST, upper 1,500 m, and full-depth heat content). Largest bar plot shows the case for the total global ocean heat content variance, smaller inner plots show regional values. Thinner dashed black lines signify a LoA of 50%. Black solid lines on the map mark the boundaries of the regions in our definitions.

of 11
While the LoA provides a useful quantification of the ultimate role of the purely passive component of the ocean, it does not describe in detail the differences between the purely passive and fully active simulations (e.g., the timing of the variance growth or its source location). To tackle this question, we consider the time-evolving variance growth for each, along with its components (Figure 3). We focus on the full-depth case, where differences between these components are greatest. Similar decompositions have been considered for surface ( Figure S1) and upper ocean ( Figure S2) cases, and exhibit similar (but less significant) behavior.
The temporal evolution of the variability reveals that the purely passive and fully active simulations differ in both magnitude and timing. As discussed in Section 2, linearity permits the decomposition of the fully active model into the sum of the purely passive component and a remaining dynamical-only component. The difference between evolving variance in the fully active model and in the purely passive model (Figure 3, solid and dotted lines, respectively) can thus be attributed to internal ocean feedbacks belonging to this dynamical-only component, which are not always constructive. Indeed, the variance in the fully active simulation is often weaker than that of its purely passive counterpart. This suggests that certain behavior is possible only in the purely passive case, and is canceled out by the dynamical-only term in the fully active simulation. This is particularly visible for heat content variance in the Indian and Southern Oceans (dominated by wind stress). There, after two decades, most of the variance growth of the purely passive Figure 3. Evolution of full-depth heat content variance in response to stochastic surface forcing in the purely passive (dotted lines) and fully active (solid lines) simulations. The difference between these lines is linked to the dynamicalonly component, which may act destructively (passive > active) or constructively (active > passive). Thinner lines show separately the buoyancy-forced (red) and wind-driven (green) components. Stars mark the point at which 50% of the final (60-year) variance is reached. 7 of 11 component stimulated by wind stress is canceled by the dynamical-only component. An example of similar behavior is provided by Cronin and Tozuka (2016), who demonstrate that Ekman transport is determined not only by wind stress and latitude (as in the classical analysis of Ekman, 1905), but also local geostrophic shear. In this perspective, Ekman transport has both a purely passive and dynamical-only component, which can act against each other.
As a measure of the rate at which the climatological variance is approached, we consider the time taken for the full-depth variance in each simulation to reach half of its final (60 year) value (Figure 3, stars). Following on from the previous discussion, the dynamical-only momentum component in the Indian Ocean acts to reduce the final variance achieved in the active case, thereby reaching half of this value in only 7 years, as opposed to 19 years in the purely passive simulation. At the opposite extreme, for the Arctic and North Atlantic, the dynamical-only contribution slows the variance evolution substantially. This effect is ultimately dominated by the buoyancy component, and is made up of similar contributions from temperature and salinity (not shown). Resultantly, in the North Atlantic, half of the final value is reached in only 3 years in the purely passive simulation, compared with 21 years in the fully active case.
The source of this continued growth in the active North Atlantic, even after the purely passive component appears to have saturated, corresponds to a change in the nature of the ocean response to buoyancy stimulation in the fully active simulation after 10 years. To determine the origin of this, we consider separately the surface distribution of the variance accumulated during the first 10 years (Figures 4a-4d) and from 10 to 60 years (Figures 4e-4h). This is determined from the elementwise computation of the variance, prior to summation, as outlined in Section 2.
In the first decade, the passive and active simulations maintain a high LoA (above 75%) and their spatial patterns are similar. Focusing on buoyancy forcing, the relatively focused region reflects the model's deep water formation site, as described in the passive tracer study of Stephenson et al. (2020). The difference between the fully active and purely passive distributions is the dynamical-only contribution (Figure 4c, contours). This corresponds to a large-scale dipole. The negative peak of the dipole overlies the positive contribution by the purely passive component, having a slight compensating effect (Figure 4c). On decadal timescales, positive contributions to variance growth in both the purely passive and dynamical-only components coincide in location, and so the two components act constructively (shading and contours in Figure 4g).
The primary difference between stimulation by wind for the fully active and purely passive components in the first decade is the intensity of the induced variance (Figures 4b and 4d). Both components are dominated by Ekman transport across a zonal band defining the region's boundary (  35 E N), with an alternating pattern at the grid scale stimulating a mechanism apparently unique to the passive ocean. The addition of the negative dynamical-only component reduces the intensity of this pattern, however. Also notable in the fully active case is a seemingly persistent (Figures 4b and 4f) stimulation of variance at the subtropical-subpolar inter-gyre interface, as well as stimulation (both positive and negative) in coastal regions of the eastern North Atlantic and Greenland Sea.

Discussion and Conclusions
We have considered the stimulation of variance in ocean heat content by surface atmospheric noise. We evaluated heat content over a range of different regions and depths in a linearized global ocean model, comparing purely passive and fully active realizations of the ocean model. In the purely passive framework, temperature anomalies either arise due to random surface heat fluxes (and can be passively transported by the mean flow), or due to random surface momentum fluxes (which redistribute existing heat). However, these resulting temperature anomalies are unable to modify the ocean circulation.
In contrast to the established techniques of using a passive tracer (e.g., Banks & Gregory, 2006;Garuba & Klinger, 2016Marshall et al., 2015;Xie & Vallis, 2012) or a slab ocean model (e.g., Clement et al., 2015;Dommenget, 2010;Dommenget & Latif, 2002;Wang & Dommenget, 2016) to investigate the role of the ocean, we have utilized a novel adjoint-based approach (Sévellec et al., 2018). The use of an adjoint model has uniquely allowed us to causally attribute heat content variance to different variables, times, and locations at the surface, by projecting onto surface sensitivity fields a representative stochastic representation of atmospheric fluxes diagnosed from a coupled climate model.
Our findings for the surface ocean (i.e., sea surface temperature) are that at least 92% of the variance in the fully active simulation is in agreement with its purely passive component. This is consistent with studies which suggest that oceanic dynamics are not needed to generate surface decadal variability (e.g., Cane et al., 2017;Clement et al., 2015Clement et al., , 2016. While our passive and active variance patterns may in some cases express a high (normalised) LoA, a purely passive model could still misestimate the amplitude of heat content variance, as the purely passive component can be offset or reinforced by the corresponding dynamical-only component in a fully active ocean. In our experiments, this discrepancy occurs over greater depth ranges (Figures 3 and S2), but the amplitudes of surface ocean responses in the passive and active cases are roughly equivalent ( Figure S1), as suggested by Clement et al. (2015) in the North Atlantic, for instance.
The dynamical redistribution of existing heat by currents arising from buoyancy anomalies has been shown in past studies to substantially impact heat uptake (e.g., Banks & Gregory, 2006;Xie & Vallis, 2012), particularly in the North Atlantic. However, we have shown that the passive redistribution of the existing heat reservoir by wind anomalies is often more important in the context of heat content variability. This leads to a driving role for the passive component over several regions and depths, in particular the deep Southern and North Pacific oceans. Nevertheless, the full-depth North Atlantic also stands out here as having an important role for ocean feedbacks, with the dynamical-only component acting to slow the growth of heat content variance. We considered the time taken to reach 50% of the variance at the end of the simulation (as a conservative estimate of the time taken to reach half the asymptotic variance), and found that the fully active model takes 7 times as long (21 years) to reach this point as the purely passive simulation (3 years) in this region. This has potential consequences for climate predictability, as the variance growth can also be seen as the accumulation of error following model initialization (Sévellec et al., 2018). The time taken to reach half of the climatological variability is often taken as a measure of the upper limit of predictability, beyond which noise dominates the predictable signal (e.g., Griffies & Bryan, 1997;Grötzner et al., 1999). The reason for this delay in the fully active North Atlantic is a shift in the response to buoyancy forcing. On sub-decadal timescales, the dynamical-only component slows variance growth, before sustaining it on 9 of 11 timescales greater than 10 years, resulting in an "S"-shaped growth curve. In exploring the spatial distribution of the components of the fully active simulation, we have observed a basin-scale dipole pattern in the North Atlantic. These patterns echo earlier sensitivity studies of the region in predecessors of our model (e.g., Sévellec & Fedorov, 2017). These studies relate North Atlantic heat content sensitivity to an ocean-only mode of variability in which heat content and AMOC anomalies feed back on each other via basin-scale thermal Rossby wave propagation (Sévellec & Fedorov, 2013) consistently with observations of the AMV (Frankcombe et al., 2009).
There are a number of considerations which are not accounted for in our approach. First, our conclusions are likely oversimplified by our use of atmospheric variability sources alone in a linear, laminar model. In a recent ensemble study at eddy-permitting resolution, Sérazin et al. (2017) suggested that a substantial portion of ocean heat content variability is intrinsic, generated by chaotic, nonlinear processes within the ocean. This suggests that we underestimate the role of the dynamical-only component by restricting it to large-scale, laminar feedbacks. This has been addressed in a separate study (Stephenson & Sévellec, 2021). In addition, the role of coupling in the stimulation of interdecadal variability is an entire field of research on its own (e.g., the review of Liu, 2012). Here, our model uses an uncoupled ocean and a stochastic representation of the atmosphere. This limits the conclusions of our work, in particular for sea surface temperature (where the surface boundary conditions have more impact). Furthermore, our stochastic representation is of limited bandwidth, effectively averaging the power spectrum of a two-decade coupled simulation. The result is a stationary (although globally coherent) white noise representation of daily-to-bidecadal atmospheric variability. We emphasize, however, that these simplifications have allowed us to use an adjoint ocean model to causally attribute the surface sources of heat content variability exactly, and with limited computational expense, an approach which offers several unique advantages of its own.

Data Availability Statement
These model outputs can be found at https://doi.org/10.5281/zenodo.4300471 while code used for our adjoint simulations and analysis may be found at https://doi.org/10.5281/zenodo.4721438.