Research ArticleOCEANOGRAPHY

A stable Atlantic Meridional Overturning Circulation in a changing North Atlantic Ocean since the 1990s

See allHide authors and affiliations

Science Advances  27 Nov 2020:
Vol. 6, no. 48, eabc7836
DOI: 10.1126/sciadv.abc7836


The Atlantic Meridional Overturning Circulation (AMOC) is crucially important to global climate. Model simulations suggest that the AMOC may have been weakening over decades. However, existing array-based AMOC observations are not long enough to capture multidecadal changes. Here, we use repeated hydrographic sections in the subtropical and subpolar North Atlantic, combined with an inverse model constrained using satellite altimetry, to jointly analyze AMOC and hydrographic changes over the past three decades. We show that the AMOC state in the past decade is not distinctly different from that in the 1990s in the North Atlantic, with a remarkably stable partition of the subpolar overturning occurring prominently in the eastern basins rather than in the Labrador Sea. In contrast, profound hydrographic and oxygen changes, particularly in the subpolar North Atlantic, are observed over the same period, suggesting a much higher decoupling between the AMOC and ocean interior property fields than previously thought.


The Atlantic Meridional Overturning Circulation (AMOC) consists primarily of a northward-flowing upper limb that is warm and saline and a southward-flowing lower limb that is relatively cold and fresh. The large amount of net heat and freshwater transports associated with the AMOC is essential for regional and global climate (1). In the background of global climate change, there have been rising concerns that the AMOC is slowing or will do so in the future. Most climate models predict a weakening AMOC in the coming century (2), while the reconstructed AMOC proxies using subpolar sea surface temperature indicate that AMOC weakening may have started since the 1950s (3). At the same time, large-scale and long-term changes in water mass properties in the entire North Atlantic have been documented. At lower latitudes, waters of Antarctic origin undergo marked changes, i.e., the Antarctic Intermediate Water (AAIW) has become warmer and saltier over the past decades (4, 5), while the Antarctic Bottom Water (AABW) has become warmer and thus lighter (6). At higher latitudes, substantial deoxygenation is identified in the upper, mode, and intermediate waters over the last five decades (7), while freshening has been observed in the overflow waters since the 1960s that compose the deep lower limb of the AMOC (8).

In the early 1990s, intense cooling and freshening of the upper layer in the subpolar gyre is attributed to a strongly positive North Atlantic Oscillation (NAO). The associated intense atmospheric forcing induced deep convection and enhanced gyre circulation that cooled and freshened the upper-layer waters. The subsequent warming and salinification in the 2000s and the following record-breaking freshening since 2012 in the subpolar gyre are also in line with the convection and circulation changes induced by the NAO (9). Those hydrographic property changes are expected, in turn, to affect circulation, e.g., by cumulatively altering the density fields (10, 11). For instance, modeling studies suggest that the cooling-induced density increase in the subpolar North Atlantic (SPNA) in the 1990s may be responsible for an enhanced overturning that can be transported downstream to the subtropics (12, 13). A question emerges with the modeled AMOC changes and the ongoing marked water mass property difference: What was the AMOC state in the past decades from an observational perspective?

Observations from the subtropical RAPID array (Rapid Climate Change-Meridional Overturning Circulation and Heatflux Array–Western Boundary Time Series) suggest a weakening AMOC between 2004 and 2012 (14), but a subsequent recovery until September 2018 has left no significant declining trend of the AMOC (15). In the SPNA, the OSNAP (Overturning in the Subpolar North Atlantic Program) record is still too short for identifying long-term changes (16). In the South Atlantic, there have also been efforts to estimate the AMOC using array observations, for instance, at 34.5°S (17). Despite the fact that these array observations have revolutionized our view on the AMOC, the relatively short records, especially in the SPNA, limit our understanding of the AMOC in previous decades.

Repeated hydrographic sections provide abundant information on hydrography and tracers (e.g., oxygen and nutrients) and, at the same time, can be used to determine circulation on a time span much longer than array observations. A classical method to determine circulation from hydrography is the box inverse method (18). The application of this method requires an enclosed oceanic box bounded by hydrographic sections and ocean boundaries. Through conservation of volume, heat, and salt, the absolute geostrophic velocity can be determined on the basis of the thermal wind relation. In this study, we use repeated hydrographic sections in the subtropical and subpolar North Atlantic from the early 1990s to the mid-2010s (tables S1 to S4) to analyze AMOC and hydrographic changes over the last three decades. The subtropical section consists of repeated occupations of the World Ocean Circulation Experiments (WOCE) line A05 at about 24.5°N. The SPNA is occupied in two parts: the Labrador Sea section that consists of the WOCE line of AR7W and the recently occupied OSNAP-West line, and the eastern SPNA section that consists of the WOCE line of AR7E and the OSNAP-East line covering the Irminger and Iceland basins and Rockall Trough (Fig. 1).

Fig. 1 Hydrographic sections and array observations.

Location of WOCE hydrographic sections A05, AR7W, and AR7E and the trans-basin OSNAP and RAPID arrays. Color shading indicates the mean absolute dynamic topography (meters) between 1993 and 2018, where regions shallower than 500 m are shaded in light gray.

