120 Years of AMOC Variability Reconstructed From Observations Using the Bernoulli Inverse

changed the perspective of AMOC variability by revealing a rich spatial and temporal structure. From 2014 onwards the RAPID-MOCHA array has been complemented by the Overturning in the Abstract Determining the long-term nature of the Atlantic Meridional Overturning Circulation (AMOC) presents a major step in understanding ocean temperature forcing, and is crucial both for placing the recently observed AMOC slowdown in the context of climate change and for predicting future climate. We present a time series for AMOC strength over the last 120 years which is derived solely from observations. Application of the Bernoulli inverse to hydrography yields the general geostrophic circulation in the North Atlantic without requiring additional sea surface height information, allowing AMOC evaluation in the pre-satellite era. AMOC varies on a multidecadal timescale, leading the Atlantic Multidecadal Variability (AMV) in sea surface temperature by 2.5 years. We consider the dynamics of the implied AMOC/AMV coupling, and hypothesize that the low-frequency variability in both parameters is driven by large-scale density anomalies circulating in the AMOC.

2 of 10 Subpolar North Atlantic Program (OSNAP) array which spans the subpolar North Atlantic (SPNA). Lozier et al. (2019) and Li et al. (2021) report early OSNAP results for the AMOC which displayed significant temporal variability which, unexpectedly, was dominated by changes in the Irminger and Iceland Basins, east of Greenland. Although these programmes give an unprecedented insight into AMOC structure and change on monthly to interannual time scales, the observational record remains too short to place AMOC variability in the context of decadal or multidecadal climate change.
Longer AMOC time series estimates have been generated using proxies such as surface temperature (Caesar et al., 2018;Rahmstorf et al., 2015), sediment cores (Thornalley et al., 2018) or coastal surface height (Mc-Carthy, Haigh, et al., 2015). Caesar et al. (2021) showed that a coherent decline is evident during the 20th century in numerous AMOC proxies, suggesting the AMOC has weakened exceptionally since the onset of the industrial era. However, proxy records are not equivalent to measurements of ocean circulation. Furthermore, there is little consistency amongst proxies on multidecadal time scales (Caesar et al., 2021;Chen & Tung, 2018), so questions remain over the dynamical links between proxies and the AMOC.
Numerical models allow AMOC hindcasts, forecasts, and sensitivity tests under a variety of forcing scenarios and have provided new insights into AMOC variability over seasonal to millennial time scales (e.g., Cheng et al., 2013;Karspeck et al., 2017;Robson et al., 2018;Zhang, 2008). Of particular interest to climate modellers is the partitioning of climate indices, such as AMV, into internally versus externally (including anthropogenically) forced components (Haustein et al., 2019). While some models suggest the AMOC and AMV are coherent over multidecadal time scales (Moat et al., 2019;Zhang, 2007Zhang, , 2008, they generally underestimate the role of the AMOC in AMV  and have largely failed to capture the observed AMV pattern under realistic external forcing (Buckley & Marshall, 2015), indicating internal variability plays a significant role. Furthermore, simulations by Arzel et al. (2018) show that internal multidecadal variability may emerge spontaneously, while the recent results of Kim et al. (2019) identify ocean dynamical changes as essential in driving the AMV. In contrast, however, simulations by Booth et al. (2012), Haustein et al. (2019), Watanabe and Tatebe (2019) recreated an AMV-like response when aerosol emissions forcing was added, while Mann et al. (2021) suggest that the AMV is the result of volcanic activity perturbing the climate system. Longer AMOC time series based on ocean observations can be used to ground-truth the both simulations and proxies, and establish dynamical links between the AMOC and AMV. For example, a recent reconstruction of the AMOC in the Nordic Seas based on hydrography (Rossby et al., 2020) indicates close correlation with the AMV over the last E 70 years and found no evidence for AMOC weakening.
The North Atlantic internal density field has been measured with increasing spatio-temporal resolution throughout the 20th century, allowing computation of the historical baroclinic streamfunction. However, the scarcity of surface height observations in the pre-satellite era (before 1993) is a barrier to quantifying the barotropic component of the flow in the ocean interior. Inverse methods (e.g., Wunsch, 1978) provide a means of determining the velocity at some reference level and thus solving for the general circulation. The Bernoulli inverse (BI) method (Killworth, 1986) generates a surface height field ( E ) and hence a surface geostrophic velocity field, based on hydrography alone. Cunningham (2000) derived a BI formulation for E  as a function of conservative temperature, Θ E , and salinity, E S, and validated it in an eddy-resolving model framework.
In this study, we apply an adapted version of the BI method from Cunningham (2000) to the EN4.2.1 (hereafter EN4) analysis data set (Good et al., 2013). The resulting E   field is combined with the hydrography data to solve for the general geostrophic circulation on the zonal section at 50°N, the approximate southern boundary of the SPNA. From these fields, the AMOC is diagnosed in neutral density space. We present our methods in Section 2, establish confidence in the BI method during the modern era in Section 3, then present the AMOC timeseries, along with associated heat and freshwater fluxes, from 1900 to present, in Section 4. We describe our approach to quantifying uncertainties in Section 5, paying special attention to the uncertainty due to sparse data coverage towards the start of the record. Finally, in Section 6, we discuss the implications of our results in the broader context of multidecadal climate variability in the North Atlantic.

Data
Gridded EN4 temperature and salinity products Good et al. (2013), using Gouretski and Reseghetti (2010) corrections, were obtained from the UK Met Office website. These consist of the all observations since 1900 optimally interpolated onto a 1-degree horizontal grid with 42 depth levels, with a monthly temporal resolution. Where observations are sparse, the fields are relaxed towards a default climatology. Each datum is accompanied by an uncertainty value and a weight factor between 0 and 1, both of which reflect the error associated with gridding a varying, sometimes sparse, network of observations. According to EN4 weight values, 50°N is latitude with the best observational coverage of the Atlantic in the first half of the 20th century. Figures S3-S6 show the temporal evolution of the EN4 Θ E and E S, their uncertainties, and their weight factors at the boundaries at 50°N. EN4 Θ E and E S values on this section were corrected for any instances of unstable stratification or sub-freezing temperatures (an occasional occurrence near the western boundary).
Gridded satellite altimetry data were obtained from the Copernicus Marine Environment Monitoring Service. AMV data were accessed via the National Centre for Atmospheric Research Climate Data Guide website (Trenberth et al., 2019). The time series consists of the area-weighted North Atlantic SST between 0 and 65°N, which is detrended by subtracting the global mean SST (Trenberth & Shea, 2006).

The Bernoulli Inverse
At each time step, we identified the intersections in Θ E S  space between each neighboring pair of EN4 profiles at 50°N. The Bernoulli Function, (Θ, , ) E B S P , where E P is pressure, was evaluated on each profile (using the function "gsw_geo_strf_Cunningham" provided in the Gibbs Seawater Oceanographic Toolbox (McDougall & Barker, 2011), and differences E B  computed at each intersection point. We assume the ocean is inviscid, such that Θ E S  intersections are linked by streamlines. As E B is conserved along streamlines up to a linear dependence on sea surface height, E , we may write: for each intersection and hence solve for the implied surface height difference, E   , between the two profile locations. In general, two adjacent profiles may have many Θ E S  intersections and the implied E   may not be unique. We therefore take the mean value of E   across all intersections.

Transport and Overturning
The meridional velocity component was computed using thermal wind shear and referenced to the surface geostrophic currents implied by E   . The calculation was carried out twice, using E   from both the Bernoulli Inverse and satellite altimetry. Following Lozier et al. (2019), a background velocity was added at each time step so that the net transport through the section was zero. The mean (±1 std) value of the background velocity was +0.30 (±0.21) cm s −1 , where northward flow is positive. We consider only geostrophic flow, therefore Ekman transport is omitted from our calculations. Transports were mapped into neutral density ( E  ) classes, integrated zonally and cumulatively integrated vertically to find the overturning stream function ( E ) in density space: The maximum value of E  is defined to be the AMOC strength for each month: Heat and freshwater transports were calculated following for example, Lozier et al. (2019): , although fails to reach many of the high-frequency extrema, corresponding to abrupt periods of exceptionally high or low overturning, implying that the BI solution may not be suitable for studying intra-annual AMOC variability. However, the method largely captures the interannual and decadal variability, with E r increasing to 0.61 when comparing time series of annual mean values. Furthermore, our reconstruction is broadly coherent with the results from the RAPID-MOCHA array, which saw ), suggesting a degree of density compensation, and display significant increasing and decreasing trends respectively, indicating that the SPNA has gained heat and increasingly exported freshwater since 1900. We note that the peaks in f E Q anomaly (i.e., reduced southward freshwater flux) in the 1970s and 2010s correspond to two major SPNA freshening events (Dickson et al., 1988;Holliday et al., 2020).

Uncertainty
The clear interannual variability throughout the AMOC time series (Figure 3a) is evidence that, even in the early 20th century when hydrographic observations are relatively sparse, EN4 contains information about the geostrophic circulation beyond the climatological default. We must, however, account for potential errors introduced into our AMOC time series due the changing number and distribution of available observations.
EN4 uncertainty values already account for the error introduced when gridding sparse data (Good et al., 2013). As our time series is informed solely by the EN4 hydrography data, we seek a simple method of translating these uncertainties into an uncertainty estimate in the AMOC. Not all data points contribute equally to the AMOC signal, with variability instead dominated by the values at the boundaries ( Figures S3-S5). Moreover, the AMOC computed in depth space, which both scales and covaries ( 0.92 E r  ) with the AMOC computed in neutral density space, is a function solely of boundary Θ E and E S values once zero net transport is imposed. Our approach is therefore to propagate the uncertainty in the boundary hydrography through each stage of the calculation for the AMOC in depth space (the explicit formulae for calculating AMOC uncertainty are left to the Supporting Information S1).
Boundary Θ E , E S and E  anomalies relative to the default climatology are shown in Figures S3-S5, while corresponding weight values, which signify the influence of observation relative to climatology, are shown in Figure S6. At the eastern boundary, these fields deviate from the climatology throughout the last century, with the exception of the two world wars. At the western boundary, however, the first 10 years of the record shows very little variability beyond the climatological default, and the variability at 1000 m depth appears  Figure S6 shows that the derived envelope of AMOC uncertainty reflects the effect of data sparsity and generally decreases as observational coverage increases over time.

Discussion
We see no significant AMOC weakening trend over the last 120 years. This is in agreement with the Rossby et al. (2020) reconstruction in the Nordic Seas, and apparently contradicts the proxy reconstructions of Caesar et al. (2018), Rahmstorf et al. (2015) and Thornalley et al. (2018). However, from the 1930s onwards we see qualitative agreement with Caesar et al. (2018), with mostly a high AMOC until the 1950s, followed by a weakening throughout the 1960s and then a lesser peak around 2000. The disagreement lies between 1900 and 1930, when we found the AMOC increased by around 4 Sv (Figure 3a) whereas Caesar et al. (2018) see a high AMOC throughout. Our uncertainty estimate exceeds 3.6 Sv over this most of this period ( Figure S6). Thus, although our results do not resolve AMOC weakening over the last century, they should not be interpreted as evidence to the contrary. Further work is required to better constrain the AMOC in the early 20th century in order to clarify this issue.
Our results indicate that the AMOC at 50°N exhibited significant multidecadal variability throughout the last 120 years, and we anticipate that AMOC variability on this time scale is non-local. For example, the decrease in AMOC strength after the millennium is consistent with RAPID-MOCHA results  for the STNA (Figure 2). The increase in the 1990s is consistent with Robson et al. (2012), who attributed SPNA warming to strong AMOC, while Jackson et al. (2016) find evidence of this strengthening phase in the STNA, again suggestive of meridional coherence.
We have encountered compelling evidence that AMOC variability at 50°N leads the AMV by E 2-3 years ( Figure 3b) and is therefore likely to be closely coupled to low-frequency North Atlantic SST variability. This result is in agreement with simulations by Moat et al. (2019), who find a 95% significant correlation between the two parameters with E 5 years lag in the SPNA. These links reaffirm the hypothesis that the AMOC and AMV are coherent (Moat et al., 2019;Rossby et al., 2020;Zhang, 2007Zhang, , 2008, and lend credence to the approach of using SST as an AMOC proxy (Caesar et al., 2018;Rahmstorf et al., 2015) over multidecadal time scales. The weakening of the AMOC and decline in Θ E Q in the last E 10 years (Figures 3a-3c) is therefore an indication that the AMV will soon enter a negative phase, as argued by Bryden et al. (2020) andFrajka-Williams et al. (2017). If the AMV is forced externally, as asserted by for example, Haustein et al. (2019) and Mann et al. (2021), this would suggest that multidecadal AMOC variability is also. Aerosol-forced AMOC variability is proposed by Menary et al. (2020) based on CMIP6 results, although with a multidecadal pattern that differs to our observational analysis.
We propose that AMOC/AMV coherence over the last 120 years is due to a two-way coupling between the two parameters. It is intuitive that high overturning, and associated meridional heat transport (Figure 4a), acts to increase SST. That the corresponding upper ocean temperature anomaly feeds back into the AMOC is perhaps less obvious. The strong anticorrelation between the AMOC and m E Q suggests that a positive northward buoyancy flux anomaly through 50°N is associated with a strong AMOC signal (since buoyancy flux is proportional to m E Q but with opposite sign). Positive buoyancy anomalies in the upper SPNA will act to inhibit downwelling, impacting the AMOC signal at 50°N at a later time. Since SPNA downwelling is the primary driver of overturning (Petit et al., 2020;Zhang & Thomas, 2021), this introduces the potential for negative feedback in the AMOC. We suggest this feedback mechanism operates on multidecadal time scales, and is therefore responsible for the dominant multidecadal mode of variability. External impulses might be responsible for exciting this behavior but, as found by Mann et al. (2021), they need not by synchronous with the AMV.
Examining the conditions north of our section ( Figure 5) shows that anomalous buoyancy advection through the 50°N section is indeed closely coupled with the density field in the eastern SPNA. This agreement between circulation and the scalar density field gives some affirmation that our dynamical approach is consistent with the evolving conditions in the wider North Atlantic. The low density since the turn of the millennium is exceptional, however, in that it does not appear to be explained by the variability in  9 of 10 From a zonally integrated perspective, our findings raise the possibility that density anomalies (relative to the long-term mean) circulate along AMOC streamlines which enter then exit the SPNA on multidecadal time scales. Multidecadal AMOC variability is intrinsically tied to the evolution the density field, while the AMV is simply the associated surface temperature signature. However, as the observed AMV time series is ∼150 E years long (Trenberth et al., 2019) it is not clear whether this behavior is stable, phase-locked by timely atmospheric feedbacks (e.g., Kostov et al., 2019), or whether the period (amplitude) of variability will drift (decrease) as the along-streamline density structure is homogenised via diffusion. Furthermore, with the relationship between density and circulation appearing to change in recent years ( Figure 5), it is unclear how the AMOC/AMV coupling will respond to a continually warming climate.