A 30-year reconstruction of the Atlantic meridional overturning circulation shows no decline

A decline in Atlantic meridional overturning circulation (AMOC) strength has been observed between 2004 and 2012 by the RAPID array with this weakened state of the AMOC persisting until 2017. Climate model and paleo-oceanographic research suggests that the AMOC may have been declining for decades or even centuries before this, however direct observations are sparse prior to 2004, giving only ‘snapshots’ of the overturning circulation. Previous studies have used linear models based on upper layer temperature anomalies to extend AMOC estimates back in time, however these ignore changes in the deep 5 circulation that are beginning to emerge in the observations of AMOC decline. Here we develop a higher fidelity empirical model of AMOC variability based on RAPID data, and associated physically with changes in thickness of the persistent upper, intermediate and deep water masses at 26°N and associated transports. We applied historical hydrographic data to the empirical model to create an AMOC time series extending from 1981 to 2016. Increasing the resolution of the observed AMOC to approximately annual shows multi-annual variability in agreement with RAPID observations, and that the downturn between 1

frequency changes, we suggest that it must represent these deeper layers. A layered model interpretation of the density structure and the associated water mass transports is shown in Figure 1b.
Here, we revisit the approach of Longworth et al. (2011) by using linear regression models to represent the AMOC, and develop the method further to include additional layers representative of the deep circulation. Section 2 describes how we trained and validated our statistical model using the RAPID dataset, and how we selected historical hydrographic data to apply 60 to the model. Section 3 describes how these hydrographic data were used to create an extended timeseries of AMOC strength from 1982 to 2016. In Sections 4 and 5, we discuss the implications of creating the longest observational timeseries of AMOC strength that incorporates variability in the deep NADW layers, and acknowledge the limitations of using an empirical model.