In the following, we shall use the repeated hydrographic sections to first examine recent hydrographic property changes in the subpolar and subtropical North Atlantic and then to perform the inverse calculation to estimate the AMOC. Putting all these results together, we shall answer the question essential to this study, i.e., do observation-based estimates show a distinctly different AMOC state in the 2000s and 2010s compared to the 1990s in the background of profound hydrographic property changes?


Hydrographic property changes

Decadal composites of salinity (see Materials and Methods) in the SPNA show that the surface layer water and the Labrador Sea Water (LSW; with a neutral density of 27.65 < γn < 27.975 kg m−3) are markedly more saline in the 2000s than in the 1990s (Fig. 2A). The regions of salinification in these water masses are accompanied by deoxygenation and warming, which can be found in the dissolved oxygen and potential temperature differences between the two decades (figs. S1A and S2A). These subpolar upper-layer hydrographic changes are in line with the different NAO phases between the 1990s and 2000s and have been discussed in previous studies (9, 19, 20). Here, we would like to emphasize two aspects that have rarely been discussed: (i) In the 2010s, signs of recent ventilation of the eastern SPNA are marked with lower salinity and potential temperature but higher oxygen concentration in the surface layer to the LSW layer (Fig. 2B and figs. S1B and S2B). However, the lower oxygen concentration and significantly higher salinity in the dense LSW layers in the Labrador and Irminger basins in the 2010s compared to the 1990s (fig. S3A) indicate that the recent convection events in the two basins are still not as strong as in the early 1990s in terms of the volume of convectively renewed water. (ii) The overflow waters (γn > 27.975 kg m−3) in the Irminger and Iceland basins are generally fresher and cooler in the 2000s than in the 1990s. The properties of the overflow waters have been linked to that of the Atlantic inflow, through the conversion from inflow to overflow in the Nordic Seas (21) and through the entrainment of the upper-layer waters in the Irminger and Iceland basins (22, 23). The decreasing salinity of the overflow waters may be a result of long-term freshening due to Greenland ice melting under global warming (8, 10). However, in the 2010s, the overflow waters in the Labrador and Irminger basins became significantly saltier compared to those in the 2000s (Fig. 2B and fig. S3C). These saltier signals may be related to the upper-layer salinification in the North Atlantic in the 2000s (9) (Fig. 2A) that had been transported to the subarctic-subpolar basins and subsequently imprint properties of the overflow waters. These changes are also evident by an increase in the oxygen concentration in the overflow waters during the 2010s, which indicates that there are more recently ventilated contributions to the overflow waters. In the subtropics, almost the entire North Atlantic Deep Water (NADW) shows freshening, while the AABW becomes warmer between the 2010s and 1990s (fig. S4). These changes resemble the characteristics noted previously (5, 6) and will not be discussed in detail here.

Fig. 2 Decadal salinity changes in the SPNA.

Salinity difference in the SPNA (A) between the 2000s and 1990s and (B) between the 2010s and 1990s. The superimposed contours are the mean neutral density surfaces calculated using all the available repeats.

AMOC calculation and validation

In addition to the water mass property changes, the hydrographic section data offer us a unique opportunity to investigate whether these changes are in association with changes in the overturning from the early 1990s to the recent period. A box inverse model is built using the repeats of the three sections occupied in nearby years (see Materials and Methods and table S4). With the subpolar sections constrained by the annual mean altimetry geostrophic velocity in each cruise year, the inverse solution is considered to represent a synoptic mean of the respective cruise year. In total, eight realizations of AMOC estimates for each of the two subpolar sections and six realizations for the subtropical section are obtained (see Materials and Methods).

The hydrography-based AMOC estimates are first validated against the direct observations by the OSNAP and RAPID array in terms of overturning streamfunction at each of the sections (see Materials and Methods). Overall, the mean overturning streamfunction of the hydrography-based estimates is in good agreement with that of the basin-wide array observations (fig. S5). Comparisons of the annual mean overturning streamfunction between the RAPID and hydrography-based estimates at 24.5°N for the four individual overlapping years reveal comparable overturning structure between the two estimates (fig. S6). It shows a clear northward upper limb in the top 1000 m of the water column, a southward returning limb between 1000 m and about 5000 m, and a northward bottom cell of AABW. The AMOC strengths of the individual years between the two estimates all agree within error bars. The mean absolute geostrophic velocity along the 24.5°N, the Labrador Sea, and the eastern SPNA sections (figs. S7 and S8) produce very robust circulation features. These circulation features are consistent with the independent observations at the respective sections, for instance, Argo float–based estimates in the Labrador Sea (24), subpolar overturning study combining hydrography and directly measured velocity data (25), and array-based western boundary current system observations in the subtropics (26). In addition, the circulation patterns are in good agreement with studies using similar hydrographic data and methods (27, 28). This justifies the use of the box inverse method to determine the AMOC in both the subtropical and subpolar regions.

Overturning in the subtropical North Atlantic

