## Abstract

The authors identify an interdecadal oscillatory mode of the North Atlantic atmosphere–ocean system in a general circulation model (GFDL CM2.1) via a linear inverse model (LIM). The oscillation mechanism is mostly embedded in the subpolar gyre: anomalous advection generates density anomalies in the eastern subpolar gyre, which propagate along the mean gyre circulation and reach the subpolar gyre center around 10 years later, when associated anomalous advection of the opposite sign starts the other half cycle. As density anomalies reach the Labrador Sea deep convection region, Atlantic meridional overturning circulation (AMOC) anomalies are also induced. Both the gyre and AMOC anomalies then propagate equatorward slowly, following the advection of density anomalies. The oscillation is further demonstrated to be more likely an ocean-only mode while excited by the atmospheric forcing; in particular, it can be approximated as a linearly driven damped oscillator that is partly excited by the North Atlantic Oscillation (NAO). The slowly evolving interdecadal oscillation significantly improves and prolongs the LIM’s prediction skill of sea surface temperature (SST) evolution.

## 1. Introduction

The Atlantic meridional overturning circulation (AMOC) has attracted enormous attention during recent decades, owing to the AMOC’s potential predictability and significant impacts on global climate. We would like to refer the reader to Buckley and Marshall (2016) for a comprehensive review. In studying the AMOC, previous studies adopted various approaches, including observations (Smeed et al. 2014), analytical models (Griffies and Tziperman 1995), hindcast ocean-only models (Eden and Willebrand 2001; Bailey et al. 2005; Böning et al. 2006), and fully coupled general circulation models (GCMs) (Kwon and Frankignoul 2012, 2014; Tulloch and Marshall 2012; Danabasoglu et al. 2012). Some consensus has been reached in the literature, for example, in terms of the contributions of the North Atlantic Oscillation and subarctic deep-water formation rates to the variability of the AMOC. However, considerable confusion still exists. While many GCMs exhibit interdecadal AMOC variability (Sévellec and Fedorov 2013; Kwon and Frankignoul 2014; Yoshimori et al. 2010), there is yet no unified understanding of the mechanisms behind this variability. Moreover, the slow equatorward propagation of circulation anomalies seen in many models (Eden and Willebrand 2001; Kwon and Frankignoul 2014; Zhang 2010), apparently inconsistent with the traditional Kelvin wave argument, is attributed either to the inappropriate model representation of Kelvin waves, the interior pathway (Zhang 2010), or to the slow advection of density anomalies (Gerdes and Köberle 1995).

Additional complexity exists with regard to the relationship between the AMOC and other processes. For example, hindcast ocean-only simulations reveal a close connection between variations in the subpolar gyre and the AMOC via Labrador Sea convection anomalies: strengthening (weakening) of the subpolar gyre and AMOC are seen as a delayed response to enhanced (suppressed) convection (Eden and Willebrand 2001; Böning et al. 2006). However, Zhang (2008) relates a weakened subpolar gyre to a strengthened AMOC, defined as the maximum meridional streamfunction at 40°N. Meanwhile, Lozier et al. (2010) derive opposing AMOC changes in the regions of the subpolar and subtropical gyres from hydrographic data. This complexity and confusion in the geographic linkages of the AMOC and its fingerprints on other variables is not unexpected considering the large spatial scale of AMOC, but it would be helpful for understanding AMOC mechanisms if we could combine and reconcile evidence from previous studies.

In this study, we exploit the scale separation of the North Atlantic atmosphere–ocean system. This allows us to invoke a linear stochastic modeling approach to extract system information in ways that are perhaps more efficient than other, more traditional, statistical approaches. The main idea is based on the “principal oscillation pattern” (POP) approach of von Storch et al. (1988) and the “linear inverse modeling” (LIM) approach of Penland and Sardeshmukh (1995), with the latter being an extension of the former. Our results show that when taking into account the oscillatory behaviors of the AMOC and gyre anomalies and the phase lag between subpolar and subtropical latitudes, resulting from the slow equatorward propagation of circulation anomalies, there is no contradiction between Zhang’s (2008) finding and the traditional expectation of an in-phase relation between the subpolar gyre and the AMOC. We also envisage gyre-specific variations in the AMOC based on a dipole pattern revealed by our analysis. Perhaps more important, our analysis confirms the ideas of Gerdes and Köberle (1995) that the slow advection of density anomalies accounts for the slow oceanic adjustment to high-latitude buoyancy forcing. Last, it is our hope that the method shown here will be applied to other GCMs that exhibit similar interdecadal AMOC variations. This would help investigating possible agreements and disagreements between models.

The paper is organized as follows. In section 2, we briefly introduce the LIM and our other methods. In section 3, we show prediction skill of the LIM constructed in various ways and discuss the implications. In section 4, we focus on one dominant interdecadal eigenmode and diagnose the oscillation mechanisms including the involved positive and negative feedbacks. In section 5, we summarize the North Atlantic atmosphere–ocean system using spectral analysis. Our conclusions are presented in section 6.

## 2. Methodology and data

### a. Methodology

Instead of directly linearizing the full GCM equations (i.e., the tangent linear model approach), we fit a linear model to the North Atlantic atmosphere–ocean system; that is, we employ the LIM method. If fitted properly, as we shall show below, LIM not only yields the best linear description of the full dynamics, but also gives insight into the statistical properties of the nonlinear part. The essence of LIM is based on scale separation: if nonlinearities decorrelate quickly enough, their effect on slow processes of interest can be approximated as stochastic noise according to the central limit theorem. If the slow processes are mostly linear, then the system can be well described by (Penland and Sardeshmukh 1995)

