## Abstract

The response of the North Atlantic meridional overturning circulation (MOC) to wind stress forcing is investigated from an observational standpoint, using four time series of overturning transports below and relative to 1000 m, overlapping by 3.6 yr. These time series are derived from four mooring arrays located on the western boundary of the North Atlantic: the RAPID Western Atlantic Variability Experiment (WAVE) array (42.5°N), the Woods Hole Oceanographic Institution Line W array (39°N), RAPID–MOC/MOCHA (26.5°N), and the Meridional Overturning Variability Experiment (MOVE) array (16°N). Using modal decompositions of the analytic cross-correlation between transports and wind stress, the basin-scale wind stress is shown to significantly drive the MOC coherently at four latitudes, on the time scales available for this study. The dominant mode of covariance is interpreted as rapid barotropic oceanic adjustments to wind stress forcing, eventually forming two counterrotating Ekman overturning cells centered on the tropics and subtropical gyre. A second mode of covariance appears related to patterns of wind stress and wind stress curl associated with the North Atlantic Oscillation, spinning anomalous horizontal circulations that likely interact with topography to form overturning cells.

## 1. Introduction

The Atlantic meridional overturning circulation (MOC) is the primary driver of poleward heat transport by the ocean. At subtropical latitudes, it is responsible for about 70% of the poleward ocean heat transport and 25% of the combined ocean and atmosphere poleward heat transport (Ganachaud and Wunsch 2000). Numerical models suggest that over the twenty-first century, the MOC will reduce in strength (Vellinga and Woods 2002) with associated reduction in the northward heat transport (Johns et al. 2011). Our ability to properly simulate, or accurately observe, a climatic trend in MOC records is impaired by our incomplete understanding of the origins of MOC variability.

The MOC in numerical models varies on a broad range of time scales, from decadal scales (Delworth et al. 1993, 2012) to interannual scales (Biastoch et al. 2008; Köhl and Stammer 2008; Zhao and Johns 2014a) and to annual (seasonal) and shorter scales (Hirschi et al. 2007; Blaker et al. 2012; Zhao and Johns 2014b). At first, processes on different time scales could be expected to linearly superpose, but numerical simulations suggest that intrinsic interannual variability of the MOC can spontaneously appear under climatological atmospheric forcing (Grégorio et al. 2015). A decade of continuous observations has confirmed that the Atlantic MOC at 26°N exhibits broadband variability (McCarthy et al. 2015), with amplitudes larger than anticipated (Srokosz and Bryden 2015). As an example, the Atlantic MOC has shown an exceptional downward linear trend of about 0.5 Sv yr^{−1} (1 Sv ≡ 10^{6} m^{3} s^{−1}) (Smeed et al. 2014), in addition to interannual variations including a year-long dramatic reduction of about 30% (McCarthy et al. 2012). At the annual time scale, the MOC at 26°N shows a substantial seasonal cycle of roughly 30% of its absolute magnitude. Prior to the 26°N moored sustained observations, the Atlantic MOC had been estimated from synoptic hydrographic surveys. From five surveys spanning 50 years, a reduction of 8 Sv was identified (Bryden et al. 2005), but this was later mostly attributed to aliasing of the seasonal variability of the MOC into longer time scales (Kanzow et al. 2010). Thus, analysis of the MOC variability is complicated by the superposition of multiple time scales of variability.

At any given latitude, the observed and simulated variability of the MOC may be induced by local or remote forcing. For example, the seasonal cycle of the MOC at 26°N is explained by coastal wind forcing off the Canary Islands and the associated heave of isopycnals by wind stress curl (Chidichimo et al. 2010; Kanzow et al. 2010). Variations in the MOC strength can also result from local adjustment to boundary waves propagating around ocean basins (Johnson and Marshall 2002; Elipot et al. 2013) or planetary waves propagating westward from the basin interior but with limited meridional extent (Kanzow et al. 2009; Zhao and Johns 2014b). The topic of local versus remote forcing of the MOC is linked to the issue of observing the MOC at a single latitude: is the measure of the MOC at a single latitude representative of large-scale MOC variability? Elipot et al. (2014) showed that the observed MOCs at 26° and 41°N (Willis 2010) were temporally coherent on near-annual time scales, yet the phases of their annual cycles were in quadrature, resulting in a null correlation (see also Mielke et al. 2013). In general, numerical simulation experiments clearly indicate that the latitudinal boundaries between tropical, subtropical, and subpolar gyres can break the meridional coherence of the MOC on various time scales (Bingham et al. 2007; Xu et al. 2014).

Numerical simulations are able to provide basinwide and consistent transport estimates at all latitudes (Bingham et al. 2007; Zhang 2010). In contrast, transport estimates at discrete latitudes from observational methods are not necessarily comparable. For the MOC, observational methods include 1) a net transport over a fixed depth range [measured from profiling floats at a nominal 3-month time resolution near 41°N (Willis 2010)], 2) the maximum of an overturning streamfunction [estimated from transbasin geostrophic shear, as near 26°N with the RAPID–MOC/MOCHA (Cunningham et al. 2007; Rayner et al. 2011)], 3) the transport of a physically coherent current near boundaries [such as the Deep Western Boundary Current near 39°N (Toole et al. 2011) and at 26°N (Meinen et al. 2013)], or 4) zonally integrated meridional transport across a partial basin width [as near 16°N (Send et al. 2011)]. In this study, we use some of the same observations in the North Atlantic, but we aim at estimating comparable oceanic transport quantities at each of these four latitudes (41°, 26°, 39°, and 16°N), applying the method of using ocean bottom pressure (OBP) gradients on the western boundary of the Atlantic’s basin (Hughes et al. 2013; Elipot et al. 2014). Next, we apply statistical methods to study the covariance between transport estimates, and investigate wind forcing as a driver of this covariance.

This paper is organized as follows. Section 2 presents a brief review of the concepts of overturning processes and observational principles. Section 3 presents the oceanic and atmospheric data used. Section 4 describes the methods used. Section 5 presents the results of analyses between the four transport time series by themselves. Section 6 presents the results on the statistical analyses between the four transport time series and the wind over the North Atlantic, and provides dynamical interpretation for the observed statistical linkage. Section 7 provides a summary and conclusions.

## 2. Overturning meridional transports: Concepts and observational principles

To investigate rapid coupling between wind forcing and overturning transports, it is useful to consider the velocity decomposition of Lee and Marotzke (1998) (see also Jayne and Marotzke 2001; Sime et al. 2006). Assuming that a time-dependence is implicit, the meridional velocity *υ*(*x*, *y*, *z*) is decomposed into three components:

where *H*(*x*, *y*) is the water depth at location (*x*, *y*). Each of these three terms can lead to an overturning, where overturning refers to a zonally integrated meridional transport that varies with depth. The first term represents velocities that are depth independent at each (*x*, *y*) spatial location, but its zonal integral can vary with depth because of varying topography and basin width. As an example, imagine a hypothetical ocean where the western half is 1000 m deep with a depth-independent velocity of 2 cm s^{−1} northward, and the eastern half is 2000 m deep with a depth-independent velocity of 1 cm s^{−1} southward. The resulting zonally averaged velocity profile will be 1 cm s^{−1} northward in the top 1000 m and 1 cm s^{−1} southward in the lower 2000 m, effectively forming an overturning circulation. The overturning transport from the first term in (1) is the so-called external mode, and is often associated with a barotropic gyre circulation. Conceptual examples of such circulations leading to an overturning are given by Lee and Marotzke (1998), Elipot et al. (2013), and Yang (2015).

The second velocity term in the square bracket of (1) leads to the so-called Ekman overturning. The first subterm in the bracket is the upper-ocean response to zonal wind stress, summing to a meridional Ekman flow distributed over a surface Ekman layer of unknown thickness.^{1} The second subterm in the bracket represents a local vertically uniform return flow that compensates the surface Ekman flow, thus forming an overturning circulation. As noted by Hughes et al. (2013), the Ekman return flow is a convenient mathematical representation that is not meant to be physically correct since it will be distributed over a range of depths. Killworth (2008) shows that the return flow in a simple linear frictional ocean model with flat bottom can vary strongly horizontally and vertically. In addition, the exact distribution may also depend on the time scales under consideration, as also shown by Jayne and Marotzke (2001) in an ocean general circulation model.