To investigate AMOC changes over the last two and half decades, we calculate the AMOC strength using the hydrography-based estimates. The AMOC strength is defined as the maximum transport in the overturning streamfunction at each section. To put the hydrography-based and array-observed results in a synthesized context, we also compute AMOC strength time series using the GECCO2 ocean state estimate (Fig. 3; see Materials and Methods) (29). Despite the fact that the RAPID observation is not assimilated (29), GECCO2 is found to reproduce the AMOC strength and variability at 26.5°N well (Fig. 3A) with a root mean square (RMS) difference of 2.8 Sv between the monthly GECCO2 and RAPID data within the overlapping period. A detailed comparison between GECCO2 and RAPID time series over different time scales (fig. S9) suggests that GECCO2 well captures the RAPID-observed AMOC variability on shorter time scales (seasonal to intra-annual, with correlation coefficient R = 0.90) and less well on interannual time scales (R = 0.70). Although GECCO2 generally shows a weakening trend from 2005 to 2012 and a succeeding recovery of the AMOC, it is still evident that GECCO2 was unable to accurately capture the observed intense weakening of the AMOC around 2009 and the subsequent rapid recovery around 2011. The mean AMOC strength with the monthly SD is 17.6 ± 3.6 Sv for RAPID and 16.8 ± 4.3 Sv for GECCO2 within the overlapping period (April 2004 to December 2014). The good agreement encourages us to use GECCO2 to extend the AMOC time series at 26.5°N back to the 1990s. The hydrography-based estimates show a mean AMOC strength of 16.3 ± 1.2 Sv (all six realizations). In spite of interannual differences, the hydrography-based AMOC strengths in the subtropics agree with the RAPID and GECCO2 results within the range of error bars in individual years (Fig. 3A). The RMS difference in the annual mean AMOC strengths between GECCO2 and hydrography is 3.4 Sv and that between RAPID and hydrography is 1.9 Sv. The consistency among the different estimates provides a mutual verification for each method in the subtropics. From January 1990 to December 2014, GECCO2 shows an insignificant declining trend of 0.10 Sv/year (based on 95% confidence level). The RAPID array observes a likely recovery of the AMOC in recent years (15), and the hydrography-based estimates show a relatively stable AMOC among the six realizations. These estimates together suggest that although we have observed a clear freshening of NADW in the subtropics between the 2010s and the 1990s, the AMOC strength in the subtropics between the two decades may not be distinctly different.

Fig. 3 AMOC strength at the subtropical and subpolar sections.

Monthly AMOC time series from trans-basin arrays (red thin lines) and GECCO2 (gray thin line) at (A) 26.5°N, (B) the eastern SPNA section, and (C) the Labrador section. The asterisks with error bars are the hydrography-based estimates in the respective years at the corresponding sections. The gray thick lines are the annual mean GECCO2 AMOC time series, with the shades indicating the SE of the annual mean values. The red thick line with shades in (A) is the annual mean RAPID AMOC time series with the SE of the annual mean. For OSNAP, the uncertainty is provided with the original data, which is determined by Monte Carlo simulations.

Overturning in the SPNA

In the SPNA, the basin-wide array observation (i.e., OSNAP, from August 2014 to April 2016) is still very short, indicating the necessity of examining the AMOC in previous years using other data sources. In contrast to the generally good agreement of subtropical AMOC estimates among GECCO2, RAPID, and hydrography, the subpolar GECCO2 AMOC estimates are substantially different from those calculated from OSNAP and hydrography. For the eastern SPNA, although GECCO2 time series does not overlap with the OSNAP record, it shows a comparable mean AMOC strength (14.2 ± 3.0 Sv; mean ± SD) with the OSNAP-East (15.6 ± 0.8 Sv) and hydrography-based estimates (15.5 ± 1.8 Sv). The RMS between the annual mean values of GECCO2 and hydrography-based AMOC is 3.2 Sv. In the Labrador Sea, in contrast, GECCO2, as in many state-of-the-art assimilation products (30), overestimates the AMOC strength (6.2 ± 4.3 Sv; mean ± SD) compared to the observations (1.6 ± 1.4 Sv for hydrography and 2.1 ± 0.3 Sv for OSNAP-West), particularly in the early 1990s (Fig. 3C). The anomalously intense overturning in the 1990s in GECCO2, similar to other assimilation models (12, 30), corresponds to anomalously strong deep convection and export of deep waters in the Labrador Sea. Note that AMOC estimates at higher latitudes among different assimilation models spread over a large range (31), and thus, one should be cautious when interpreting the assimilation results at high latitudes.

For the hydrography-based estimates, only the last two repeats of the Labrador section and the eastern SPNA section in 2014 and 2016 partly align with the OSNAP records in time. The hydrography-based AMOC strength successfully captures the array-observed increment of the AMOC at the OSNAP-East section and the decrease at the OSNAP-West section between 2014 and 2016. This encouraging agreement between the hydrography-based estimates and the OSNAP observation, together with the very robust mean circulation features (fig. S6 and table S6), could be largely attributed to the application of altimetry geostrophic velocity to constrain the inverse model (see Materials and Methods for details). As a result, it is reasonable to assume that the hydrography-based estimates are able to realistically reproduce the observed mean AMOC in the individual cruise years, and thus the changes of the AMOC in the SPNA among the repeats, for the altimetry geostrophic velocity is available throughout most of the studied period.

On interannual time scales, the hydrography-based estimates show variations of up to about 5 Sv in the eastern SPNA and about 2 Sv in the Labrador Sea. On decadal time scales, despite the application of the annual mean altimetry geostrophic velocity as constraints for the hydrography-derived AMOC estimates, the sparsely sampled hydrographic sections (eight repeats) throughout the past three decades can hardly be used to infer or exclude a trend of the AMOC in the SPNA due to the risk of aliasing effect (32). However, it is worth noting that the mean AMOC strengths in the 2010s in both the Labrador Sea (2.2 ± 1.4 Sv) and the eastern SPNA (15.8 ± 1.8 Sv) are not distinctly different from their mean AMOC strength in the 1990s (1.4 ± 1.4 Sv and 15.9 ± 1.8 Sv, respectively), according to the error bars determined as the SD of all the eight realizations. This small difference in the decadal mean AMOC strength calculated from hydrography appears to lack any clear connectivity to the profound decadal hydrographic variations in the upper-layer waters.