Model data
the sum of these four components, a depth-dependent external transport (T ext ) is added to ensure mass is conserved and that there is zero net flow across the section. The assumption of zero net flow holds on timescales longer than 10 days .
The internal geostrophic transport, T int , is estimated from dynamic height profiles. These are created by merging data from 80 individual moorings to create four profiles for each of the western and eastern boundaries and the western and eastern sides of the Mid-Atlantic Ridge. For example, at the western boundary, most data comes from instruments deployed on the WB2 mooring, but additional data from the deeper, more eastern, WBH2 and WB3 moorings are used to cover the full depth. This results in vertical profiles of temperature and salinity with sparse resolution, which are then vertically interpolated using a monthly climatology for each location. As RAPID data are vertically interpolated over a pressure grid, depth will henceforth 85 be reported in decibars (dbar) rather than metres. These four merged and interpolated temperature and salinity profiles are used to calculate dynamic height (referenced to 4820 dbar), which is then extrapolated to the surface using a seasonal climatology.
The strength of the AMOC is then the maximum of T amoc integrated over depth from the surface, i.e., the maximum of the transport streamfunction; and the AMOC depth (z amoc ) is the depth of that maximum at each time step, usually around 1100 m. The upper mid-ocean (UMO) transport (T umo ) is defined as the mid-ocean transport (T mo ) integrated between the surface 90 and the AMOC depth (Equation 2), where the mid-ocean transport is the sum of the internal, external and Western Boundary Wedge transports. This net southwards UMO transport, which includes the southward gyre recirculation and the northwards Antilles Current, is the main contributor to the AMOC (Smeed et al., 2014).
2.2 Developing the model 95 For use in the regression models, absolute salinity, conservative temperature and in situ density were calculated from the gridded in situ temperature and practical salinity data created during the RAPID calculations. As all AMOC transports are filtered during the RAPID calculation, the same Butterworth 10-day, low-pass filter was also applied to the salinity, temperature and density data; the filtered data were then averaged from a 12-hourly resolution to a monthly mean. Anomalies of these data and RAPID-estimated transports were created by subtracting the mean between 27 May 2006 and 21 February 2017, and these 100 monthly mean anomalies were used to train all our regression models. We revisited the linear regression made by Longworth et al. (2011), an ordinary least-squares (OLS) regression between the conservative temperature anomaly at 400 dbar and the thermocline transport anomaly. Longworth et al. (2011) suggested that increasing southwards thermocline transport (a negative transport anomaly) causes more warm water to recirculate close to the western boundary, and so the temperature at any particular depth within the thermocline will increase (a positive temperature 105 anomaly). They chose 400 dbar as both the mid-point of the thermocline layer, and a depth at which every profile had a sensor deployed within 50 m. We used 132 monthly mean conservative temperature anomalies at 400 dbar from the RAPID western boundary profile, compared to their 39 historic CTD profiles. Our regression however, shown in the scatter plot in Figure 2, explained only 20% of the variance of the thermocline transport anomaly, compared with 53% for the original result based on the shorter time series. To investigate whether the regression fit could be improved by using the temperature anomaly from a 110 different depth, we used an algorithm to repeat the same regression using the conservative temperature anomaly every 20 dbar from 220 to 800 dbar, and report at which depth the highest explained variance, or adjusted R 2 value, was found. The highest explained variance found by the algorithm was 51% for the regression using the western boundary temperature anomaly at 780 dbar. Changing the regression to use the density anomaly makes very little difference, increasing the adjusted R 2 to 0.54, with the anomaly again at 780 dbar. We expanded this one layer model by creating multiple linear regression models using two, three and four explanatory variables to represent two, three and four layers respectively, reflective of the water mass and circulation depth structure. We used RAPID-defined layers: the thermocline between the surface and 800 dbar; Antarctic Intermediate Water (AAIW) between 800 and 1100 dbar; Upper North Atlantic Deep Water (UNADW) between 1100 and 3000 dbar; and Lower North Atlantic Deep Water (LNADW) between 3000 and 4820 dbar (Figure 1b). UMO transport (T umo ), as the main contributor to the AMOC, was 120 preferred to thermocline transport as the dependent variable for all subsequent regression models. Variability of the AMOC is dominated by the western boundary (Elipot et al., 2014;Bryden et al., 2009);and Frajka-Williams et al. (2016) showed strong positive correlation between UMO transports and isopycnal displacements on the western boundary around 820 m, and negative correlation between LNADW transport and isopycnal displacements on the western boundary between 1500 m and the bottom. Western boundary density anomalies were thus chosen as the independent variables representing each layer, with 125 the exception of AAIW. Since the seasonal cycle of the AMOC is driven largely by eastern boundary density anomalies, which have a sub-surface maximum around 1000 dbar , the AAIW layer was represented by an eastern boundary density anomaly between 800 and 1100 dbar. The multiple linear regression equation (Equation 3) shows the four explanatory variables -three western and one eastern boundary density anomalies -representing the thermocline (ρ z1 wb ), AAIW (ρ z2 eb ), UNADW (ρ z3 wb ), and LNADW (ρ z4 wb ) layers. The superscripts z1, z2, z3, and z4 indicate the anomaly depths to be 130 identified by algorithm.
T umo (t) = α ρ z1 wb (t) + β ρ z2 eb (t) + γ ρ z3 wb (t) + ζ ρ z4 wb (t) Initially a model with two explanatory variables representing the thermocline (ρ z1 wb ) and AAIW (ρ z2 eb ) layers was implemented, then variables representing the UNADW (ρ z3 wb ), and finally the LNADW (ρ z4 wb ) were added. All possible combinations of density anomaly depths within each layer were used to run the regression to find which combination gave the maximum 135 adjusted R 2 value. For example, for the two-layer model, each western boundary density anomaly every 20 dbar between 220 and 780 dbar was combined in turn with each eastern boundary density anomaly every 20 dbar between 800 and 1080 dbar, and regressed against the UMO transport anomaly. Due to the number of iterations required by adding the fourth variable, the western boundary density anomalies representing the UNADW and LNADW in this model were chosen every 100 dbar rather than every 20.