The final term of (1) leads to a baroclinic (i.e., vertically sheared) meridional flow, with . The velocity *υ*_{sh} consists mostly of a thermal-wind sheared velocity that is balanced by the zonal density gradient but could also include non-Ekman ageostrophic flow. In a numerical model, Lee and Marotzke (1998) find that Ekman overturning dominates the meridional overturning of the Indian Ocean. In a coupled climate model, Sime et al. (2006) find that the contributions to the MOC of each term of (1) in the Atlantic Ocean on seasonal and interannual time scales depend on the latitude under consideration.

Let us now consider the meridional geostrophic velocity *υ*_{g} from the zonal pressure gradient:

where *f* is the Coriolis parameter, and *ρ* is the water density. The zonal integral of this equation gives the geostrophic meridional mass transport per unit depth:

where *p*_{E} and *p*_{W} are the OBP on the eastern and western boundaries, respectively. Thus, *T* is given by the difference between OBP on each side *x*_{W} and *x*_{E} of an ocean basin. Overturning, by definition, is a measure not of the net flow across a given latitude, but of compensating meridional flows at different depths, meaning a zonally integrated flow that has vertical shear. Thus, to capture an overturning transport, it is not so much absolute OBP signals that are needed but rather the vertical OBP gradient alongside boundaries [see Bingham and Hughes (2008) for an extended discussion of this point]:

The sheared transport ∂*T*/∂*z* can then be formally separated into two contributions: one arising from the western boundary OBP gradient, and one from the eastern boundary OBP gradient, independently of the interior velocity field. In an ocean basin with vertical sidewalls, the vertical pressure gradient is proportional to density anomalies through the hydrostatic relation. In the presence of sloping boundaries, horizontal geostrophic velocities near the boundaries are also needed to obtain the full vertical pressure gradient (Hughes et al. 2013).

The appropriateness of using the OBP gradient method to estimate overturning was demonstrated in an ocean general circulation model (OGCM) of the North Atlantic by Bingham and Hughes (2008). They found that the western boundary OBP gradient integrated to form a layer transport representative of the MOC explained more than 90% of the interannual variability of transports calculated directly from the model velocity fields. The dominance of the western boundary OBP variance is due to more energetic flow on the western boundaries and westward accumulating variability associated with Rossby waves and eddies. From observational data, Elipot et al. (2014) found that the dominant signal of the MOC near 26° and 41°N is the geostrophic overturning, which is itself dominated by the western boundary contribution. They further demonstrated that OBP gradient time series on the western boundary, integrated within appropriate depth ranges to form transport quantities, captured a large fraction of the variability of the MOC. In particular, at 26°N, the equivalent of the western boundary OBP gradient integrated relative to and below 1000 m is representative of the variability of the MOC at semiannual, and longer, time scales.

Of the three terms in (1), the first and last terms are primarily geostrophic. In the second term, *υ*_{e} is the result of a frictional process, but the compensation term (the integral) is assumed to be geostrophic. The overturning transport estimated from vertical pressure gradients following boundaries as in (3) should therefore capture overturning transports arising from all but the *υ*_{e} contribution. In this study we investigate the covariance of western boundary pressure gradient contributions to overturning transports at four different latitudes, with respect to the wind forcing on a basin scale. Because our transport time series are only a few years long, and because of the nature of the methodologies applied, we investigate near-instantaneous velocity responses of the oceanic circulation, which we expect will be manifested in the first two terms of (1). The baroclinic ocean response to wind forcing, manifested in the third term, is mediated from the ocean interior by westward propagating planetary waves, and is delayed by months or years until it reaches the western boundary to influence the geostrophic shear estimated from the western boundary pressure gradients. For example, the North Atlantic Oscillation (NAO) atmospheric pattern drives a response in the North Atlantic Ocean characterized by anomalous horizontal circulations at the boundary between subtropical and subpolar gyres (Visbeck et al. 2003). Eventually, these velocity responses project onto the western boundary pressure, and thus influence the overturning. Instead, the mechanisms of adjustment considered here are typically deemed barotropic, as they are communicated by fast propagating barotropic waves within the ocean interior and around ocean basins boundaries (O’Rourke 2009).

## 3. Oceanic and atmospheric observations

### a. Oceanic overturning transport time series

#### 1) Derivations of transport time series at RAPID WAVE Line B, Line W, and RAPID–MOC/MOCHA

We study the basin-scale covariance of the North Atlantic MOC by considering the western boundary contribution to zonally integrated meridional transport relative to and below 1000 m, from observations at four different latitudes. The four mooring arrays from which data are used are shown in Fig. 1: Line B of the RAPID Western Atlantic Variability Experiment (WAVE) array near 42°N (Elipot et al. 2013), the Woods Hole Oceanographic Institution Line W near 39°N (Toole et al. 2011), RAPID–MOC/MOCHA near 26.5°N (Cunningham et al. 2007), and the Meridional Overturning Variability Experiment (MOVE) array at 16°N (Send et al. 2011). The common length of the transport time series from these four arrays is 1325 days (3.6 yr), so we are limited to studying processes acting on time scales less than three years and seven months (i.e., from seasonal to interannual time scales).

Elipot et al. (2013) applied (3) to derive western boundary contributions to zonally integrated meridional transport relative to and below 1000 m from Line B and Line W, two arrays separated by about 1000 km along the western boundary. The two resulting time series called *T*_{W} (39°N) and *T*_{B} (41°N) were shown to be coherent and almost in phase for all time scales from 3 months to 3.6 years. At shorter time scales, they were still coherent but with group delay estimates implying a propagation speed of 1 m s^{−1} between the two latitudes, consistent with expectations for baroclinic coastally trapped wave speeds. Elipot et al. (2014) showed subsequently that these two time series were representative of the Atlantic MOC as captured by Argo float data analyses near 41°N (Willis 2010), on semiannual time scales and longer.

A third time series of overturning transport below and relative to 1000 m, called *T*_{26}, was derived by Elipot et al. (2014) from RAPID–MOC/MOCHA, and shown to be strongly coherent and out of phase with the MOC strength, defined from the same array as the maximum of the vertically integrated streamfunction (Kanzow et al. 2010). The overturning transport *T*_{26} captured most of variance of the MOC at periods longer than 2 years. At periods of 6 months to 2 years, *T*_{26} captured most of the western boundary contribution to the geostrophic variance of the MOC.

No propagating signals were detected from the latitudes of Line B and Line W to 26°N, and while *T*_{B} and *T*_{W} were coherent with *T*_{26} on semiannual and longer time scales, there was a 90°-out-of-phase relationship resulting in a null correlation. The reasons for the coherence between Line B and Line W and 26°N was unclear.

#### 2) Derivation of the deep overturning transport time series at the MOVE array

The mooring array of the MOVE experiment located near 16°N is designed to capture the deep meridional flow in the western basin of the North Atlantic, between Guadeloupe in the Antilles to the west and the Mid-Atlantic Ridge to the east. The details of the instrumentations and moorings, as well as transport calculations and analyses can be found in Kanzow et al. (2006, 2008) and Send et al. (2011). The volume transport at the MOVE array is calculated by combining the unreferenced interior mass transport between an eastern tall density mooring (M1) located west of the mid-Atlantic ridge and a western tall density mooring just east of Guadeloupe (M3), with the volume transport estimated by direct velocity measurement (mooring M4) between mooring M3 and the continental rise between M3 and Guadeloupe. Based on water masses boundary considerations, absolute transport is derived by referencing geostrophic velocities to zero at 4950 m (Send et al. 2011).

Here we use data from moorings M3 and M4 only to derive a western boundary contribution to the overturning transport relative to and below 1000 m. First, we calculate the vertical shear of the interior transport with the east boundary density profile set to constant values where the results here are independent of the choice of constant value. Second, vertical profiles of cross-sectional velocity are calculated by linear interpolation and constant extrapolation at each time step from a discrete number of current meters on moorings M3 and M4. Those profiles are multiplied by nominal cross-sectional areas to form profiles of transport per unit depth at each mooring, which, when summed, provide a total transport profile per unit depth in the western wedge. This transport profile is differentiated in the vertical to obtain the transport shear in the wedge which is then added to the interior shear to estimate the total western boundary transport shear. This shear is then integrated from zero at a reference level of 1000 m downward to 4000 m to obtain *T*_{M}, the western boundary contribution to overturning transport relative to and below 1000 m. Note that the *T*_{M} daily time series derived here is anticorrelated (*ρ* = −0.14 with a *p* value of 0.15) with the North Atlantic Deep Water (NADW) transport time series of Send et al. (2011) for the 8 February 2002–23 June 2009 period. This may seem surprising but cross-spectral analysis (not shown) reveals that the absolute value of coherence phase between those two time series is mostly greater than 90° for time periods shorter than about 8 months (corresponding to anticorrelation at those time scales) but becomes less than 90° for longer time periods (corresponding to positive correlation). This implies that the two time series convey similar transport tendencies at longer time scales.