The hydrography-derived mean AMOC strength in the Labrador Sea (1.6 ± 1.4 Sv) and in the eastern SPNA (15.5 ± 1.8 Sv) over the studied period is similar to that of the first 2-year OSNAP-West (2.1 ± 0.3 Sv) and OSNAP-East (15.6 ± 0.8 Sv) results (16), respectively. Previous hydrography-based studies at individual sections in the SPNA show an overturning strength of 2.0 to 2.5 Sv at the Labrador section (24, 27), 16.5 ± 2.2 or 18.4 ± 3.4 Sv at about 59.5°N in the eastern SPNA (33, 34), and 11.4 to 18.5 Sv at the OVIDE (Observatory of Interannual and Decadal Variability in the North Atlantic) section that connects Greenland and Iberian Peninsula (35, 36). The hydrography-based AMOC estimates in the SPNA in this study are in line with the previous independent estimates. However, rather than focusing on only one side of the SPNA, our hydrography-based study simultaneously considers the western and eastern basins of the SPNA, allowing us to examine whether the partition of the AMOC between the western and eastern SPNA is robust over a time scale longer than the OSNAP record. The results indicate that over the past three decades, the eastern SPNA rather than the Labrador Sea has always been playing a dominant role in determining the AMOC strength in the SPNA as a whole. The relatively small portion of the Labrador Sea overturning can be attributed to a density-compensating effect in the Labrador Sea that minimizes the diapycnal transformation of water masses (37). We provide here observational evidence for the stable partition of the AMOC between the western and eastern subpolar regions on decadal time scales despite marked water mass property changes over the same time, which broadens the OSNAP results (16) to a longer-term context. In addition, despite previously noted bias in the Labrador Sea overturning, it is clear that GECCO2 is in line with observations for the eastern SPNA dominating the total subpolar overturning (fig. S10 and Fig. 3, B and C).


Despite profound upper-layer water property changes in the SPNA, our results demonstrate a comparatively stable overturning in the region during the past three decades. This circulation pattern is consistent with the weak variability in the North Atlantic Current (NAC) on decadal time scales (fig. S11) that is the primary upper-limb source of the subpolar overturning, and with a stable overflow transport in the AMOC lower limb during the last decades (23, 38). Surface buoyancy fluxes induce large variations in water mass transformation in the subpolar basins (39). However, previous studies have shown that the deep waters recirculate in the subpolar basins subsequent to their formation for up to decades before exporting to the subtropics (and thus affecting the AMOC), mainly due to the existence of the interior pathways (13, 40). This is consistent with the water mass transformation analysis (39), which shows variations of about 2 Sv in the subpolar AMOC on decadal time scales in response to the surface buoyancy forcing at high latitudes. In addition, there exists a density-compensating effect between temperature and salinity changes in the subpolar basins, which is most prominent in the Labrador Sea (27, 37) and, to a lesser extent, in the eastern subpolar gyre (41). Those characteristics of the SPNA all disfavor a marked shift of the AMOC state (42). Last, our results support the subpolar AMOC’s weak response to the Labrador Sea changes, as suggested by recent OSNAP observations (16). During the years of strong deep convection in the Labrador Sea and Irminger Sea (i.e., the early 1990s and 2014–2015) (19, 20), the AMOC in the western and eastern SPNA was not stronger compared to periods of weak convection.

Our results are in line with other independent observation-based estimates showing a stable AMOC downstream of the SPNA since the 1990s. Across the subpolar-subtropical boundary, the reconstructed AMOC time series based on satellite altimetry and hydrography along the OVIDE section hardly reveals any long-term trends from 1993 to 2017 (35, 43). In the subtropics, the altimetry-reconstructed AMOC time series at 26.5°N shows an insignificant declining trend from 1993 to 2014 (44), consistent with GECCO2. Together with the newly updated RAPID time series that shows a likely recovering AMOC (15), the altimetry-reconstructed AMOC time series supports our conclusion that the AMOC states in the subtropics were not distinctly different between the 1990s and 2010s. This, in addition, implies that altimetry sea surface height may also be used to infer AMOC changes on decadal time scales before array observations.

Overall, our results indicate that the connectivity between AMOC changes (2) and ocean interior property changes (7, 9) is weaker than previously thought. Apart from the density compensation effect of temperature and salinity discussed earlier, it is possible that hydrographic property changes especially in the SPNA are more related to heat and freshwater redistributions by the subpolar gyre circulation (45). Besides, there may be a lagged response of the overturning to property changes over longer time scales that are not resolved by those hydrographic records. Thus, further investigation will be needed to explore those possibilities. Our study highlights the necessity of continuous and combined basin-scale hydrographic surveys, array observations, as well as AMOC simulations and reconstructions to comprehensively access the connectivity between hydrographic and overturning changes in the North Atlantic across a range of latitudes and time scales.


Water mass property differences