140
The OLS regressions were checked against a number of assumptions that should hold for a linear regression model to be fit for purpose, for example, a known issue for those models based on time series is auto-correlation of residuals. We found that all our OLS models, whether using simple or multiple linear regression, showed autocorrelation of residuals, indicated by Durbin-Watson values between 0 and 1. An alternate linear regression model was used, the generalized least-squares model with auto-correlated errors (GLSAR), which models autocorrelation of residuals for a given lag (McKinney et al., 2019). For 145 each of our models, autocorrelation was significant for a lag of one month. 7 https://doi.org/10.5194/os-2020-71 Preprint. Discussion started: 14 August 2020 c Author(s) 2020. CC BY 4.0 License.

Evaluating the model
The simplest model selected by algorithm, regressing UMO transport anomaly on the western boundary density anomaly at 780 dbar gives an adjusted R 2 value of only 0.49, shown in the top time series plot ( Figure 3a). This plot also shows relatively large model prediction intervals (orange shaded area), which give the range of UMO transport anomalies that we have 95% 150 confidence will occur for that combination of boundary density anomalies. However, adding the eastern boundary density anomaly at 1020 dbar (ρ 1020 eb ) increases the maximum adjusted R 2 to 0.74 ( Figure 3b) and reduces the prediction intervals. Adding the ρ z3 wb and then the ρ z4 wb western boundary density anomalies increases the adjusted R 2 value of the model to 0.76 and 0.78 respectively (Figure 3c and d). Although increasing layers from two to three to four doesn't increase the adjusted R 2 greatly, it does reduce the standard error of the regression, from 1.85 Sv for the single layer model, to 1.32, 1.27 and 1.23 155 Sv for the two, three and four-layer models. The algorithm also selects slightly different density anomaly depths for these last two regressions: 720, 980 and 1200 dbar when three variables are used; and 740, 980, 1200 and 3000 dbar when all four are included. As using explanatory variables to represent all four layers gives a regression model that explains the greatest variance in the UMO transport anomaly and has the lowest standard error, it was selected to apply to the historical hydrographic data.
The resulting multiple linear regression equation (Equation 4) shows the depths chosen for each of the four density anomalies 160 by the algorithm, and the coefficients for each.
T umo = 40.5 ρ 740 wb − 98.6 ρ 980 eb + 46.7 ρ 1200 wb + 46.8 ρ 3000 wb , The model was tested against RAPID data not used to train it, made available following the most recent expedition and cov- values well (r=0.75). It also shows that when the eastern boundary density anomaly at 980 dbar is replaced with a monthly climatology, the trends and variability are also well captured (r=0.71). The reason that a climatology was tested will be discussed in subsection 2.5 when we describe the selection of historical hydrographic data.
Our model was trained on monthly mean density anomalies, but was to be used with hydrographic data from much shorter periods of a day or two. To evaluate how well these 'snapshot' profiles represented the longer periods, we simulated them by 170 randomly selected 20 single points from the 7961 available 12-hourly values from the most recent RAPID data. These were applied to the model and the predicted UMO compared to the observed monthly mean UMO for the same time. Bootstrapping the model prediction showed that around 65% of the observed UMO values were within the prediction interval of the corresponding model UMO. The standard deviation of the error (predicted UMO minus observed UMO) was 2.8 Sv.
The CTD profiles to be applied to the regression model to estimate UMO transport occur at irregular intervals, so to allow 175 comparisons between periods, we calculated the mean transport anomalies for a given time window. Additionally, we calculated the weighted rolling mean for the transport anomalies with the RAPID annual cycle removed, using a Gaussian distribution  comparison. Since the UMO is a transport specific to RAPID, we also estimated the AMOC by adding the model-derived To evaluate the co-variability of the density anomaly selected to represent each water mass transport and the observed UMO transport anomaly, we determined the coherence between them using a multi-taper spectrum following Percival and Walden (1998). This method reduces spectral leakage while minimising the data loss associated with other tapers. The number of tapers anomaly at 740 dbar, which represents the thermocline water mass transport, shows significant coherence at periods of around 65, 75 and 95 days; and then at all periods from 120 days to 4 years. The highest significance occurs for periods of around 65 days and 2-3 years (Figure 5a), and is in phase at all periods (Figure 5b). The eastern boundary density anomaly at 980 dbar, representing the AAIW water mass transport, shows significant coherence for periods between 100 and 120 days, 150 to 200 195 days, and between 280 days and 1.5 years, with the strongest coherence at around 160 days and just over 1 year (Figure 5c).
The coherence for this variable is mostly out of phase by 180° (Figure 5d), which is consistent with its negative coefficient. The western boundary density anomaly at 1200 dbar, representing the UNADW water mass transport, shows significant coherence for periods between 67 and 80 days, and around 90 days, and to a lesser extent around 200 days (Figure 5e), For the first two coherence peaks it is out of phase by up to -90°( Figure 5f); at the 200 day peak it is in phase. Finally, the western boundary 200 density anomaly at 3000 dbar, representing the LNADW water mass transport, shows significant coherence only for periods at just under 140 days, and between just over 1 year, or 400 days (Figure 5g), and just over 2.5 years, or 990 days. For this latter period, the co-variability is out of phase by nearly 90°, showing that increasing density at 3000 dbar is related to a weakening of the southward UMO transport (Figure 5h).