where **X** is the state variable, including all relevant variables we are interested in, and ** ξ** is the noise vector. It has been known for long that there indeed is a separation of time scale in the atmosphere–ocean system. For example, by taking the annual mean of our data, almost all atmospheric processes and most oceanic eddy processes are relegated to the stochastic noise part

**, while the linearizable advection would be left in the linear dynamics part described by .**

*ξ*The operator matrix can be determined from covariance statistics as

where (τ) is the covariance matrix for lag *τ*. The covariance matrix of the LIM system would then evolve as

which, compared with the covariance matrix for lag *τ* from observations or GCMs, could be used as a test for the LIM. The quantity is the so-called Green function. The noise covariance matrix = 〈*ξξ*^{T}〉*dt* can be determined from the fluctuation–dissipation relationship

The first several eigenvectors of represent the dominant noise patterns, referred to as the noise empirical orthogonal functions (noise-EOFs). Note that the noise term ** ξ** is white in time but not in space as fast processes can have large spatial scales.

If the LIM is consistent with the original system, one can hope that the dominant linear dynamics and the noise statistics approximate the slow and fast processes of the system, respectively. The degree of consistency between the LIM and the full system can be measured by the cross-validated forecast skill, following Newman et al. (2011).

The dynamics of the linear system can be diagnosed by applying eigen-decomposition to the operator , whose eigenvectors **u**_{i} are principal oscillation patterns (von Storch et al. 1988) or empirical normal modes (hereinafter eigenmodes or simply modes; Penland and Sardeshmukh 1995), while corresponding eigenvalues *λ* = *ζ* + *ηi* give exponential decaying times (−1/*ζ*) and characteristic periods (2π/η). If no eigenvalues are degenerate, then **X** can be uniquely represented as the linear combination of the **u**_{i} values; that is, eigenmodes act as the basis of the linear system:

where the *i*th column of is the eigenmode **u**_{i} while the *i*th component of **d**(*t*) contains the corresponding coefficients. Modes that have a longer decay time scale will dominate over modes that decay quickly if they are equally excited by noise, and manifest themselves in the summation of all modes, namely the **X** time series. This is why less damped eigenmodes attract more attention. The quantity **d**(*t*) can then be obtained by projecting **X** onto ^{−1} or equivalently, the adjoint patterns , which form a biorthogonal set with .

To diagnose related fields **Y** that are not included in constructing the LIM, we calculate the so-called associated correlation pattern *υ*_{i} as the regression of **Y**(*t*) onto *d*_{i}(*t*) (von Storch et al. 1988),

where the error **e**(*t*) is to be minimized.

Besides LIM, we also examine the spectral properties of two kinds of oscillators: a linearly driven damped oscillator and a delayed oscillator, which are the two main conceptual prototypes for atmosphere–ocean coupling (Griffies and Tziperman 1995; Battisti and Hirst 1989). The linearly driven oscillator *x* can be described by

where *f* is stochastic white forcing and *F*(ω) its Fourier coefficients, *γ* is the damping coefficient, and ω_{0} is the natural frequency of *x*. The phase lag *φ* between the oscillator *x* and the force *f* is then (Taylor 2005)

and the amplitude of the oscillator *A* varies as (Taylor 2005)

Coefficients such as *γ* and ω_{0} in this simple linear system can be easily determined by combining the LIM technique and Hilbert transform (see appendix A).

As to the delayed oscillator, we adopt the separated form as in Sun et al. (2015):

where *x* denotes SST or other oceanic variable, denotes atmospheric forcing, *c* is the heat capacity, *a* is the forcing strength, *b* is *x*’s time scale, and *τ* is the delay time between *x* and .

Last, all spectral results are derived from using the multitaper method (MTM; Mann and Lees 1996) and further smoothed in frequency domain. Significance levels for coherency are determined from an *F* test while that for a spectrum are determined from a red noise null hypothesis and a χ^{2} test.

### b. Data

In this study, we use a 4000-yr-long GFDL CM2.1 preindustrial control run. The CM2.1 atmosphere component has a resolution of 2° × 2.5° with 24 vertical levels. The ocean component has a horizontal resolution of 1° × 1°, which becomes progressively finer approaching the equator, where the meridional resolution is ⅓°. The ocean has 50 vertical levels in total, 22 of which are evenly spaced within the top 220 m. There is a moderate drift in the first five centuries, evident from Delworth et al. (2006) and also from Delworth and Zeng (2012). However, our conclusions about the least damped oscillatory mode are not noticeably influenced by the trend. Hence, we do not detrend the data. All data used in this study are annual mean anomalies.

#### LIM setup

Our state variable **X** includes the SST, barotropic streamfunction (hereafter gyre), meridional streamfunction (hereafter AMOC), and surface salinity (SSS) defined for the region bounded by 25°N, 65°N, 80°W, and 20°E, with Hudson Bay and the Mediterranean Sea masked. SST and SSS are included to show possible atmospheric influences. We do not explicitly include atmospheric variables into **X** as they mostly behave like stochastic noise on the annual time scale and can therefore deteriorate . SST is defined here as surface temperature above sea ice when sea ice is present, and adopting the temperature below sea ice does not change our conclusions. The two streamfunctions are included to represent the oceanic circulation: while the AMOC is known to be important for the North Atlantic atmosphere–ocean system (Delworth and Mann 2000), the horizontal gyre circulation is closely linked with AMOC due to the east–west density gradient in the subpolar gyre region [see Kwon and Frankignoul (2014) for an elegant summary] and hence are also included. While some studies (Zhang 2010; Kwon and Frankignoul 2014) calculate AMOC streamfunction in the density space rather than in the depth space, this approach mostly increases the AMOC magnitude at high latitudes in GFDL CM2.1, which does not change our main conclusions. Therefore, we define AMOC in the depth space here.

The training lag to compute from Eq. (2) is 1 yr. The choice of 1 yr is motivated by the fact that the involving lagged covariance calculations must be performed at lags smaller than half the smallest oscillation period of the system to give accurate results.^{1}

To exclude possible nonlinear, small-scale processes, SST and SSS data and *u*, υ data for the gyre streamfunction computation are first interpolated to a 4° × 4° grid. Velocity data *υ* used to compute the AMOC streamfunction are not spatially smoothed. Because various fields with different units are involved, each anomaly field is normalized by its domain-averaged climatological root-mean-square amplitude; that is, our results are all nondimensional.

To reduce the dimension of the system, we project data onto the empirical orthogonal function (EOF) space and truncate to retain only the leading principal components (PCs). Hence **X** includes the first several PCs of SST, SSS, gyre, and AMOC. To facilitate interpretation, we present results transformed back from EOF to grid space. We retain 9, 8, 9, and 9 EOFs for SST, gyre, AMOC, and SSS, which explain 85%, 83%, 97%, and 73% of total variance, respectively. The numbers of retained PCs are chosen by maximizing LIM’s SST prediction skill at 5-yr lead time. We choose a 5-yr lead time simply because we expect the least damped oscillatory mode to influence the prediction skill significantly while 5 yr is about a quarter of the oscillatory mode period (see below) where a sinusoidal oscillation starts to diverge significantly from a linear fit line. Moderate variations in adopted PC numbers do not change our main conclusions.

To illustrate ocean dynamics, we calculate associated correlation patterns for the mixed layer depth (MLD) and for the upper ocean heat content, the upper ocean in situ density, and the upper ocean *υ* velocity over the top 1000 m. In examining the role of the atmosphere, we focus on the North Atlantic Oscillation (NAO), defined as the first EOF of annual mean sea level pressure over the region bounded by 20°N, 80°N, 90°W, and 40°E. We also use an AMOC strength index, defined as the maximum meridional streamfunction at 45°N (denoted as 45°AMOC).

## 3. How well does LIM represent GCM dynamics?

As SST receives influence from both above and below, SST prediction skill could serve as a critical test for approximations of the North Atlantic atmosphere–ocean system. We therefore examine the SST prediction skill to test how well LIM represents the system. In Fig. 1 we compare the LIM prediction skill of North Atlantic SST against an autoregressive (AR-1) model. The increase in skill of LIM relative to AR-1 is concentrated in the subpolar gyre region, implying that different mechanisms control subpolar versus subtropical SST variations. We also note that the increase in skill does not vary monotonically with lag but reaches a maximum at around 5-yr lead time (Fig. 1c), which suggests that it is the inclusion of relatively low-frequency information, probably ocean dynamics, that contributes to the skill improvement. We explore this conjecture by comparing with a different LIM, constructed from SST only (referred to as SST-LIM). The fact that the oscillation in lagged autocovariance of North Atlantic SST is better captured by LIM than by SST-LIM (Fig. 1f) confirms that including ocean dynamics, contained in the two streamfunctions of **X**, allows LIM to capture the interdecadal variability in SST and make better long-term predictions.

Zanna (2012) employed HadSST2 to examine the prediction skill by a SST-LIM and found that the skill is worst in the subpolar gyre region and much better in the subtropical North Atlantic. We repeat the calculation of Zanna (2012) with our CM2.1 data and find that the skill still concentrates north of 30°N. This discrepancy indicates that LIM’s skill is state-dependent; our control simulation has very different dominant EOFs than the ones derived from HadSST2, probably due to the external forcing in observations.

Next, we focus on investigating the oceanic dynamics that are responsible for the interdecadal SST variations. As implied by Fig. 1d, SST shows maximum variance in the region bounded by 40°N, 65°N, 80°W, and 30°W. Later we compute SST averaged over this region between GCM and LIM and denote this single variable as our SST index.

## 4. Least damped oscillatory eigenmode

There are a total of 35 eigenmodes generated from this LIM, 14 of which are complex while the remaining 7 are real (2 × 14 + 7 = 35). As eigenmodes summed together must equal **X**, which is real, complex eigenmodes occur as conjugate pairs. Counting each complex conjugate pair as a single mode, we get 21 real/complex modes, the first four of which are listed in Table 1. While complex eigenmodes represent oscillatory modes, real ones represent exponentially decaying modes. The least damped mode has a decay time scale of 26 yr, which is very different from the extremely long decay time scale of 1200 yr in Tziperman et al. (2008). This difference suggests that the 3D temperature and salinity fields studied in Tziperman et al. (2008) contain more information than the gyre and AMOC streamfunctions do. Table 1 also lists contributions of the first four eigenmodes to the SST index variations and 45°AMOC variations. Note that the first four modes all contribute much to the SST index and 45°AMOC variations.

We find that the optimal initial condition mostly projects onto the adjoint patterns of the first four eigenmodes, while the subsequent evolution mostly projects onto the first four eigenmodes. The fact that eigenmodes are not orthogonal to each other confirms that the system is nonnormal and that interference among eigenmodes with different decay time scales is important (Tziperman et al. 2008). However, as we demonstrate below, recognizing the importance of nonnormality does not exclude a particular eigenmode from describing important physical mechanisms.

Among these 21 eigenmodes, we are most interested in eigenmode 2 (hereafter mode 2), not only because it is the least damped oscillatory mode but also because it contributes remarkably to both the SST index variations and the 45°AMOC variations (with correlation of 0.45 and 0.73 respectively; see Fig. 2), revealing the connection between the deep ocean and low-frequency SST variations, and possibly implying atmosphere–ocean interactions. The associated eigenvalue suggests that mode 2 has an estimated period of around 20 yr whereas spectral analysis suggests a probable range of 15–30 yr.

### a. Robustness test

As this section focuses on interpreting the spatial patterns of mode 2, it is necessary to first show that its spatial patterns are robust. We start with testing the robustness against changing the number of PCs included in **X**. We choose the PC number combinations that maximize SST explained prediction variance for a chosen SST PC number (ranging from 2 to 10) by varying gyre, AMOC, and SSS PC numbers. Then we calculate spatial correlations between eigenmodes in those combinations and eigenmodes in the [9, 8, 9, 9] combination. It turns out that for the least damped oscillatory mode, all spatial pattern correlations are greater than 0.95 and are highly significant (not shown), indicating that mode 2 is robust. The decay time scale, period, and contributions to the SST index and 45°AMOC variations of mode 2 also barely change. We also test the robustness of eigenmode patterns to sampling. When we evenly segment the whole 4000-yr data into four segments of 1000 yr each, we find that mode 2 in each segment is nearly identical to that in the other segments (not shown).

### b. Physical interpretation

The spatial pattern of mode 2 is shown in Fig. 3, where we see that mode 2 evolves from 0° to 270° as ** α** →

**→ −**

*β***→ −**

*α***, where mode 2 (**

*β***u**and λ = ζ +

*ni*are the eigenvector and eigenvalue of the complex mode; see section 2a for details). As the complex mode is only unique up to an arbitrary rotation, we choose the normalization parameter

*c*by requiring that

*α*^{T}

**= 0,**

*β*

*α*^{T}

**= 1, and**

*α*

*β*^{T}

**≥ 1 (Penland and Sardeshmukh 1995). That is, we require**

*β***and**

*α***to be 90° apart (i.e., in quadrature).**

*β*Most notable is the AMOC pattern, revealing a slow equatorward propagation of the high-latitude AMOC anomalies over a 10-yr period (note the progression from Fig. 3c to Fig. 3g to Fig. 3k). As the positive AMOC anomalies strengthens and propagates to subtropical latitudes, a negative anomaly forms at high latitudes (Fig. 3k), which then initiates the other half cycle of the opposite sign. This slow propagation between 60° and 35°N suggests that velocity anomalies propagate with density anomalies at the speed of mean advection. In the gyre component of mode 2 one can see a similar equatorward propagation: accompanying the high-latitude positive AMOC anomaly, there is a cyclonic gyre anomaly south of Greenland (cf. Figs. 3b,c), which strengthens over time and propagates to the gyre boundary after 5 yr (Fig. 3f) and then to subtropical latitudes another 5 yr later, when the Gulf Stream shifts to the south in link with the gyre anomaly (Fig. 3j). Simultaneous with the phase flip in AMOC anomalies, there are also anticyclonic gyre anomalies in high latitudes at phase 180°.

A similar slow, equatorward propagation of circulation anomalies was revealed in an early GFDL ocean primitive equation model by Gerdes and Köberle (1995), who studied the oceanic adjustment to a prescribed Iceland Sea buoyancy forcing over 10 yr of simulation. The slow adjustment was shown to result from the density anomalies being advected equatorward slowly along the North American coast, which induces velocity anomalies within the North Atlantic Deep Water (NADW) flow (confirmed by a model tracer propagation). A similar result was reproduced in the GFDL CM2.1 model more recently (Zhang 2010). That the circulation adjusts with the advection of density anomalies explains the close correspondence between the gyre anomalies and the AMOC anomalies: they both incorporate NADW flow anomalies, which responds to the deep density anomalies. Note that the slow equatorward propagation of NADW flow anomalies can be readily confirmed here in its associated correlation pattern (not shown). It was shown that in the absence of background velocity, the anomalous circulation induced by the density anomalies can on its own advect the density anomalies southward, a mechanism termed “self-advection by density perturbations” (Gerdes and Köberle 1995). However, in the GCM or real-world conditions, the mean velocity is typically much larger than the anomalous velocity and hence ultimately determines the propagation speed of density and circulation anomalies.

We would like to point out here that the in-phase relation between the AMOC and the subpolar gyre via NADW flow anomalies is consistent with other model studies of oceanic response to enhanced Labrador Sea convection (Böning et al. 2006; Eden and Willebrand 2001). The physical significance of mode 2 is therefore established by its similarity with GCM perturbation experiments and is further confirmed by its similarity with patterns of the gyre and the AMOC streamfunction when regressed onto the 45°AMOC index (not shown). However, none of these earlier studies realized that the slow propagation is part of an interdecadal oscillation, whose mechanism, as we shall show below, is mostly embedded in the subpolar gyre circulation.

To clarify the oscillation mechanism, we start from explaining what processes are responsible for the phase flip (a negative feedback) and the growth of anomalies (a positive feedback), two essential processes for any oscillations. To illustrate the negative feedback that accounts for the phase flip we calculate the associated correlation patterns [Eq. (6)] for upper ocean heat content, upper ocean density, upper ocean υ velocity, and mixed layer depth (Fig. 4). The anticyclonic gyre anomaly in the eastern subpolar gyre at 90° phase (Fig. 3f) clearly corresponds to the collocated positive ocean heat content anomaly (Fig. 4f) and the negative density anomaly (Fig. 4g), which then propagate around the subpolar gyre and toward the gyre boundary (note the progressions from Fig. 4f through Fig. 4j to Fig. 4n, and from Fig. 4g through Fig. 4k to Fig. 4o). Comparing the upper ocean heat content anomaly and the upper ocean υ velocity anomaly at 90° phase (Figs. 4e,f), the heat content anomaly seems to result from anomalous poleward advection due to the enhancement/shift of the North Atlantic Current. Kwon and Frankignoul (2014) also identified the density anomalies in the eastern subpolar gyre as the critical factor for the phase flip in the AMOC ~20-yr periodicity in their CCSM3 model.

To better understand the heat content/density anomalies, we modify our state vector **X** such that it explicitly includes upper ocean heat content as an evolving component and we inspect the simultaneous evolution of gyre anomalies and heat content anomalies of the same eigenmode in Fig. 5. We find that the least damped mode in this newly constructed system is essentially the same as the one discussed above, suggesting that ocean heat content is so strongly coupled with the gyre or AMOC anomalies that it does not provide linearly independent information about the mode that interests us. The evolution starts from a specific initial condition (explained below) and then evolves under the Green function (τ). Note that there is no noise input during this evolution; that is, . As we want the evolution to heavily project onto mode 2, we design the initial condition as follows: we first calculate the optimal initial condition, defined as the first right singular vector of the Green function (τ_{0}), which corresponds to the largest singular value (Penland and Sardeshmukh 1995), and then regress the optimal initial condition onto all adjoint patterns and only use the part that regresses onto the adjoint pattern of mode 2, which highly correlates with the optimal initial condition (*r* = 0.95). For the same reason, we pick τ_{0} to be 7 yr, around a quarter of the mode 2 period (the optimal initial condition is almost unchanged for τ_{0} ranging from 4 to 10 yr). Furthermore, we only retain the gyre and AMOC components of the initial condition because including other components does not change the evolution much. It is clear that in Fig. 5 gyre anomalies and upper ocean heat content anomalies closely parallel with each other and the phase change happens before year 13.

To further demonstrate the importance of the heat content anomaly for the phase change, we modify to eliminate the influence from ocean heat content onto other variables. In Fig. 6, we compare the evolution of gyre and AMOC anomalies with ocean heat content influence turned on or off. Without the influence of ocean heat content, the phase flip disappears and the circulation anomalies decay without much propagation.

The positive feedback that is responsible for anomaly growth seems to be consistent with the temperature feedback proposed by Levermann and Born (2007); that is, temperature anomalies induce circulation anomalies via geostrophic adjustment, which then enhance the initial temperature anomalies via anomalous advection, forming a positive feedback. The mutual enhancement between ocean heat content anomalies and gyre anomalies is evident from Fig. 5: the initial heat content anomaly associated with the North Atlantic Current anomaly at year 1 induces an anticyclonic gyre anomaly at the eastern subpolar gyre region at year 4, whose northward branch then enhances the heat content anomaly; this mutual enhancement operates until year 10 when the heat content anomaly is advected to the subpolar gyre center by the mean gyre circulation. Note that the positive feedback stated here accounts for the growth of anomalies seen in Fig. 3f and 3j while the negative feedback described above accounts for the phase flip seen in Fig. 3f or 3n. The fact that the gyre and AMOC anomaly evolves much weaker during the evolution without heat content influence (Fig. 6) supports the view that the heat content anomalies are also part of the positive feedback. Associated with the enhancement and propagation of heat content anomalies around the subpolar gyre, convection in the Labrador Sea gradually enhances (cf. Figs. 4d and 4h). This explains the growth of AMOC anomalies at high latitudes from 0° to 90° (cf. Figs. 3c and 3g).

Last, we would like to point out that in the above discussion it is implied that density anomalies in the subpolar gyre region are dominated by temperature, while salinity anomalies partially counteract it (not shown but can be inferred from the opposing sign of SSS and SST in Fig. 3). In summary, anomalous advections by North Atlantic Current give rise to temperature anomalies in the eastern subpolar gyre region, which then grow via mutual enhancement between temperature anomalies and circulation anomalies (i.e., the positive feedback discussed above); meanwhile, the temperature anomalies are advected by the mean circulation (this advection could be aided by the anomalous circulation but the mean circulation should dominate, as discussed above) from the eastern subpolar gyre region toward the gyre center region; by the time that temperature anomalies reach the gyre center region, the related gyre anomalies induce temperature anomalies of the opposite sign in the eastern subpolar gyre region and the other half of the cycle then begins (i.e., the negative feedback discussed above). The ~10-yr half period is thus set by the mean advection of heat content (density) anomalies around the subpolar gyre, starting from the eastern subpolar gyre and ending at the subpolar gyre center. As the heat content anomalies reach the convection site (Fig. 4, black contour), they induce convection anomalies, which reinforce the gyre anomalies and initiate the AMOC anomalies (Figs. 3b,c), which then propagate equatorward along with the density anomalies.

Meanwhile, other processes probably help with enhancing or maintaining the temperature or circulation anomaly against dissipation during the propagation. For example, the bottom vortex stretching around the Grand Banks due to downslope motion (Zhang and Vallis 2007) seems to intensify the cyclonic gyre anomalies (Fig. 3f); shifts of the Gulf Stream should intensify the temperature anomalies and hence the gyre anomalies (Fig. 3j).

The above discussions suggest that mode 2 is an ocean-only mode. However, atmospheric forcing could be important by providing external energy. We also notice that mode 2 is associated with a strong SST anomaly in the Labrador Sea (Fig. 3a) accompanying the high-latitude AMOC anomaly (Fig. 3c), which suggests an influence of the atmosphere. Because of the difficulty of including atmospheric variables into LIM, we cannot directly examine cause and effect and exclude the possibility of a coupled atmosphere–ocean oscillator via LIM. Therefore, in the next section, we use the spectral analysis to examine the role of the atmosphere, represented by the NAO.

### c. References to previous work

First, we note that the oscillation between the monopole and dipole patterns of upper ocean heat content (Figs. 4b,f) is similar to that described in Sévellec and Fedorov (2013). They summarize the oscillation mechanism as “geostrophic self-advection” in combination with the mean meridional temperature gradient, similar to the Rossby wave propagation mechanism except that density here plays the role of potential vorticity. Their geostrophic self-advection is similar to the aforementioned “self-advection by density perturbations” in Gerdes and Köberle (1995). However, as argued above, with the density anomalies propagating along the cyclonic mean circulation (Fig. 5), the mean advection should dominate over anomalous advection if the mean velocity is relatively large. In light of this, different conclusions of these two studies seem to result from different mean strengths of the subpolar gyre: with a too weak subpolar gyre in Sévellec and Fedorov (2013) (see their Fig. 1), the anomalous advection dominates. Nonetheless, both studies attribute the generation of density anomalies to anomalous advection. In addition, the phase correspondence between density and AMOC anomalies is similar between the two studies: dipole density anomalies in the subpolar gyre region correspond to an enhanced AMOC in the middle to high latitudes. Therefore, these two studies seem to reveal similar mechanisms but differ in whether mean advection or anomalous advection dominates, resulting from different strengths of the mean circulation.

Zhang (2008) noticed that in the same GCM as this study, AMOC anomalies at 40°N closely correspond to a north–south dipole in gyre anomalies and 400-m subsurface temperature anomalies, with the weakening and warming subpolar gyre corresponding to a strengthening AMOC. This is consistent with our mode 2 at phase 180° (Figs. 3j,k and 4j). Meanwhile, if one defines the AMOC by subpolar latitudes, the 180° phase of mode 2 then yields the traditional in-phase relation between the AMOC and the subpolar gyre. The contradiction hence stems from a dipole pattern in AMOC anomalies. This dipole pattern in AMOC at phase 0° and 180° is also reminiscent of the observed opposing changes in AMOC in subpolar versus subtropical gyres (Lozier et al. 2010), which implies that the AMOC changes seen in hydrographic data could result from internal variations rather than anthropogenic trends. In addition, our mode 2 reproduces the relation between AMOC strength and shifts of the Gulf Stream found in earlier studies of the same GCM (Zhang and Vallis 2007; Joyce and Zhang 2010). The associated surface signals—that is, the cooling (freshening) confined in the Northern Recirculation Gyre region and the simultaneous appearance of warm (saline) anomalies south of Greenland (Zhang and Vallis 2006)—are also captured by mode 2 (Figs. 3i,l).

The heat content evolution here is similar to that in the preindustrial run of a related model, GFDL ESM2M [see Fig. S3 in MacMartin et al. (2016)] and not dissimilar to their 4 × CO_{2} run except for the latitudinal extent. While MacMartin et al. (2016) point out that the strong coupling between the subpolar gyre and AMOC is important in their 4 × CO_{2} run, we argue that a similar mechanism acts in our CM2.1 preindustrial run, and possibly also in their preindustrial run. The contradiction that different gyre strength and hence different advection time scales in different CO_{2} scenarios yield similar oscillation period could probably be reconciled by considering that a stronger subpolar gyre also has a larger extent and thus possesses a similar advection time scale as a weaker one.

Yoshimori et al. (2010) employed regression analysis to explore the oscillatory behavior of the AMOC in CCSM3 and also found the southward propagation of AMOC anomalies. However, there are two main differences between their and our conclusions. First, they found the density anomalies that lead to anomalous deep convection are mostly due to salinity anomalies that are advected from the Labrador Sea toward the Irminger Sea by the anomalous gyre circulation. This is very different from our finding that the mean advection of temperature dominated density anomalies induces the anomalous convection in Labrador Sea and an anomalous AMOC. This discrepancy may result from different climate states: their CCSM3 run is much colder and fresher than our CM2.1 preindustrial control run, which thus causes a higher sensitivity to salinity. Second, they found that atmosphere–ocean coupling is what triggers the phase change of the oscillation, different from our ocean-only mode conclusion. The comparison between Yoshimori et al. (2010) and our study highlights the advantage of LIM over regression analysis: while regression analysis can only reveal relationships between pairs of variables, LIM can extract a dynamically consistent system of many variables, which then can be easily extended to other GCMs.

Now we return to the nonnormality effect. Tziperman et al. (2008) found the importance of modal interference in explaining thermohaline circulation anomaly growth and concluded that the thermohaline circulation or AMOC cannot be interpreted as a simple harmonic oscillator. However, as noted above, the eigenmodes extracted in our study are very different from those extracted in their study. Our mode 2 is the second least damped mode and hence gains more importance than the interdecadal oscillatory mode in their study. The optimal growth in our case starts with cancellation among multiple eigenmodes, as in their study, but ends with a pattern heavily projecting onto mode 2, which is not the case in their study. Therefore, the contrasting conclusions based on the same GCM are due to the fact that the two studies use different variables to construct the LIM and hence study different subsystems of the North Atlantic atmosphere–ocean system.

## 5. Spectral analysis and the NAO–AMOC relationship

In this section, we employ a totally independent method, spectral analysis, to provide more insights into the nature of the oscillation in the North Atlantic system. Figure 7 compares the spectral and correlation characteristics between 45°AMOC and NAO with that between a linear oscillator and its white forcing. Besides the full 45°AMOC time series (Fig. 7, left column), we also look at the mode 2 contribution to 45°AMOC (Fig. 7, right column). As shown in Fig. 7a, at decadal and longer frequencies, the coherency between NAO and 45°AMOC is statistically significant and the phase lag between 45°AMOC and NAO (labeled “phase-model”) agrees with that of a linear oscillator *x* driven by white forcing *f* [labeled “phase-spring”; see Eq. (8); for the determination of parameters see appendix A]. In particular, for the linear driven oscillator *x* the phase lag ϕ(ω) at the natural frequency ω_{0} should be −90°, which agrees well with the AMOC–NAO −90° phase lag at the AMOC spectral peak (~15 yr). The phase based on a delayed oscillator mode with memory time scale of 10 yr (labeled “phase-delay”) is also shown in Fig. 7a (for discussion of other AMOC memory time scales, see appendix B) and clearly show less agreement with the phase calculated from GCM. In light of the linear oscillator theory, the spectrum of AMOC (Fig. 7b) looks quite similar to the amplitude response of the oscillator *x* to white forcing *f*, which reaches its maximum near the natural frequency ω_{0}. The lag correlation between AMOC and NAO (labeled “cor-model” in Fig. 7c) is also very similar to that between a linear oscillator *x* and the white forcing *f* (labeled “cor-spring” in Fig. 7c). The maximum correlation is obtained when the NAO leads the AMOC by approximately one-quarter period (−4 yr), consistent with a −90° phase lag at resonance. Thus, we are encouraged to interpret the AMOC response to the NAO as that of a linear damped oscillator driven by white forcing rather than a delayed oscillator. Last, we emphasize that the spectral results between NAO and 45°AMOC that from mode 2 alone (Figs. 7d–f) agrees with above interpretation, namely that the phase, coherency, spectrum, and correlation between the NAO and the 45°AMOC from mode 2 reproduce well that between the NAO and the entire 45°AMOC at interdecadal frequencies. This supports the idea that the ~20-yr period oscillation seen in the North Atlantic climate system in this GCM mostly arises from mode 2. This reproduction is also a very encouraging support for LIM as information about the NAO is only indirectly included in the LIM through the fluctuation–dissipation relationship [Eq. (4)].

We emphasize that the agreement between the NAO–AMOC relation and the linearly driven damped oscillator indicates that the NAO serves mostly as an “external” force to the AMOC oscillator instead of being an internal component of the oscillator; that is, the AMOC mode is more likely an ocean-only mode rather than an atmosphere–ocean coupled mode. This prompts the question of why a white-noise NAO can excite this interdecadal oscillation. Results here indicate that the answer is resonance. The AMOC would respond to the NAO on all frequencies but the response is much higher near the AMOC’s natural frequency, which thus exhibits itself in the AMOC time series and the AMOC–NAO lag correlation. Delworth and Zeng (2016) also examined the “resonance” idea between NAO and AMOC in the same GCM by adding extra NAO-related heat flux forcing of various frequencies and showing that the AMOC response amplitude does maximize when the external forcing is around a period of 20 yr. However, they did not identify mode 2 as an eigenmode of the system. Also, they could not identify the NAO as “required” to maintain the system, which we shall demonstrate below with the noise-EOFs.

We regress SST, gyre streamfunction, AMOC streamfunction, and SSS onto the NAO to show direct evidence of a NAO influence. In Fig. 8, the evolution from lag = 1 yr to lag = 13 yr bears a similarity to mode 2 from phase 0° to 270°. However, the lag = 0 pattern lacks obvious connection with the patterns at subsequent lags and it is not similar to the mode 2 adjoint pattern, that is, the optimal excitation pattern of mode 2 (see Fig. 6, I.C. column). The SST pattern shown in Fig. 8a is clearly part of the well-known “tripolar pattern” (map begins at 25° so only two poles are visible) that results from the NAO and related turbulent heat flux (Cayan 1992a,b), while Fig. 8b resembles the intergyre gyre resulting from the NAO-induced wind stress curl, described by Marshall et al. (2001) and also seen in Eden and Willebrand (2001). The dipole AMOC anomalies in Fig. 8c are also driven by NAO-related surface wind stress via Ekman drift (Eden and Willebrand 2001). The quasi-tripolar pattern in Fig. 8d is the salinity response to NAO-induced Ekman advection (Mignot and Frankignoul 2003). In short, the lag = 0 pattern is due to the immediate response of SST, SSS, and streamfunctions to NAO-induced wind stress and heat flux anomalies. However, this direct influence is transient while the influence that projects onto the deep convection is long-lasting. In other words, the mechanism proposed by Marshall et al. (2001) seems not active in this GFDL CM2.1 mode, in which the NAO coupled with the wind stress–driven intergyre gyre pattern yields an interdecadal oscillation.

Now we would like to direct the reader’s attention to the dominant noise pattern, the first noise-EOF (defined as the first EOF of noise covariance ), which is remarkably similar to the NAO regression pattern at lag = 0 yr (Fig. 9). This suggests the dominant role of the NAO in driving the North Atlantic atmosphere–ocean system. Although the importance of the NAO is quite a constant in the literature (Barrier et al. 2014), the new contribution here is a self-consistent dynamical description of AMOC and gyre as an oscillating eigenmode of the system driven by rapidly varying forcing, which is “objectively” identified with the NAO.

## 6. Conclusions

Utilizing linear inverse modeling (LIM), we have identified an interdecadal oscillatory mode of the North Atlantic atmosphere–ocean system that accounts for 20% of the SST-index variations (SST averaged over region bounded by 40°N, 65°N, 80°W, and 30°W) and 53% of the 45°AMOC variations. The oscillation mechanism is mostly embedded in the subpolar gyre: density anomalies in the eastern subpolar gyre, resulting from anomalous advection of North Atlantic Current, propagate around the subpolar gyre following the mean gyre circulation; when density anomalies reach the subpolar gyre center around 10 yr later, the associated gyre anomalies induce density anomalies of the opposite sign in the eastern subpolar gyre and initiate the other half of the cycle. While the density anomaly reaches the Labrador Sea deep convection region, it brings about AMOC anomalies too. This oscillatory mode also incorporates a slow equatorward propagation of the circulation anomalies, along with density anomalies, from subpolar latitudes to subtropical latitudes, evident in both the AMOC and the gyre streamfunction. The slow equatorward propagation is consistent with an early study of a related oceanic GCM (Gerdes and Köberle 1995), which has been reproduced with GFDL CM2.1 more recently (Zhang 2010). Combining the oscillation and the equatorward propagation, it is readily seen that signals originating from subpolar latitudes can propagate into subtropics and exert influences there. Although the influence could be concealed by local wind influences (Lozier 2010), the long-term potential predictability resulting from the oscillation still argues for the monitoring of the subpolar gyre (Lozier et al. 2017).

LIM, aided by associated correlation pattern calculations, allows us to extract dynamically consistent signals in multiple variables. Our analysis reproduces many results that have been seen in other studies yet remain unconnected. Further, through spectral analysis, we have shown that this mode is more probably an ocean-only mode excited by atmospheric forcing and, in particular, that this mode can be approximately viewed as a linearly driven, underdamped oscillator with NAO providing at least part of the “external” forcing. The dominant role of NAO is confirmed by our “noise-EOF” calculation.

When regressing SST, streamfunctions, and SSS onto the NAO, we see a clear gap between patterns of zero lag and patterns of nonzero lags. We argue that this is an indication of two different processes: a direct, fast response to NAO-related wind stress and heat flux forcing is followed by a slow response to NAO-induced Labrador Sea buoyancy forcing. It is the latter process that has potential predictability owing to the existence of an oscillation and the slow evolution.

Last, we want to point out the advantages of LIM over other, more traditional statistical techniques and perhaps sometimes over GCM experiments, as well as the limits of LIM. If scale separation and the linearity of slow processes are quite well satisfied, LIM can reliably extract dominant modes of slow processes, as demonstrated by mode 2 in this study. Meanwhile, some statistical characteristics of the fast, nonlinear processes can also be estimated via noise-EOF calculations and other techniques shown in Penland and Hartten (2014). Moreover, one can easily modify to experiment with dominant processes without the limits of computation cost. Last, because of the additive noise assumption in our LIM, we cannot examine the possible feedback from the ocean onto the atmosphere, which is a goal of future work.

## Acknowledgments

Funding for B.Z. and T.R. was provided by the National Science Foundation under Grant 1446292. Partial support also came from the Office of Global Programs. We especially thank Matthew Newman and reviewers of the manuscripts for grateful suggestions that help to improve the work significantly. We appreciate Tara McQueen for her help with acquiring data.

### APPENDIX A

#### Determination of *γ* and ω_{0}

Consider , where *ζ* is complex. Then *x* = ℜ(ζ) and *x*_{H} =(ζ), where *x*_{H} is the Hilbert transform of *x*.

Let *y* = *ζ* ′, then . Thus,

The linear operator has eigenvalues if , where . Then the eigenvector **u**_{1} corresponding to is

where is a normalization constant. Similarly, the other eigenvector is

Therefore,

and, absorbing the normalization constants into and ,

The preceding yields

The above operator can be determined by employing LIM and then used to derive *γ* and ω_{0}. For example, if *x* is 45°AMOC from GFDL CM2.1, LIM solves

hich yields γ = 4.8 × 10^{−9} s^{−1}, ω_{0} = 1.26 × 10^{−8} or a decay time of 7 yr and a period of 16 yr. If *x* is 45°AMOC from mode 2, LIM yields a decay time of 15 yr and a period of 18 yr. Considering the uncertainty, the internal period decided here is consistent with the period estimated from mode 2 (see text).

### APPENDIX B

#### Phase of Delayed Oscillator

Consider the phase lag between *x* and as in

Write the differential equation in finite difference form,

where Δ*t* is the time step and superscript denotes the time in the sense that *t* = *n*Δ*t*. Take the Fourier transform of the above equation and denote the Fourier transform of *x* as (*x*) and complex conjugate of *x* as *x*^{★}:

Denote cross-spectrum as :

Therefore, the phase lag follows:

Thus, (ω) is a function of only *cb*/Δ*t*. Note that if both the real and imaginary part of are negative, should be in the range of [−π/2, −π]. Note that *cb*/Δ*t* is actually the decay time of *x* in units of Δ*t*. When *x* denotes an oceanic variable, *cb*/Δ*t* should be much greater than 1. Hence,

Figure B1 shows the plot of (ω) for *cb*/Δ*t* equal to 5, 10, 20, 50, and ∞.

## REFERENCES

_{2}

*Classical Mechanics*. University Science Books, 786 pp.

## Footnotes

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

^{1}

This is analogous to Fourier analysis, which only resolves frequencies smaller than the Nyquist frequency. To minimize this “Nyquist lag” issue, one must often choose the smoothing window as the training lag, which in our case is one year.