To examine how water mass properties in the different sections have changed over the last two and half decades, we calculate decadal differences of properties (i.e., salinity, potential temperature, and dissolved oxygen). Because all the section data are on a uniform vertical grid, we first linearly interpolate all the repeats of each section (tables S1 to S3) onto a uniform horizontal grid of 0.2° longitude resolution. We then average the interpolated repeats of each section corresponding to each decade to obtain a composite of salinity, potential temperature, and dissolved oxygen for each decade. The decades are the 1990s, 2000s, and 2010s. Last, we calculate the decadal differences of each property at each section based on the decadal mean composite values (Fig. 2 and figs. S1 to S3). Note that at 24.5°N, there is only one repeat (2004) available for the 2000s. Therefore, only the differences of salinity, potential temperature, and dissolved oxygen between the 2010s and 1990s is calculated. For the Labrador section, the repeats of 2014 and 2016 follow the OSNAP-West line, which is quite different from the AR7W line. Therefore, the repeats of 2014 and 2016 are excluded from the water mass property difference calculation. For the eastern SPNA section, all the repeats follow a nearly identical cruise track in the Irminger Basin from Greenland to the Reykjanes Ridge (centered at about 31°W); east of the Reykjanes Ridge, the cruise track of the repeats in 1992, 1995, 1996, and 1997 strongly deviates from the rest of the repeats. Therefore, for the calculation of the 1990s composite, the Irminger Basin is represented as the mean of all the repeats (1991-1, 1991-2, 1992, 1995, 1996, and 1997), while the basins east of the Reykjanes Ridge are represented only by the two repeats occupied in 1991 (1991-1, 1991-2). As a result, the 1990s composite east of the Reykjanes Ridge should be viewed as a representation of the early 1990s.

Box inverse model

To determine the overturning transport along the two sections, a box inverse model is applied to an oceanic box formed by the 24.5°N section, the eastern SPNA section, the Labrador section, and ocean boundaries. All sections are occupied multiple times since 1990s. It is important to perform a box inverse model using hydrographic sections measured in nearby years, because the inverse solution for one section is sensitive to the baroclinic and horizontal shear structure of the other section (4). We therefore combine the repeats of the three sections measured in nearby years and carry out eight inverse model calculations (see table S4 for the combinations). Note that each of the 1992 and 2015 repeats of the 24.5°N section is used twice to increase the number of AMOC estimates in the SPNA.

A box inverse model is formulated in the framework of property (e.g., volume or mass, salt, and heat) conservation that combines thermal wind relation determined from hydrography and surface-layer Ekman transports. With the aid of direct velocity measurements and well-known circulation features [e.g., Florida Current (FC)], the set of conservation equations can be solved to obtain time–mean absolute transports with formal error bars (46, 47). Here, we will briefly introduce the inverse model configuration in this study with respect to conservation equations, constraints, reference level, reference velocity, and weightings. For a detailed description of the inverse model configuration with equation expressions, please refer to the previous studies (4, 47).

Conservation equations. The oceanic box is vertically divided into 17 layers in neutral density coordinate based on the characteristics of water masses (table S5). Conservation equations are formulated for volume in the whole box and in each of the 17 layers, while heat and salt anomaly conservation equations are only applied in layers 9 to 17, whose density surfaces do not outcrop over the year. Salt anomaly is used instead of salinity, because variation of salinity throughout the water column is small compared to the mean salinity. Directly using salinity values will give the salt conservation equation large weighting and result in unrealistic solutions. Salt anomaly is calculated as the difference between in situ salinity values and an area–mean salinity of individual repeats. The conservation equations take into account not only the lateral transports across the sections but also the dianeutral fluxes of volume, salt anomaly, and heat. The conservation equation of a property in a layer can be expressed asΣn=1Nhmhm+1(vnrel+vnref)CnvΔxndz+[wmCAmCmh+wm+1CAm+1Cm+1h]0(1)where hm is the depth of upper boundary of the mth layer; vnrel and vnref refer to the relative and reference geostrophic velocity at the nth station pair [a pair consists of two adjacent CTD (Conductivity, Temperature, Depth) stations], respectively; Cnv represents the vertical areal mean concentration of property C between the nth pair of CTD stations and the upper and lower boundaries of a layer; Cmhrepresents the horizontal areal mean concentration of a property C on the mth layer surface; Δx is the horizontal distance between the two stations of a station pair; and wmC is the dianeutral flux for property C at the upper boundary of the mth layer. Note that we formulate a dianeutral flux for each of the properties, namely, volume, salt anomaly, and heat, which represents a combined effect of dianeutral advection and mixing. This has been shown to be an effective parameterization for representing net dianeutral fluxes of volume, salt, and heat in inverse models (4749); Am is the horizontal area of the upper boundary of the mth layer. Equation 1 states that the volume of a layer is conserved through horizontal advection term and vertical flux term.

In Eq. 1, except vref and w, all the other quantities are known. A set of linear conservation equations for each property in the whole box and layers can be expressed in a matrix form asAx=b(2)where matrix A consists of layer areas of station pairs multiplied by property concentration and layer surface areas multiplied by property concentration, vector x consists of the unknown reference geostrophic velocity and dianeutral velocity, and vector b consists of divergence of properties due to the relative geostrophic velocity and Ekman transport.