205
The importance of the LNADW in the AMOC decline compared to the UNADW  suggested an additional linear regression model between the LNADW transport anomaly and two independent variables; a western boundary density anomaly at a depth within the LNADW layer, and the Ekman transport anomaly. The Ekman transport is included in the model as Frajka-Williams et al. (2016) found that LNADW transport showed a deep baroclinic response to changes in Ekman transport. We applied a similar algorithm to the UMO regression model, repeating the regression for the deep density anomaly 210 every 20 dbar between 3000 dbar and 4820 dbar, and reporting the maximum explained variance.
For the LNADW linear regression, the optimal model with the highest adjusted R 2 value of 0.75 for a western boundary density anomaly at 3040 dbar, close to the boundary between the upper and lower North Atlantic Deep Water layers at 3000 dbar. The resulting linear regression has a standard error of 0.94 Sv (Equation 5), and the coefficients show that a positive anomaly in LNADW transport is associated with both a negative density anomaly and negative Ekman transport anomaly.

215
This means a reduction in the deep southwards LNADW flow is linked to lower density water at 3040 dbar, which we would expect with a reduction in overturning. The inverse relation between LNADW and Ekman transports reflects the statistically significant inverse correlation (r = -0.58) found by Frajka-Williams et al. (2016) between the two transports.
This model was also tested using the RAPID data between February 2017 to November 2018, and the model-predicted 220 LNADW transport anomaly shown in Figure 4b shows that it compares well with the RAPID-observed equivalent (r=0.73).
The hydrographic profiles selected for use in the UMO empirical model were also applied to this LNADW model, and the four-year and weighted rolling means calculated using the same windows.