### b. Other data

We investigate the forcing of the overturning transports by the wind. We use the 10-m wind data from the Cross-Calibrated Multiplatform (CCMP) ocean surface wind vector product (Atlas et al. 2011), obtained from the NASA Physical Oceanography Distributed Active Archive Center (http://podaac.jpl.nasa.gov). The resolution of this product is 0.25° × 0.25° at 6-h intervals, and the region used is 0°–60°N, 0°–80°W in the North Atlantic. A 1.25° two-dimensional Gaussian smoothing window is applied at each time step and then subsampled every 0.5° to reduce the volume of the data. To match the spectral content of the transport time series, a third-order type-I Chebyshev filter with a cutoff frequency of 1 cycles per day (cpd) is applied to the time series of wind stress at each grid point, in both forward and reverse directions to ensure zero-phase distortion of signals. The wind time series are then subsampled at 12-h intervals.

We also analyze changes of the geostrophic surface circulation as revealed by absolute dynamic topography (ADT) data produced by SSALTO/Data Unification and Altimeter Combination System (DUACS) and distributed by AVISO (http://www.aviso.oceanobs.com/duacs/). Specifically, we used the merged, delayed-time, reference ADT map product at 7-day interval on a ⅓° Mercator grid. Note that we use the products before the update of April 2014. We also use the mean dynamic topography product Centre National d’Études Spatiales–Collecte Localisation Satellites (CNES-CLS09) version 1.1 (v1.1) (Rio et al. 2011).

## 4. Statistical methodologies

### a. Analytic signal and analytic correlation

We use the analytic transform (Gabor 1946) in our analyses because, as we will show in our results, this transformation conveys phase and phase difference information from temporal time series (Jacovitti and Scarano 1993; Marple 1999). It also forms the basis of the analytic eigen method described next. When *x*(*t*) is a real-valued time series, its complex-valued analytic extension *x*_{+}(*t*) is

where is the Hilbert transform of *x*(*t*):

Here, is the Cauchy principal value integral, * represents the convolution operator, and .

The analytic correlation between two zero-mean time series *x*(*t*) and *y*(*t*) is defined as the correlation between their respective analytic transforms (Jacovitti and Scarano 1993; Marple 1999):

where *E*[⋅] is the expectation or time average operator and (⋅)* is the conjugation operator. It is relatively straightforward to show that the (zero lag) analytic cross covariance is equal, up to a real factor, to the frequency integral of the cross-spectrum of *x*(*t*) and *y*(*t*). Thus, the phase of the analytic covariance, like the phase of *ρ*_{+}, is a power-weighted sum of all phases of the cross-spectrum, and will be dominated by the phases of the cross-spectrum in the frequency bands where this one has the largest power (see appendix A).

### b. Analytic extension of singular value decomposition analysis

The singular value decomposition (SVD) method is used in climate sciences to decompose the cross-covariance patterns between two real-valued scalar field variables, a left one and a right one, into statistical modes potentially revealing linear couplings between the two fields (Preisendorfer and Mobley 1988). This is also known as maximum covariance analysis (MCA; von Storch and Zwiers 2002). When the left and right fields are the same, the SVD method reduces to the empirical orthogonal function (EOF) method. A variant of the EOF method exists when the single-field variable components have undergone the analytic transform [(4)] and thus become complex-valued variables. The method is then known as complex (Barnett 1983; Horel 1984) or Hilbert (von Storch and Zwiers 2002; Hannachi et al. 2007) EOF analysis.

To the best of our knowledge, the variant of the SVD method when distinct left and right field variables have both undergone the analytic transform has not been described before, and it is named here the analytic SVD (ASVD) method. Under specific conditions, such as when signals of interest have a clear and unique periodicity, the ASVD method can be equivalent to a SVD method when one of the two fields has been lagged in time (e.g., Czaja and Frankignoul 1999) because the analytic covariance (or correlation) integrates the cross-spectrum ( appendix A). Here, the modes that will be revealed by our analyses do not have a single periodicity, and their spectra are generally red. Thus, the phase information cannot be readily interpreted as a temporal lag. Yet, the time evolution of the phase of the principal component (PC) time series of these modes still indicate a cyclic and oscillatory character of the explained variance.

The algebra necessary to conduct the ASVD analysis is standard, yet care needs to be taken because the data are complex valued (e.g., Schreier 2008). To establish our conventions, appendix B describes the ASVD method in detail. Here we note two points of importance. First, the coupling coefficient of a given mode, which measures the strength of the linear relationship between the left and right field variables for that mode, is the analytic correlation (6) between the complex-valued PC time series of the left field and the complex-valued PC time series of the right field. By construction, the coupling coefficient is real valued, and thus the PC time series are “in phase” on average. It is the patterns of the phase of the left and right singular vectors for that mode (i.e., the spatial patterns) that determine the phase lags between the individual components within each field, and between the left and right fields. The second point of importance is that we choose to decompose the wind stress (a bivariate field variable) into its rotary components (clockwise and counterclockwise) (Lilly and Olhede 2010), rather than into its Cartesian components (zonal and meridional). The reason for this choice is that applying ASVD onto Cartesian components intertwines geometric and temporal phase information of the bivariate variables, making them difficult to extricate. In contrast, ASVD applied to rotary components leads to relatively tractable elliptical modes of variance with distinguishable geometry and phase information; in particular the geometry of the variance ellipses of a given mode is the same as the geometry of the instantaneous hodographs of the vector anomalies [see Elipot and Beal (2015) for details].

### c. Spectral model and estimates

For the purpose of simulation, we fit a Matérn model to the observed transport time series *T*_{j,t} for *j* = 1, …, 4. The Matérn model (Matérn 1960) is more commonly applied to spatial data (Stein 1999) but is also reasonable for time series analysis (Sykulski et al. 2016). The spectral density of the model is

for which the parameters are usually interpreted as follows. The parameter is an overall energy level, *α*_{3} determines the smoothness or differentiability of the process, and *α*_{2} determines the range or correlation decay.

We estimate the parameter ** α** = (

*α*

_{1},

*α*

_{2},

*α*

_{3}) of the Matérn spectrum by maximizing the Whittle likelihood (Whittle 1953):

where

and is the sample mean of *T*_{j} (Sykulski et al. 2016). The sum over the indices *k* corresponds to the frequency bands achievable from the *N* data points time series. The first Slepian data taper is *h*_{0,t} (Walden 2000), used to remove leakage in the Fourier transform. A single taper for the estimation of ** α** is used because the objective of its usage is to minimize spectral leakage rather than to estimate the spectrum. The maximization of is achieved by applying the standard Nelder–Mead optimization method (Press et al. 1988). The optimum values for each transport time series are listed in Table 1.

We also estimate the auto- or cross-spectral density function of our quantities *T*_{j} by a multitaper estimate, which is formed from individual orthogonal Slepian tapers *h*_{k,t}; each individual tapered estimate is written as

A spectral estimate is formed by averaging across tapers and so we obtain (Walden 2000)

### d. Bootstrapping

Throughout this study, the Matérn spectrum model for each transport time series *T*_{j}, is used to assess the significance of the various statistics estimated from the observational data. We use a parametric approach, coupled with phase scrambling, to bootstrap whole time series (Theiler et al. 1992; Davison and Hinkley 1997, p. 408). From the Matérn model parameters obtained for each *T*_{j}, simulated replicated time series are generated as follows. The Fourier transform of a simulated time series corresponding to *T*_{j} is generated with a random phase for each discrete frequency *ν*_{k} as

where is the Matérn model for *T*_{j}, and where *Z*_{1}(*ν*) and *Z*_{2}(*ν*) are two zero-mean unit-variance Gaussian random sequences of length (*N*/2) − 1, the number of frequencies sampled, coupled with two real-valued unit variance Gaussian random sequences at *ν* = 0 and *ν* = 1/2 just multiplied by . To make the generated time series real valued, the sequence is extended to negative frequencies using Hermitian symmetry of the Fourier transform. The simulated time series is then obtained by taking the inverse Fourier transform. To avoid periodic sequences a series of twice the length of the data is generated, and half the series subsequently discarded. This operation is repeated 10^{4} times to obtain a pool of simulated time series. Typically, the statistical analyses in this study (correlation, coherence, complex empirical orthogonal function analysis, singular value decomposition) are repeated over these simulated realizations, and the distributions of the statistics from the simulations are used to assess the significance of the statistics calculated from the real observations.

## 5. Results: Relationship between transport time series

### a. Standard and analytic correlations

For analyses, we consider the original 12-hourly time series for their overlapping time period, from 22 August 2004 to 8 April 2008, and also the time series after a 3-month third-order Butterworth low-pass filter is applied forwards and backward to prevent phase distortion (Fig. 2). The 3-month cutoff corresponds to the minimum time scale at which Elipot et al. (2013) detected significant coherence between *T*_{B} and *T*_{W} but time delays not significantly different from zero. In addition, Elipot et al. (2014) found that *T*_{W} exhibited some coherence with *T*_{26} for periods longer than 2 months (although *T*_{B} exhibited significant coherence with *T*_{26} only at periods longer than 15 months). In both cases, the phase in coherent bands was found to be near −90°, implying that the overturning transports at Line B and Line W led the transport at 26°N.

Here we conduct cross-spectral analyses with the new time series *T*_{M} to find that it exhibits significant coherence at the 95% confidence level with *T*_{26} and *T*_{W} only in a few marginal frequency bands corresponding to periods longer than 2 months (not shown). As a consequence of weak coherence, the only significant correlation between the time series at 12-hourly resolution is found between *T*_{B} and *T*_{W} at 0.18 (Table 2; see also Elipot et al. 2013). The correlations of *T*_{26}, and of *T*_{M}, with the other three time series are indistinguishable from zero. The correlation between *T*_{B} and *T*_{W} increases to 0.59 for the 3-month low-pass-filtered time series, yet all other correlations remain near zero.

The realization that a specific phase organization may exist between the four time series prompts us to calculate the complex-valued analytic correlation *ρ*_{+}. The analytic correlation between all pairs of transports for the 12-hourly and 3-month low-pass-filtered time series is reported in Table 2, displaying the absolute values and complex arguments, or phases in degrees. We find that the transport time series adjacent in latitude all have modest, yet significant, analytic correlation with absolute values between 0.19 and 0.27 for the 12-h time series. Between *T*_{B} and *T*_{W} the analytic correlation phase is −15.5°, suggesting again that *T*_{W} slightly lags *T*_{B}. Between *T*_{W} and *T*_{26} the phase is −98.6° and between *T*_{26} and *T*_{M} the phase is −69.5°. The absolute values of analytic correlations for the 3-month low-pass-filtered time series are increased but overall the organization of the phases does not change much. Examining Table 2, there seems to exist an overall pattern of correlation between these time series after accounting for phase lags. Furthermore, the arrangement of these phases suggests that there could be an underlying common signal or forcing pattern at the source of these correlations.

### b. Analytic EOF analysis

To investigate whether the analytic correlations between transport pairs are representative of a common mode of variability, we apply the ASVD method (section 4) to the transport time series. Since the left and right fields for analysis are here identical, it is effectively an analytic EOF (AEOF) method, which is also known as complex or Hilbert EOF analysis (Barnett 1983; von Storch and Zwiers 2002). Because there are four transport time series, the analysis produces four modes explaining all the variance. Using our bootstrapping method to assess significance, we find that only the first mode, hereafter AEOF1, is significant at the 95% confidence level (see Table 4), explaining 36% of the variance. The eigenvector for AEOF1 is displayed on a complex plane (see Fig. 3a) scaled to represent mode anomalies in Sverdrup. AEOF1 causes typical transport anomalies between 5.6 Sv (at Line W) and 3.2 Sv (at Line B and MOVE), which are of the same order of magnitude as the standard deviations of the transport time series (5.1, 6.6, 5.6, and 7.7 Sv for *T*_{B}, *T*_{W}, *T*_{26}, and *T*_{M}, respectively). The projection of the four time series onto AEOF1 results in the first analytic PC (APC1), plotted in Figs. 3b and 3c. [The variance explained by this first mode for each time series is listed in Table 5. This mode explains the most variance for *T*_{26} (55.7%) and the least for *T*_{B} (19.1%).]

Here we choose to represent AEOF1 when the phase of the component for *T*_{B} is 180°, that is when the anomaly for *T*_{B} is southward (Fig. 3a). At those times, the phase of *T*_{W} is separated by approximately −12° from the phase of *T*_{B}, and the phase of *T*_{26} is separated by approximately −91° from the phase of *T*_{W}. Thus, *T*_{26} is approximately in quadrature phase from *T*_{B} and *T*_{W} for this mode. In addition, the phase of *T*_{M} is separated approximately by 52° from the phase of *T*_{26}, making *T*_{M} separated by about 156° from the phase of *T*_{B}. Thus, the overall picture is one of *T*_{B} and *T*_{W} in phase, and both of them in quadrature phase with *T*_{26}, and out of phase with *T*_{M}.

The time variability of this mode is given by the complex time series APC1. The phase of APC1 (Fig. 3c) follows a mixed annual to semiannual cycle, with higher-frequency variability superimposed. The amplitude of APC1 (Fig. 3b) has annual and semiannual modulations (this is more evident for the 3-month low-pass-filtered version of the PC) but also a pronounced near-monthly variability. The spectrum of APC1 is consequently red and broadband, which means that we cannot assign a single frequency to the time variability of the mode (Fig. 3d). The energy is mostly contained at low frequencies where the spectral power levels off at periods longer than 3 months. Yet, the first-moment of the spectrum—equivalent to the energy-weighted average frequency—is 1/27.7 cpd, which indicates that variability on monthly time scales is important (also indicated by a significant peak near the 34-day period).

The AEOF analysis identifies a coupling between the transport time series, not only pairwise as the analytic correlations already showed, but also between all of them, modulated in amplitude from one year to the next and also at higher frequencies, with a temporal phase that is loosely locked to an annual-to-semiannual cycle. It is tempting to interpret the phase of the eigenvector for AEOF1 as a signal propagation, as is typical in complex EOF analyses (Barnett 1983). However, this would be valid for narrowband signals only, which is not consistent with the spectrum of AEOF1 (Fig. 3d). Instead in section 6 we interpret the pattern of AEOF1 as a rapid adjustment, or response, of the meridional overturning between 16° and 41°N to basin-scale wind forcing.

### c. Fits to annual and semiannual cycles

To characterize further the seasonal variability in the transport time series, we conduct least squares fits of annual and semiannual frequency models *T*_{j}(*t*) = *A*_{j} cos(2*πνt* + *ϕ*_{j}), with *ν* = 1/365.25 and 1/182.625 cpd. The results (amplitude, phase, and amount of variance explained) are listed in Table 3 and the corresponding curves are drawn in Fig. 4. The sums of the fits for each oceanic transport time series are also included in Fig. 2. The sum of annual and semiannual cycles explain less than 20% of the variance of the 12-hourly time series, except for *T*_{M} at 16°N for which 27.6% of the variance is accounted for. When time scales shorter than 3 months are filtered out, these cycles explain between 40% and 50% of the variance of *T*_{26} and *T*_{M}, about 29% of the variance of *T*_{W}, and about 19% of the variance of *T*_{B}.

At the annual frequency, *T*_{B} and *T*_{W} are in phase, with a maximum overturning (maximum negative anomaly) at the beginning of May, and a minimum overturning at the beginning of November (Fig. 4a). For *T*_{26}, the maximum overturning occurs at the beginning of August, and for *T*_{M} in mid-October. The phase arrangement of the annual cycle is close to the phase arrangement of the AEOF1 mode described earlier (Fig. 3). These transport time series are representative of the western boundary contribution only to the overturning, yet near 26° and 41°N they exhibit the same approximate phasing as identified in the conventional MOC time series which include the variability of the eastern boundary (maximum overturning in summer, minimum in winter) (Kanzow et al. 2010; Mielke et al. 2013).

When the annual and semiannual cycles are summed, the overturning is maximum for *T*_{B} and *T*_{W} around the beginning of July and minimum in October. For *T*_{26}, the sum of the two cycles exhibits two similar minimum overturning at the beginning of April and in October, and a maximum in mid-July. For *T*_{M}, the sum of the two cycles predominantly peaks with a maximum overturning at the end of August and a minimum in May.

## 6. Results: Relationship to wind stress and wind stress curl

In this section, we investigate the relationship between the overturning transports and the wind over the North Atlantic. Figure 5 shows the mean and standard deviation fields of the filtered wind stress (Figs. 5a,b) and wind stress curl (Figs. 5c,d) for the overlapping period of the transport time series, from 22 August 2004 to 8 April 2008. The mean wind stress exhibits an anticyclonic circulation over the subtropical gyre, with westerlies north of 35°N and the trade winds to the south. Accordingly, the wind stress curl is negative over the subtropical gyre away from coastal areas, and positive over the subpolar gyre. South of 20°N the wind stress curl is mostly positive apart from over the eastern equatorial Atlantic. The variance of wind stress increases from south to north. South of 25°S the wind variance ellipses are generally oriented along the mean wind stress direction, showing the steadiness of the trade winds. In contrast, to the north of 25°N, the variance ellipses are more isotropic with no clear orientation. Like the pattern of the mean curl, the pattern of the standard deviation of the curl is not purely zonal, but exhibits a southwest–northeast tilt.

### a. Correlation patterns

Inspired by the results of the AEOF analysis, rather than considering the standard correlation, we consider the analytic correlation between transports and wind stress. Our convention is such that if the correlation is due to a narrowband oscillatory signal, a negative phase indicates that the signal propagates from *x* to *y* or equivalently that *x* precedes *y* in time.

The analytic correlations between the transport time series and both components of the wind stress ** τ** = (

*τ*

_{x},

*τ*

_{y}) and its curl

**k**⋅ ∇ ×

**are displayed in Fig. 6. The first, and striking, result is that the strongest correlation with any wind stress variable does not occur at the respective latitudes of the overturning transports. Rather, common correlation patterns appear to be associated with large spatial scales of the wind stress over the North Atlantic.**

*τ*The four overturning transport time series exhibit weak but significant analytic correlation with *τ*_{x} in near-zonal large patterns between 15° and 35°N, with phases between 0° and −90°. In addition, *T*_{B} and *T*_{W} are significantly correlated with large areas of *τ*_{x} north of 45°N, with phases between 90° and 180°. The series *T*_{W}, *T*_{26}, and *T*_{M} are significantly correlated with large regions of *τ*_{x} south of 15°N with phases between −135° and −45° for *T*_{W} and *T*_{26}, and phases between 90° and 180° for *T*_{M}. In summary, the patterns of analytic correlation with *τ*_{x} are similar for all transport time series, except that the pattern for *T*_{M} is shifted in phase.

The patterns of analytic correlation with *τ*_{y} are roughly oriented southwest to northeast (Fig. 6, middle), characteristic of the meridional structure of weather regimes (e.g., Barrier et al. 2014). Considering the region of the domain north of 20°N, for *T*_{B} and *T*_{W}, the southeast part of the domain exhibits significant analytic correlation with *τ*_{y} with a phase between −90° and 0°, and the north and northwest parts of the domain exhibit significant correlation with a phase between 0° and 90°. The series *T*_{26} and *T*_{M} also exhibit patterns of significant correlations, located in the center and in the western parts of the domain, with phases about 180° apart. The phases of this dipole for *T*_{M} are shifted by approximately −90° compared to *T*_{26}. South of 20°N, *T*_{W}, *T*_{26}, and *T*_{M} all exhibit significant correlation with *τ*_{y} but with 90° phase differences from *T*_{W} to *T*_{26} and to *T*_{M}. If one considers together the analytic correlations with both *τ*_{x} and *τ*_{y}, and shift all phases by 180°, a positive overturning anomaly (negative transport anomaly) at *T*_{B}, *T*_{W}, and *T*_{26} corresponds to an approximately in-phase large scale anticyclonic anomaly of the wind stress over the whole North Atlantic basin. This apparent pattern would be valid for *T*_{M} but with a −90° phase shift.

The correlation patterns between the transports and wind stress curl (Fig. 6, bottom) are less striking than with the wind stress components, with smaller areas with significant correlation. This may result from the added noise due to the spatial derivatives calculated for curl. Even so, there is a marked dipole pattern in the tropics for all transport time series, with centers south and north of 10°N from 90° to 180° out of phase. The phases of these dipoles are common between *T*_{B} and *T*_{W} but shifted by approximately −45° for *T*_{26} and an additional −45° for *T*_{M}, for which this dipole is broader. Another noticeable pattern of correlation for *T*_{B}, and to a lesser extent for *T*_{W}, is another dipole outside the tropics, with a center of action with phases −90° located over the eastern Atlantic at 40°N and another center with phases shifted by about 90° over the eastern Atlantic near 20°N. Interestingly, *T*_{26} is significantly correlated and in phase with a broad region of wind stress curl located above the Gulf Stream after it separates from the west coast of North America.

These geographical patterns of analytic correlation suggest a common, basinwide response of the overturning transports to the large-scale wind and wind stress forcing. This common response is further investigated next.

### b. Singular value decomposition analysis of transport covariance with the wind stress and wind stress curl

We conduct ASVD analyses between a left field constituted of collocated time series of wind stress curl and wind stress decomposed into its rotary components and a right field constituted of the four oceanic transport time series. All time series are normalized by their respective standard deviations so that the analyses are based on correlations, which equally weight all data. The total number of modes that can be considered is limited by the minimum number of individual components in one of the two coupled variable fields under study, here four for the transports. The statistical significance of each mode is assessed by repeating the ASVD calculation for the cross-correlation matrices formed between the original wind stress time series and the 10^{4} sets of simulated transport time series, and by calculating the probabilities of obtaining singular values as large as those obtained using the real transport time series (Table 4). We find no singular value as large for the first two modes with the simulated data, and thus deem these first two modes to be significant. We interpret the coupled pattern emerging from the ASVD analyses as being representative of the response of the overturning transports to wind stress forcing.

#### 1) Seasonal mode

The first mode, ASVD1 (Fig. 7), is characterized as an annual, or seasonal, mode of variability since its APC1 time series exhibit a 360° phase progression over a year (Fig. 7e). The annual cycle is less evident for the absolute values time series (Fig. 7d), although there is a tendency for APC1[∇ × ** τ**,

**] to be larger in late summer (August) of each year. The correlation between the two APC1 time series (0.51) indicates a strong coupling between the wind stress pattern and the overturning pattern for this mode.**

*τ*The wind patterns for this mode are shown in Figs. 7a and 7b for the wind stress curl and the wind stress vectors, respectively. In Fig. 7b, the geometry and typical magnitude of the wind stress pattern are indicated by variance ellipses [as drawn, they are also instantaneous hodographs; see Elipot and Beal (2015)] and the relative importance of this mode on the total wind stress variance at each pixel is given by the homogeneous correlation map (color shading). The anomalies associated with this mode are relatively strong over the equatorial region (south of 15°N), corresponding to an oscillation of the trade winds (Fig. 5). There they explain a sizable fraction of the total variance as the homogeneous correlation is generally greater than 0.5. The mode anomalies are also strong above the subtropical gyre to the west and to the northeast, as anomalous circulation cells of opposite signs, although the pattern only captures a small fraction of the total variance of the wind stress in those regions. The instantaneous wind stress anomalies are also shown in Fig. 7b (green vectors) at times when APC1[∇ × ** τ**,

**](**

*τ**t*) = 1 (i.e., with zero phase), which approximately occur in the middle of each calendar year (Fig. 7e). At such times, the wind stress anomalies consist of an anticyclonic circulation over the western subtropical gyre and a cyclonic circulation in the northeast corner of the domain, and also consist of weak forcing to the east and over the equatorial region. At later times, when the phase of APC1[∇ ×

**,**

*τ***] progresses by 90°, the instantaneous wind stress anomalies also rotate by ±90° depending on the polarity of the ellipses (cyclonic or anticyclonic). At these times, the wind stress anomalies are relatively weak in the west and north parts of the subtropical gyre, but are relatively large in the entire north equatorial region.**

*τ*The wind stress curl anomalies for this mode (Fig. 7a) consist mostly of a relatively strong zonally elongated dipole with centers at about 5° and 19°N, with phases consistent with the wind stress vector anomalies just described. The pole near 19°N has a phase near −90° while the pole near 5°N has a phase near 90°, implying a differential Ekman pumping forcing over the tropical region at one-quarter and three-quarters of the cycle of this mode. To the north, the impact of the curl for this mode is much weaker (correlation near 0.1–0.2) and exhibits a 180°-out-of-phase dipole between the center of the subtropical gyre and its northeast corner.

The overturning response to this mode is shown (Fig. 7c) with colored arrows the size of which correspond to the standard deviations of the response, and the directions of which correspond to the phases. The response is such that *T*_{B} and *T*_{W} are approximately in phase near ±180°, which implies a negative transport anomaly below 1000 m and hence a strengthening of the MOC at these latitudes. The response for *T*_{26} is offset compared to the two northern latitudes, with a phase near 135°, and the response for *T*_{M} is even farther offset with a phase near 60°. This arrangement of phase indicates primarily that the response at 16°N is instantaneously of opposite sign to the responses at the other three latitudes. The phases also indicate that within a phase cycle the response exhibits a strengthening of the overturning first occurring near 42°N, progressing south to eventually reaching 16°N, one-third of a cycle later. The magnitude of the transport response increases from north to south, from 1.5 Sv for *T*_{B} to 6.3 Sv for *T*_{M}. The amount of variance of the transport time series explained by this mode (Table 5) also increases from north to south, at 9.2% for *T*_{B} to over 50% for *T*_{26} and *T*_{M}.

We hypothesize that the results for the transports are representative of an Ekman overturning in response to large-scale patterns of wind stress forcing, varying on seasonal time scales. In an OGCM, Jayne and Marotzke (2001) showed how, at 30°N in the Pacific basin, the surface meridional Ekman transport anomalies are almost exactly compensated instantaneously by a transport calculated as the top-to-bottom vertical integral of the model velocities (after removal of near-surface Ekman velocities). To some extent, this type of barotropic adjustment was confirmed after the first year of observations of the meridional transport components at 26°N (Kanzow et al. 2007). In the model of Jayne and Marotzke (2001), the spatial structure of the seasonal variability (defined as average January conditions minus average July conditions) of the overturning streamfunction is well reproduced by a near-surface Ekman layer and a depth-independent (but still horizontally varying) meridional velocity return flow field equal to the opposite of the surface Ekman transport divided by the ocean’s depth, as in the second term of Eq. (1). The time scales associated with the Ekman overturning are very short, on the order of an inertial period for the spinup of Ekman transports, and on the order of a day (the time needed for barotropic waves to traverse a basin) for the barotropic adjustment of the depth-independent response (Jayne and Marotzke 2001; Willebrand et al. 1980). As a consequence of the spatial structures of wind forcing, and potentially of the geometry of ocean basins, Ekman overturning cells develop within basins with large-scale meridional structures that are quite distinct from the mean overturning cells (Jayne and Marotzke 2001; Sime et al. 2006). Furthermore, the vertical structure of these Ekman cells are such that an overturning transport between 1000 and 4000 m, relative to 1000 m, will not only oppose in direction the surface Ekman transport but will also exhibit substantial shear, with an amplitude depending on latitude [see as an example Fig. 4 of Jayne and Marotzke (2001) and Fig. 4 of Sime et al. (2006)]. To reiterate, even if the near-instantaneous Ekman overturning at a given latitude manifests itself as a vertically uniform return velocity at depth, the resulting deep transport on seasonal time scales may still be vertically sheared, and thus may constitute an overturning detectable by pressure gradients on basins’ boundaries.

To test the hypothesis that ASVD1 for transports corresponds to an Ekman overturning like just described, we calculate the meridional Ekman transport as a function of latitude from the instantaneous zonal wind stress anomalies shown in Fig. 7b (we use the analytic transform of the zonal wind stress for this mode, and hence the result is an analytic meridional transport which contains phase information). North of 5°N, we find that the magnitude of such Ekman transport is typically less than 0.4 Sv so it does not match in magnitude the overturning response for ASVD1. Yet, we plot the Ekman response with arbitrary constant value (Fig. 7c) and observe that the phases of the Ekman transport indicate a general pattern of northward transport between 10°N and approximately 40°N, and a southward transport between 40° and 50°N. At one-quarter cycle later for this mode, the phases of the meridional Ekman transport is rotated by 90° counterclockwise (not shown), implying little Ekman transport between 10° and 50°N but some northward transport near 10°N. We expect that a direct response of the deep overturning transports would generally be 180° out of phase with the Ekman transports. This is not exactly what we observe but it would suffice to displace southward by about 5° of latitude the overturning transports to make this picture consistent. It is possible that the northward tilt of the gyre boundary to the east, as well as the complicated bathymetry of the North Atlantic, are responsible for the mismatch between meridional Ekman flows induced by zonal stress and deep overturning transports. We still conclude that our limited observations are consistent with an Ekman type of overturning, set up on seasonal time scales.

#### 2) Modal response to the NAO

The second coupling mode, ASVD2, between wind stress and overturning transport (Fig. 8) is associated with the pattern of the NAO for the wind stress. This is demonstrated by the significant correlation (*ρ* = 0.51) between the 30-day low-pass-filtered real part (Re) of APC2[∇ × ** τ**,

**] and the NAO index [obtained from the NOAA Climate Prediction Center (http://www.cpc.ncep.noaa.gov); Fig. 9], which is slightly larger when only October–April data are used for the computation (**

*τ**ρ*= 0.58). The correspondence with the NAO is further demonstrated in Fig. 10 which shows composites of the wind stress vector and wind stress curl anomalies for positive and negative NAO daily indices on one hand, and positive and negative Re{APC2[∇ ×

**,**

*τ***]} on the other hand, for the time period of analysis. The two sets of composite anomaly maps clearly show similar patterns. In positive phases of the NAO and ASVD2, compared to the mean (Fig. 5), there are positive wind stress curl anomalies on the southern edge of the subtropical gyre and negative anomalies on the northern edge. The wind stress anomalies consist of an anomalous anticyclonic circulation centered above the northeast corner of the subtropical gyre and westward to southwestward anomalies on the southern edge of the subtropical gyre and above the north equatorial Atlantic. These wind stress anomaly patterns result in a northern shift of the mean wind stress pattern for that period.**

*τ*The oceanic overturning transport response for ASVD2 (Fig. 8c) is typically weak for *T*_{26} (0.5 Sv standard deviation), with an absolute phase less than 90° implying a northward transport anomaly, and thus a weakening of the overturning, for the wind patterns displayed in Figs. 8a and 8b. In contrast, the response is relatively strong for the other three transports and all have absolute phases larger than 90°, which implies a common strengthening of the overturning at these latitudes. At times when APC2[*T*] has phase of zero and an amplitude 1, the response of *T*_{B} is 3.9 Sv with 180° phase, the response of *T*_{W} is 3.3 Sv with 152° phase, and the response of *T*_{M} is about 3.9 Sv with −149° phase. To investigate if this mode can correspond to an Ekman overturning, the phase of the predicted meridional Ekman transport from the zonal wind stress associated with ASVD2 is calculated and displayed in the same panel. The Ekman transports consist of flows with a positive northward component between approximately 13° and 40°N (i.e., with an absolute phase less than 90°). Thus, a strengthening of overturning transports for *T*_{B}, *T*_{W}, and *T*_{M} is consistent with a compensation response to the Ekman near-surface flow, but is inconsistent with the very weak positive overturning response at 26°N.

The correlation between APC2[∇ × ** τ**,

**] and the NAO index suggests a possible alternate mechanism for the forcing of deep overturning transports. NAO positive periods, like Re{APC2[∇ ×**

*τ***,**

*τ***]} positive periods, are associated with a negative wind stress curl anomaly centered above the intergyre region between 35°N to the west and 50°N to the east, expected to spin an “intergyre gyre” anomalous anticyclonic circulation (Marshall et al. 2001). Alternatively, such anomalies can be seen as a meridional displacement of the mean circulation. We verify that this is the case for our observation period by calculating the weighted difference of ADT between positive and negative periods of Re{APC2[∇ ×**

*τ***,**

*τ***]} (Fig. 11) after removing at each grid point a fit to a sinusoidal function with annual frequency to minimize the impact of steric seasonality. The resulting ADT is depressed in the southern part of the subtropical gyre, and generally lifted in the northern part. Compared to the mean, this corresponds to a northward shift of the subtropical gyre, implying a possible spinning up on short time scales of a barotropic Sverdrup circulation (e.g., Pedlosky 1979). Ekman pumping induced by wind stress curl is balanced by meridional geostrophic flows, which have been demonstrated from GRACE observational data to project onto the OBP of the midlatitude North Atlantic (Piecuch and Ponte 2014). A barotropic circulation would have no overturning impact in a flat bottom ocean with vertical walls. Yet, in the real ocean, the difference of topography can induce vertically sheared zonally integrated transport, and thus might be responsible for our observations (typically vertically uniform flow over deep regions and western boundary return flow over shallow regions) (Elipot et al. 2013; Yang 2015). The formal dynamical link between barotropic gyre circulation and the MOC has been shown to be via the torque of the OBP, a term arising both in the vertically integrated vorticity equation and the vorticity balance of the MOC (Yeager 2015). The amount of variance explained by ASVD2 is the strongest for the two northern latitudes (see Table 5), which is consistent with expectations of NAO-type atmospheric patterns affecting regions outside of the tropics on interannual time scales.**

*τ*## 7. Summary and conclusions

The aim of this study has been to assess the meridional coherence of the MOC from an observational standpoint, and to identify the forcing of coherent variability. For this, we have derived comparable transport time series from observational arrays at four different latitudes. These transports are determined from western boundary pressure gradients, leading to the calculation of the western boundary contribution to meridional overturning transport below and relative to 1000 m (Fig. 2). At 41°, 39°, and 26°N, these time series were shown to be representative of the MOC on semiannual and longer time scales. The resulting time series overlap by only 3.6 years, limiting this study to subannual-to-interannual time scales of variability.

Over their common length, the time series spectra (Fig. 12) do not reveal any outstanding common periodicity. Yet, a simple fit to sinusoidal oscillations with period of one year (Fig. 4) suggests that the western overturning is at its maximum at the beginning of May at 41° and 39°N, is maximum at the beginning of August at 26°N, and is maximum mid-October at 16°N. While the sinusoidal fits are no proof of coherent variability, using an analytic EOF analysis we find that the four time series do covary significantly between the annual and semiannual periods (Fig. 3). This mode of variability explains a sizeable portion of the variability at individual latitudes (Table 5), and the arrangement of the phases is such that *T*_{B} and *T*_{M} are approximately in phase, *T*_{26} is in quadrature phase from *T*_{B} and *T*_{W}, and *T*_{M} is further offset to be nearly out of phase from *T*_{B} and *T*_{W}.

To investigate a possible common forcing for this mode of overturning covariance, we have considered the analytic correlation between each transport time series and winds between 10° and 60°N in the North Atlantic. We identified striking common patterns of correlation with geographic centers that are not necessarily at the same latitudes as the transport time series. The application of analytic correlation also highlights the need to properly account for phase information. Applying a newly extended method of SVD analysis, which we have here called the analytic SVD or ASVD, we identified two significant modes of covariance (Figs. 7 and 8). The first mode is a near-annual mode of oceanic overturning, which we have interpreted to be an Ekman overturning in response to a large-scale pattern of wind forcing. The second mode is related to NAO-like patterns of winds over the North Atlantic Ocean (Figs. 9–11), and we interpret the overturning response as being the result of a barotropic Sverdrup circulation, which, when it interacts with topography, projects onto the overturning transports. This second mode had a center of action at the boundary between the subtropical and subpolar gyres, forming the so-called intergyre gyre. In summary, the ASVD analysis with the wind stress and wind stress curl is able to explain more than 50% of the variance of each individual transport time series when the contributions from the first two modes are summed (Fig. 13 and Table 5). The impact of the first seasonal mode is the strongest for the two southernmost overturning time series, and the impact of the second NAO-like mode is the strongest for the two most northerly time series.

A limitation of the SVD method is that the patterns of transport and of wind stress are designed to be orthogonal, providing a constraint on the structure of second and higher modes which limits their ability to represent natural modes of variability, which may not share the same orthogonality properties. Another approach is to use the method of weather regimes, which circumvents the caveat of orthogonality by clustering data to extract recurrent and quasi-stationary patterns. Barrier et al. (2014) used this method in a forced ocean model and also found that the MOC underwent a fast wind-driven response in the form of Ekman overturning cells, spanning wide ranges of latitudes, and delineated by the latitudes of Ekman transport convergence and divergence. Despite the limitations of modal analysis, we have been able to extend the standard SVD by using phase information and applying analytic methods. Considering the relative phases was key to explain a common response of the overturning at a discrete set of latitudes.

Another limitation comes from the real nature of the observations. The hypothesized fast wind-driven barotropic response, which we believe can explain our observed modes (Eden and Willebrand 2001), is likely obscured by the baroclinic response that occurs on longer, noninstantaneous time scales, which should eventually modify the fast barotropic response (Anderson and Killworth 1977). Finally, our study ignores the eastern boundary contribution to the variability of the overturning which has been shown to be important on annual, or again seasonal, time scales (e.g., Zhao and Johns 2014a).

Despite the strength of having comparable time series representative of MOC processes, an important restriction is the limited time span of the time series used. While the patterns and ocean responses identified here are statistically significant, longer time series could improve the physical interpretation of the ocean response. In fact, *T*_{B} has been extended through the continuation of the RAPID WAVE RAPID-Scotian Line for the time period 2008–14 (Hughes et al. 2013), and observations at Line W have also continued through 2014, but the data are not yet available. Both the RAPID–MOC/MOCHA and MOVE arrays are still ongoing, and Fig. 2 shows the continuation of *T*_{26} and *T*_{M} through 2011. An interesting and noticeable feature is that both *T*_{26} and *T*_{M} show a low-frequency increase in the last half of 2009 followed by a decrease in the first half of 2010, corresponding to the exceptional decrease of the AMOC at 26.5°N (McCarthy et al. 2012; Srokosz et al. 2012), suggesting a meridional coherence of this event between 26° and 16°N. This in-phase relationship between these two latitudes does not appear to correspond to any of the two ASVD modes identified in this study, where *T*_{26} and *T*_{M} are not in phase. The exceptional downturn at 26.5°N was primarily due to a combination of anomalously negative Ekman transport, combined with an intensification of the southward return flow in the upper midocean, reflected partly into the deeper layer (McCarthy et al. 2012), and also captured by *T*_{26} (Elipot et al. 2014, their Fig. 2). Whether the same processes occurred at 16°N and can be explained by a meridional coherent response to atmospheric forcing requires further investigation.

## Acknowledgments

The authors thank Miguel Angel Morales Maqueda and Richard Williams for insightful discussions. This work was supported by the U.K. Natural Environment Research Council as part of the RAPID program. Eleanor Frajka-Williams was funded by a Leverhulme Trust Research Fellowship. Sofia Olhede was supported by an EPSRC leadership fellowship EP/I005250/1. MOVE contributions were made under Award NA10OAR4320156 from the Climate Observations Division, National Oceanic and Atmospheric Administration, U.S. Department of Commerce. Previously, MOVE was funded by the German Bundesministerium für Bildung und Forschung (Grants 03F0246A and 03F0377B). MOVE data are freely available via the international OceanSITES program.

### APPENDIX A

#### Analytic Covariance and Correlation

Consider two zero-mean continuous variables *x*(*t*) and *y*(*t*), and their two analytic transforms *x*_{+}(*t*) and *y*_{+}(*t*), respectively. Since analytic variables are complex valued, the definition of the cross-covariance function between *x*_{+}(*t*) and *y*_{+}(*t*) is the expectation of the Hermitian product, which is the product of the complex conjugate of the first variable and of the second variable (other conventions may be chosen but one of the two variables needs to be conjugated):

From the Wiener–Khinchine theorem, the expression above can be rewritten as

where is the cross-spectrum of *x*_{+}(*t*) and *y*_{+}(*t*), which can be obtained from the cross-spectrum *S*_{xy}(*f*) of *x*(*t*) and *y*(*t*) (e.g., Bendat and Piersol 1986, chapter 13):

The cross-spectrum *S*_{xy}(*f*) is a complex-valued function of frequency *f*, which can be written by convention (Jenkins and Watts 1968)

which defines the coincident or cospectrum *L*_{xy}(*f*) and the quadrature or quad-spectrum *Q*_{xy}(*f*), as well as the amplitude cross-spectrum and the phase cross-spectrum

Thus, assuming that the function is absolutely continuous for *f* ≥ 0

and, at zero lag,

The analytic cross-correlation coefficient at zero lag between *x*(*t*) and *y*(*t*) is

Since and are real valued and correspond to variances, the phase of is identical to the phase of . Thus, according to (A9), this phase is a power-weighted sum of all phases of the cross-spectrum of *x*(*t*) and *y*(*t*). Figure A1 is an illustration of this, showing a cross-spectral analysis between the transport *T*_{B} and the zonal component of wind stress at the location 31.125°N, 37.875°W where the analytic correlation between these two quantities is the largest (Fig. 6). The phase of the analytic correlation is −51.55°, which is the phase of the sum of the complex-valued cross-spectrum from the zero frequency up to approximately 0.02 cpd, the range of frequencies where the cross-spectrum has the most power.

### APPENDIX B

#### Analytic Singular Value Decomposition Analysis

We describe here the analytic SVD method (ASVD). We consider a first, or left, complex-valued field variable {*x*_{t,j}} observed at *M* locations (*j* = 1, 2, …, *M*) and *N* discrete times (*t* = 1, 2, …, *N*). This field variable is complex valued because the analytic transform (4) has been applied to each time series. If the field variable is complex valued for another reason than calculating the analytic transform, then the interpretations of the mathematical method presented here are not quite valid. Each location *j* defines an *N* × 1 data column vector,

where Δ*t* is the time interval of the time series and (⋅)^{T} is the transpose matrix operation as per usual. The *M* column vectors are subsequently combined in a *N* × *M* data matrix

We also consider a second, or right, complex-valued field variable {*y*_{t,k}}, observed at the same *N* discrete times, and at *P* locations (*k* = 1, 2, …, *P*). The *P* locations of the right field are not necessarily equal in number to, or coinciding in space with, the *M* locations of the left field. Thus, we have a second data matrix of dimensions *N* × *P*

constructed analogously to . Without further loss of generality, it is hereafter assumed that *P* ≤ *M*. Assuming that all time series have zero mean for simplicity, the *M* × *P* cross-covariance matrix between field variables {*x*_{t,j}} and {*y*_{t,k}} is

where *E*[⋅] the expectation operator, and (⋅)^{H} is the conjugate transpose matrix operation. The (*j*, *k*) component of is

where (⋅)* is the conjugate operator and the cross-covariance function at zero lag between *x*_{j}(*t*) and *y*_{k}(*t*). Note that in practice the sample cross-covariance matrix is

for which the (*j*, *k*) entry is

The truncated SVD decomposition of the cross-covariance matrix (B4) is

where is an *M* × *P* matrix, and **Λ** and are *P* × *P* matrices. If we write

then the *k*th column vector **u**_{k} is the singular vector for {*x*}, also called the left (spatial) pattern, and the *k*th column vector **v**_{k} is the singular vector for {*y*}, also called the right pattern, both for the *k*th ASVD mode. In this analytic case, and are both complex-valued matrices. The matrix is the identity matrix, and the matrices and are unitary matrices, which means that their columns are pairwise orthonormal:

The matrix **Λ** is strictly diagonal, and on its diagonal are found the *P* real-valued and positive singular values *λ*_{k} usually arranged in decreasing order. The real-valued *k*th ratio

is the squared fraction covariance (SFC) of mode *k*, usually expressed in percentage, and is interpreted as the amount of the total cross covariance that is captured by each coupling mode, characterized in space by the singular vectors.

The singular vectors provide statistical spatial patterns leading to coupled modes of covariance between the left and right fields. These patterns are modulated in time by the expansion coefficients time series or analytic principal components (APC) time series. For each mode *k*, the complex-valued *a*_{k}(*t*) and *b*_{k}(*t*) APC time series for the left and right fields, respectively, are found in the column vectors obtained by projecting the data matrices onto their respective singular vectors:

and

Those *P* vectors are combined in the *N* × *P* PC matrices:

The APC time series can be written using polar notations

where *α*_{k}(*t*) and *β*_{k}(*t*) are absolute value, or positive amplitude time series, and *χ*_{k}(*t*) and *ϕ*_{k}(*t*) are phase time series, defined by

and similarly for *β*_{k}(*t*) and *ϕ*_{k}(*t*). In some specific cases,

can define instantaneous frequencies and for mode *k*. For the ASVD method, note that *a*_{k}(*t*) and *b*_{k}(*t*) are analytic time series, which implies that their Fourier components are null for negative frequencies. Using (B13), (B8), and the unitary property (B10), direct calculations yield

Since **Λ** is diagonal, nonnegative, and real, this result implies that for a given mode *k*, the APC time series of the left and right fields are in phase on the time average. Additionally, as in standard (non-complex valued) SVD analysis, it implies that an APC time series of the left field for a given mode is uncorrelated with all the APC time series of the right field for the other modes. The strength of the coupling for mode *k* between the two fields is measured by the correlation coefficient

which is thus real valued.

In conclusion, the data matrices, that is the reconstructed variability for any mode *k* ≤ *P*, are obtained by multiplying the *k*th APC time series by the conjugates of the *k*th singular vectors:

Thus, as we do in this study, it is advantageous to represent the spatial structure of a given ASVD mode by displaying the conjugate of a singular vector and the corresponding complex-valued APC time series.

One way of presenting results from SVD analyses in general is to compute the homogeneous covariance vectors or “maps” between each field variable and its respective APC time series. Using (B12), the *M* × 1 and *P* × 1 homogeneous covariance vectors for the left and right fields for mode *k* are (Bretherton et al. 1992)

where and are the autocovariance matrices of the left and right field, respectively. The homogeneous covariance vectors become homogeneous correlation vectors when the left-hand sides of (B20) are calculated after normalizing each column of and , as well as normalizing **a**_{k} and **b**_{k}. The homogeneous correlation vectors can also be calculated from the right-hand sides of (B20) if they are respectively divided by and , and and are correlation matrices.

Alternatively, one can choose to compute the covariance vectors between each field variable and the APC time series of the other field, which are called the heterogeneous covariance vectors or maps. By using (B12), noting that , and the orthogonality property (B10) of the singular vectors, the *M* × 1 and *P* × 1 heterogeneous covariance vectors for the left and right field for mode *k* are found to be

The heterogeneous covariance vectors become heterogeneous correlation vectors when the left-hand sides of (B21) are calculated after normalizing each column of and , as well as normalizing **b**_{k} and **a**_{k}. The heterogeneous correlation vectors can also be calculated from the right-hand sides of (B21) if they are respectively divided by and and the *λ*_{k} are the singular vectors of the cross-correlation matrix . In conclusion, the left heterogeneous covariance vector is proportional to the left singular vector, and the right heterogeneous covariance vector is proportional to the right singular vector. Representing graphically the heterogeneous covariance (or correlation) vectors has the advantage of showing both the pattern of singular vectors and the strength of the linear relationship between the two fields. Note that if the coupling coefficient (B18) for a given mode is strong, the homogeneous and heterogeneous maps can appear very similar. For the case of an EOF analysis where ≡ the homogeneous and heterogeneous maps are the same.

The components of the heterogeneous covariance vectors (B21) for the left and right field variables have the same phases as the components of the singular vectors of the left and right field variables. This means that the phase patterns of the singular vectors of the left (right) field variable show the time-average phases between the left (right) field variable and the APC time series of the right (left) field variable. In contrast, (B20) show that there can exist any average phase of covariance between the left (right) field variable and the left (right) APC time series.

In section 6 we conduct an analytic SVD analysis between the transport variables (right field) and the wind stress vector (left field), which is a bivariate variable. We apply the method described above but we decompose the bivariate field variable into its time-domain rotary components (Lilly and Olhede 2010), as opposed to its Cartesian (zonal and meridional) components, to ultimately reconstruct elliptical modes of motions of the wind stress. This reconstruction is described in the appendix of Elipot and Beal (2015).

## REFERENCES

*Random Data: Analysis and Measurement Procedures.*3rd ed. John Wiley & Sons, 594 pp.

*Bootstrap Methods and Their Application.*Cambridge University Press, 582 pp.

*Spectral Analysis and Its Applications*. Holden-Day, 525 pp.

*Geophysical Fluid Dynamics*. 2nd ed. Springer, 710 pp.

*Principal Component Analysis in Meteorology and Oceanography.*Elsevier, 425 pp.

*Numerical Recipes in Fortran 77.*Cambridge University Press, 933 pp.

*Interpolation of Spatial Data: Some Theory for Kriging.*Springer Verlag, 247 pp.

*The North Atlantic Oscillation: Climatic Significance and Environmental Impact*,

*Geophys. Monogr*., Vol. 134, Amer. Geophys. Union, 113–146.

*Statistical Analysis in Climate Research*. Cambridge University Press, 484 pp.

## Footnotes

This article is licensed under a Creative Commons Attribution 4.0 license.

© 2017 American Meteorological Society.

This article is licensed under a Creative Commons Attribution 4.0 license.

^{1}

At the RAPID–MOC/MOCHA the Ekman transport calculated from wind data is evenly distributed over the upper 100 m of the ocean.