Reference level and reference velocity. A box inverse model essentially solves the conservation equation by finding the appropriate reference geostrophic velocity and dianeutral velocity in the context of the least square method. A good choice of a reference level, therefore, sets favorable initial conditions for the inverse model to retrieve more realistic solutions (48). In the absence of a reference geostrophic velocity, an ideal choice for the reference level is a “level of no motion,” which, however, hardly exists in the ocean. A best practice may be to look for a density surface that separates water masses flowing in opposite directions. For the 24.5°N section, the neutral density surface of γn = 28.141 kg m−3 is selected as the reference level, which separates the southward-moving lower NADW and northward-moving AABW (28). For the eastern SPNA section, the surface of γn = 27.65 kg m−3 (equivalent to σθ = 27.55 kg m−3) is selected as the reference level that separates the northward-flowing NAC and the southward-moving LSW. For the Labrador section, the deepest common depth between the adjacent CTD profiles is used as the reference level (25) due to the fact that there is no apparent density surface separating water masses moving in opposite directions in this section. For all the sections, when the maximum depth of the CTD profiles is shallower than the selected reference neutral density surface, the deepest common depth of CTD pairs is used as the reference level depth. Relative geostrophic velocity is calculated at the middle point of any CTD station pair. Over steep bottom topography, the deepest common depth of a station pair may not reach the bottom depth at the middle point of the station pair. In such cases, the relative geostrophic velocity below the deepest common depth of the station pairs is set to equal the value at the common depth.

An initial reference velocity can be added to the inverse model as an a priori guess for the reference velocity. Given its least square nature, the box inverse model seeks for solutions satisfying the conservation equations but not far away from the a priori guess (48). Therefore, a properly chosen initial reference velocity facilitates solving the inverse model accurately. For the 24.5°N section, it has been shown by previous studies that zero reference velocity on the chosen isoneutral surface, combined with particular constraints (see the next section for details), is sufficient to retrieve a sensible vertical and horizontal circulation pattern and AMOC strength (4, 28, 50).

For both sections in the SPNA, the initial reference velocity is estimated using the daily level 4 altimetry surface geostrophic velocity from the Copernicus Marine Environment Monitoring Service (hereafter altimetry velocity). The altimetry velocity is first extracted according to the cruise track of each repeat. Annual mean values and the corresponding SD in the calendar year of each repeat are calculated along the cruise track. On the basis of the thermal wind relation, relative geostrophic velocity profiles are calculated relative to the reference level using the hydrographic data. Last, the initial reference velocity is calculated as the difference between the annual mean altimetry velocity and the averaged relative geostrophic velocity in the upper 100 m. The initial reference velocity serves as an a priori guess that sets the best possible initial state for the box inverse model, particularly with respect to the NAC and the boundary currents, e.g., East and West Greenland Current and Labrador Current. The inverse model adjusts these non-zero initial reference velocities in accordance with the imposed constraints and uncertainties. Note that because there is no altimetry velocity available before 1993, long-term mean (1993–2018) initial reference velocity is used for the repeats occupied before 1993.

Constraints. The inverse model formulation allows additional information on well-known circulation characteristics to be integrated as constraints that guide the inverse model to seek for physically reasonable solutions. A transport time series of the FC is estimated using submarine telephone cables across the Straits of Florida (fig. S12) (51). In this study, we constrain the FC transport using the annual mean values of the corresponding year of the repeats at the 24.5°N section. Note that the long-term (1982–2017) FC transport is very stable over interannual to decadal time scales with a mean and SD of 31.1 ± 2.4 Sv (52). Using this long-term mean value to constrain the inverse model will not result in significant changes in the inverse solution compared to using annual mean values for the individual repeats. In addition, we constrain the Deep Western Boundary Current (DWBC) transport of the 24.5°N section as −26.5 ± 13.5 Sv (“−” represents southward flow) (4, 28). This value is obtained on the basis of a variety of DWBC estimates (26, 53, 54); the uncertainty is given by half of the constraint value due to the large variations of DWBC transports among the different estimates. Only for the 2011 repeat, we applied an additional constraint to the upper NADW water transport as –12.3 ± 2.3 Sv. This value is calculated as the annual mean and SD of the original RAPID upper NADW transport time series in 2011. Without this additional constraint for the 2011 repeat, the obtained upper NADW transport would be unrealistically strong (−18.1 Sv) together with an extremely weak lower NADW transport (−0.9). After the application of this constraint, the upper NADW transport in 2011 is –14.0 Sv with an automatically corrected lower NADW transport of –6.1 Sv without further constraint. The resultant overturning streamfunction is consistent with the RAPID overturning structure in 2011 (fig. S5C). The discrepancy may be attributed to the hydrographic structure in the DWBC region observed in the 2011 repeat. The hydrographic structure combined with a zero initial reference velocity sets an initial DWBC transport that is extremely strong in the upper deep-water layer but very weak in the lower deep-water layer. This initial condition for the inverse model is too far away from the RAPID-observed overturning structure. As a result, the inverse model would be unable to obtain a solution close to the array observation without this constraint. Note that the application of this constraint results in a marginal increase in the AMOC strength by 0.7 Sv from 16.5 ± 1.5 Sv to 17.3 ± 1.4 Sv, which is not significant.