Selecting historical hydrographic data
As Longworth et al. (2011Longworth et al. ( ) documented, between 1980 and 2017 there are many more hydrographic profiles close to the 225 western boundary than the eastern. Since the eastern boundary density anomaly shows strong seasonal variability (Pérez-Hernández et al., 2015;Chidichimo et al., 2010), replacing it in the model with a monthly climatology from the RAPID eastern boundary data, as described in the previous subsection 2.3, allows us to use all available western boundary profiles, although at the cost of losing a little of the explained UMO transport variance.
Historical hydrographic data were obtained from the World Ocean Database 2018 (WOD2018) (Boyer et al., 2018), and 230 from the datasets processed during Longworth et al. (2011), with duplicate data removed. Data were selected initially based on a region defined by latitude and longitude, from 24°N to 27°N, and from 75°W to 77.5°W. CTD profiles were then grouped by date, with a date group defined as separated by 3 days or more. From each group, we selected profiles based on similarity of

Discussion
Although the AMOC has been well-observed since 2004 by RAPID, before this, estimates of AMOC transport were restricted to approximately decadal transatlantic sections. It has been estimated that a time series of at least 60 years is necessary to detect long-term change in the AMOC due to anthropogenic global warming (Baehr et al., 2008), so extending the AMOC record into the past is crucial. Although proxies have been used to extend AMOC estimates earlier, these use one, or at most two, 315 layers to represent AMOC dynamics (Longworth et al., 2011;Frajka-Williams, 2015). However, Baehr et al. (2007) showed that deep density measurements were important in reducing the length of the time series required to detect anthropogenic change, and single layer models neglect this. In this study, we showed that an empirical regression model applied to historical hydrographic data could be used to improve the resolution of the UMO and hence AMOC transport estimates compared to the sparse transatlantic sections. In addition, by representing the deep return layers of the AMOC, the model could capture low 320 frequency changes missed by other proxy models.
To develop this empirical model of the AMOC, we regressed UMO transport on western boundary density anomalies within each of the thermocline, UNADW and LNADW layers, and an eastern boundary climatology within the AAIW layer, using an algorithm to select the best depth for each density anomaly. The selected model was then applied to historical hydrographic CTD profiles to predict UMO and hence AMOC transport strength between 1981 and 2016, at approximately annual resolution.

325
This resolution is sufficient to show decadal and multi-year variability, with a model uncertainty of around ±2.5 Sv.
There is no overall trend in either AMOC or UMO as estimated by the model, but four-year means, following Smeed et al.  -1992, 1993-2003 and 2004-2014 to -18.3, -18.7 and -18. by RAPID well, with decreasing northwards AMOC transport and decreasing deep southwards return flow balanced by an increase in southwards gyre recirculation .
Although this model increases the temporal resolution of AMOC estimates, the resolution is still coarse compared to RAPID and the time intervals between profiles are inconsistent. The longest period where no interval is greater than 1 year is October 1988 to July 1994. There are only 2 intervals longer than 2 years: September 1981 to April 1985;and February 1998 to 350 April 2001. The longest interval is 1328 days, the mean is 210 days. Although this resolution is sufficient to show multi-year variability, as shown by the four-year means, the length of some of the sampling intervals and their inconsistency means the model cannot show interannual variability reliably.

Conclusions
In conclusion, this study shows that the dynamics of the AMOC can be represented by an empirical linear regression model 355 using boundary density anomalies as proxies for water mass layer transports. More than one layer, represented by boundary density anomalies, are required to capture lower frequency changes to UMO transport. Deep density anomalies combined with Ekman transport are successful in reconstructing LNADW transport, the deepest limb of the AMOC in the subtropical North Atlantic. Previous proxies for AMOC or UMO at 26°N that rely on single layer dynamics (e.g. Frajka-Williams (2015); Longworth et al. (2011)) cannot capture this low frequency variability. This is also the case for similar reconstructions at other 360 latitudes, for example Willis (2010). Single layer dynamics are also fundamental to estimates of the AMOC that use fixed levels of no motion such as the MOVE array (Send et al., 2011) or inverted echo sounders (see McCarthy et al. (2020) for details).
We have shown the importance of the inclusion of deep density measurements in AMOC reconstructions and believe these to be key to identifying the fingerprint of anthropogenic AMOC change (e.g. Baehr et al. (2008)).
Our model, applied to historical hydrographic data, has increased the resolution of the observed AMOC between 1981 and 365 2004 from approximately decadal to approximately annual, and in doing so we have shown decadal and four-yearly variability of the AMOC and its associated layer transports. The result is the creation of an AMOC timeseries extending over 3 decades, including for the first time deep density anomalies in an AMOC reconstruction.
Our model has not revealed an AMOC decline indicative of anthropogenic climate change (Stocker et al., 2013) nor the long-term decline reported in SST-based reconstructions of the AMOC (Caesar et al., 2018). It has accurately reproduced the . However, a weakened AMOC was not the primary cause of these anomalies Holliday et al., 2020). Whether a restrengthened AMOC will ultimately have a strong impact on Atlantic climate such as was believed to have occurred in the 1990s (Robson et al., 2012) remains to be seen.