A net southward transport of −1.6 ± 0.5 Sv comes from the Arctic through the Davis Strait to the Labrador Sea along the continental shelf shallower than 500 m (27, 55). About half of this amount returns to the Arctic through the Nordic Seas in the eastern SPNA, leaving a net southward transport of about −0.8 ± 0.6 Sv through the North Atlantic, which is equivalent to the amount of the Bering Strait transport from the Pacific (16, 56). Neglecting surface freshwater fluxes due to evaporation-precipitation and river runoff, we applied the following net transport constraints for each of the sections, i.e., at the 24.5°N section, a net −0.8 ± 0.6 Sv southward transport is applied; at the eastern SPNA section, a net 0.8 ± 0.6 Sv northward transport is applied; and at the Labrador section, a net −1.6 ± 0.5 Sv southward transport is applied. Using a zero net transport at a section will result in a marginal change on the final overturning strength. For instance, setting the net −1.6 ± 0.5 Sv at the Labrador section to zero for all the repeats and keeping the net transport value at the other sections unchanged will lead to an increase in the mean Labrador overturning from 1.6 to 1.8 Sv.

For all the sections, the Ekman transport is prescribed in the first layer of the inverse model. The Ekman transport is estimated using the monthly wind stress of the National Centers for Environmental Prediction Climate Forecast System Reanalysis (NCEP/CFSR). The annual mean value of the cruise year is used. In addition, transformation of density between outcropping layers in the oceanic box due to surface heat and freshwater fluxes is also prescribed following these authors (57, 58). The air-sea heat and freshwater fluxes used for the transformation calculation are also from the monthly NCEP/CFSR data. After the transformations of density across the outcropping layer interfaces are calculated, the volume convergence in a density layer is the difference in transformation between the upper and lower boundaries of the density layer.

Weightings and error estimate. Because of the excessive number of unknowns compared to equations, a box inverse model is generally an underdetermined system of equations. The solution and error estimation depend on the imposed weightings and constraints. The weightings consist of two parts: the row weighting matrix (R) and the column weighting matrix (C). Following previous studies (47), the row weightings are determined by the reciprocal of the transport uncertainty for each layer and the imposed constraints. The transport uncertainty is set to gradually decrease from 8.2 Sv for the surface layer to 0.5 Sv for the bottom layer, and for the whole box, the uncertainty is set to 15.9 Sv (28, 50). The column weightings are determined from the uncertainty in the reference velocity and in the dianeutral velocity. For the Labrador and eastern SPNA sections, the uncertainty in the reference velocities is directly determined as the SD of the altimetry-based reference velocity divided by a square root of 2. The square root of 2 is used as a scaling factor to avoid assigning exceptionally large weightings to the regions with large variability in the surface velocity. For the 24.5°N section, the initial reference velocity is 0, and the uncertainty of the reference velocity is set to be 0.02 cm s−1, except in the boundary current regions where the uncertainty is set to be 0.04 cm s−1 to take into account the relatively larger variability (28, 50). The uncertainty for the dianeutral velocity is set to in the order of 10−5 m s−1, which gives the best conditioning of the inverse model (48, 49).

After the constraints and weightings are applied, the inverse model can be solved using the singular value decomposition method (18). The a posteriori errors of the solution can be estimated through a Gauss-Markov estimator with an error covariance matrix P asP=CCAT(ACAT+R)1AC(3)

Validation of the inverse solution against array observations

To validate the inverse calculation, we calculate the overturning streamfunction for each of the sections and compare with the direct observations by the OSNAP and RAPID array. It is particularly important to calculate the overturning streamfunction in density space for the SPNA, where northward-flowing light upper-limb waters are transformed into southward-flowing dense lower-limb waters. The strong inclination of the isopycnals in the SPNA (Fig. 2) biases the overturning in depth space. The hydrography-based overturning streamfunction for each repeat in density space is calculated as follows Ψ(γn)=γnminγnxwxev(x,γn)dxdγn(4)where γn is neutral density surfaces; the superscript min indicates an integration from the minimum neutral density; xe and xw are the position of the eastern and western boundaries of the corresponding sections, respectively; and v(x, γn) is the cross-sectional velocity as a function of horizonal position and neutral density. For the subtropics, on the other hand, the relatively flat isopycnals (fig. S3) make it marginally different when the overturning is calculated in density or depth space. To be consistent with the RAPID calculation (59), the overturning streamfunction is computed in depth space at 24.5°N. The calculation is similar to Eq. 4, except that the cross-sectional velocity in depth space v(x, z) is integrated.

Through the application of the inverse method, total mass conservation has been achieved for all the realizations (fig. S5). At 24.5°N, the hydrography-based estimates and RAPID both show a typical double-cell overturning structure with northward-moving water in the upper 1000 m, southward-moving NADW between 1000 m and about 5000 m, and again northward-moving AABW directly above the bottom. The hydrography-based overturning streamfunction is also compared with the RAPID-derived annual mean overturning streamfunction in the four individual overlapping years, namely, 2004, 2010, 2011, and 2015 (fig. S6). Note that RAPID started operation in April 2004. To resolve a full annual cycle, the 2004 mean RAPID streamfunction is calculated as the average between 1 April 2004 and 31 March 2005. In general, the hydrography-based estimates and RAPID show comparable overturning structures and strengths. Differences between the two estimation results can be seen mainly in the upper ocean (upper 1000 m). For instance, the comparatively weaker hydrography-derived overturning strength in 2004 (fig. S6A) is likely induced by the observed hydrographic structures in the upper interior ocean, because the primary northward components of the upper ocean (i.e., the FC and the meridional Ekman transport) are already constrained using the observed annual mean values. An interannual AMOC time series reconstructed from altimetry data (44) also indicates that the AMOC was in a relatively weak state around 2004 compared to a period directly afterwards. This indicates that the period from January to April 2004, which was not covered by RAPID, may be responsible for the difference in the 2004 overturning streamfunction.

At the eastern SPNA section, the mean maximum overturning transport is found at about σθ = 27.55 kg m−3. This density surface separates the northward-flowing warm and saline North Atlantic water from the net southward-flowing LSW and overflow waters (fig. S5). At the Labrador section, the overturning transport generally reaches the maximum at about the 27.68 kg m−3 isopycnal, which is consistent with the OSNAP-West array result. Below the maximum overturning, the mean southward flow in the hydrography-based estimates is relatively weak compared to that of OSNAP-West. Differences in overturning between the two estimations are, to some extent, expected, because most of the realizations (six of eight) follow the AR7W line, which deviates from the OSNAP-West array position. Note that in the Labrador Sea, the OSNAP-West streamfunction shows a southward transport of about −2.6 Sv for water lighter than 26.4 kg m−3, while the hydrography-based streamfunction shows a southward transport of −1.3 Sv. This difference is mainly due to the fact that OSNAP-West calculation includes the complete southward transport over the shallow Canadian continental shelf. The transport value is adopted as a climatological mean from the output of a model simulation (16). However, the westernmost part of the shelf waters is usually not fully covered by hydrographic measurements and is only constrained in this study as a net southward transport of −1.6 ± 0.5 Sv representing the transport from the Arctic through the Davis Strait.

The mean velocity sections at 24.5°N, the Labrador Sea, and the eastern SPNA sections also closely resemble the observed horizontal circulation features by studies in the corresponding basins (figs. S7 and S8 and table S6) (2428). Overall, the good agreement between the hydrography-based estimates and the basin-wide array observations justifies the use of the box inverse method to determine the AMOC in both subtropical and subpolar regions.

GECCO2 AMOC time series

The GECCO2 ocean state estimate is used to extend and synthesize the hydrography- and array-based AMOC estimates. To obtain the GECCO2 AMOC time series, overturning streamfunction using GECCO2 data must be first calculated. Similar to the calculation method for the hydrography-based results, the GECCO2 overturning streamfunction is also computed in density space in the SPNA and in depth space at 24.5°N. For the calculation in density space for any time step, the neutral density field along the corresponding section is first calculated. Then, the cross-sectional velocity profile at each grid point of the section is vertically binned to a uniform density grid with 0.02 kg m−3 bin width according to the density of the velocity profile. Last, the binned velocity profiles of the section are integrated in the same manner as Eq. 4. By repeating the calculation above for each time step of GECCO2, overturning streamfunction of each time step is obtained. The GECCO2 AMOC time series is then the maximum transport of the overturning streamfunction of all time steps. The calculation in depth space is directly performed using the depth coordinate of GECCO2.


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank R. J. Greatbatch and N. P. Holliday for comments on an early version of the manuscript. We also thank two anonymous reviewers for thorough examination and comments that are very helpful for improving this manuscript. Funding: This study is supported by the National Key R&D Program of China (2019YFA0606701); the National Natural Science Foundation of China (41731173); the Strategic Priority Research Program of the Chinese Academy of Sciences (XDB42000000 and XDA20060502); Key Special Project for Introduced Talents Team of Southern Marine Science and Engineering Guangdong Laboratory (Guangzhou) (GML2019ZD0306); Innovation Academy of South China Sea Ecology and Environmental Engineering, Chinese Academy of Sciences (ISEE2018PY06); the Leading Talents of Guangdong Province Program; and Independent Research Project Program of State Key Laboratory of Tropical Oceanography (LTOZZ2005). J.K. acknowledges funding by H2020 grant 727852 (Blue-Action) and 862626 (EuroSea). F.L. acknowledges funding from the Physical Oceanography Program of the U.S. NSF (grant OCE-1948335). Author contributions: Y.F. and F.L. initiated the work and collected all the data. Y.F. performed data analysis/visualization and drafted the initial manuscript. All the authors discussed the results and reviewed and edited the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: The A05 section data are provided by the CLIVAR and Carbon Hydrographic Data Office (CCHDO; at The eastern subpolar North Atlantic section data are available through CCHDO and the British Oceanographic Data Centre ( The Labrador section data are available through CCHDO and Pangaea ( GECCO2 data are available through the Integrate Climate Data Center of University Hamburg ( Argo data were collected and made freely available by the International Argo Program and the national programs that contribute to it (; The Argo Program is part of the Global Ocean Observing System. All data products from OSNAP, including those that support the main findings of this paper, are publicly available at Data from the RAPID AMOC monitoring project are funded by the Natural Environment Research Council and are freely available at The Florida Current cable and section data are made freely available on the Atlantic Oceanographic and Meteorological Laboratory web page ( and are funded by the DOC-NOAA Climate Program Office—Ocean Observing and Monitoring Division. The altimetry geostrophic velocity data are available from the European Copernicus Marine Environment Monitoring Service (
View Abstract

Stay Connected to Science Advances

Navigate This Article