1 Introduction

The last millennium is a key period for paleoclimate studies. The relative abundance of high resolution proxy data for this period with respect to earlier ones has facilitated reconstructions of relevant climate variables—most notably temperature and hydroclimate parameters—on different spatial scales (e.g., Mann et al. 2008; Cook et al. 2010; Pages 2k Consortium 2013). Estimations of past forcing factors have additionally allowed simulating the climate evolution of this period with models of different complexity (e.g., Crowley 2000; Bauer et al. 2003; González-Rouco et al. 2003b; Servonnat et al. 2010; Jungclaus et al. 2010), thus making systematic comparisons between simulations and reconstructions possible (e.g., Jansen et al. 2007; Fernández-Donado et al. 2013; Stocker et al. 2014). Climate simulations have also been used to explore the importance of external forcing and internal climate variability to explain the most relevant climate excursions of the last millennium, such as the Medieval Climate Anomaly (MCA, 950–1250) and the Little Ice Age (LIA, 1450–1850) (e.g., Miller et al. 2012; Pages 2k Consortium 2013).

Due to its large inertia and prominent role in the meridional heat transport, the Atlantic Meridional Overturning Circulation (AMOC) has been invoked as a fundamental element to understand past and future climate changes, especially in the North Atlantic. Reconstructions of past AMOC variability are thus essential to clarify the history of climate during the past millennium. The first attempts in this line aimed at reconstructing the volume transport of the Florida Current (FC), the southern expression of the Gulf Stream (Lynch-Stieglitz et al. 1999; Lund et al. 2006). This reconstruction led to the inference that the AMOC strength was about 10 % weaker during the LIA than it is today, thus supporting previous suggestions that a reduced northward oceanic heat transport might have contributed to the anomalously cold conditions over Europe at that time (Broecker 2000). A recovery of the surface AMOC during the nineteenth and twentieth centuries following a minimum in the LIA was also reconstructed from \({}^{14}\hbox {C}\) shallow-water sea-shells off North Iceland (Wanamaker et al. 2012).

Subsequently, a variety of techniques have been proposed to infer AMOC variability in an indirect manner. Some reconstructions simply rely on the use of geological proxy archives to some extent sensitive to changes in the AMOC or its components (e.g., Wanamaker et al. 2012; Dylmer et al. 2013). Others are based on model-based covariances between the AMOC and other oceanic variables that are easier to quantify, like Labrador Sea deep densities (Robson et al. 2014) or mean surface temperatures in the warming hole region (Rahmstorf et al. 2015). Yet these reconstructions generally disagree with the one in Lund et al. (2006). Rahmstorf et al. (2015), for example, showed instead an AMOC weakening during the industrial period as compared to the LIA, in agreement with results from instrumental records (Dima and Lohmann 2010) and CMIP5 historical model simulations (Drijfhout et al. 2012; Cheng et al. 2013; Jungclaus et al. 2014). These latter however generally disagree on the mechanisms underlying the AMOC variability, as well as on the AMOC sensitivity to external forcing (e.g., Swingedouw et al. 2010; Servonnat et al. 2010; Ortega et al. 2012; Hofer et al. 2011). Thus, both AMOC variability and its influence on the climate of the last millennium are still far from being clear.

Understanding the origin of the discrepancies between these reconstructions of the AMOC past variability is essential. Discrepancies might arise due to differences in the ability of the underlying methods to infer the actual variability of the oceanic circulation; for instance, some techniques might only be valid at some particular frequency range, or might largely be influenced by local climate factors. Climate simulations have recently emerged as a powerful tool to assess the validity of reconstruction techniques by providing a physically consistent pseudo-reality in which the different reconstruction steps can be reproduced, and the outcome compared to the actual target (e.g., González-Rouco et al. 2003a; Von Storch et al. 2008).

Here, we assess whether a reconstruction of the oceanic circulation can certainly be obtained by using the thermal-wind transport, and hence zonal density gradients, as predictor. This reconstruction technique constitutes the basis of earlier studies by Lund et al. (2006) and Lynch-Stieglitz et al. (2009). It directly quantifies the geostrophic transport of an oceanic current by applying the thermal wind equation to density estimates at its lateral boundaries. The advantage of this approach is that it is obtained from well-known physical relations independent of the models, in contrast to covariances based on model simulations; besides, it has already been proved valid for short, idealized model-based experiments (Hirschi and Lynch-Stieglitz 2006) and is the base for current in-situ AMOC measurements (e.g., Kanzow et al. 2007), though it is still uncertain whether this approach can be extended for estimations of variability on longer (multidecadal) timescales.

The analysis is conducted in a transient climate simulation of the last millennium with the ECHO–G atmosphere-ocean general circulation model (AOGCM) driven by changes in solar irradiance, volcanic forcing and greenhouse gases (GHGs), in which we perform pseudo-reconstructions of the AMOC and the Florida Current (FC). We thus aim to explore, first, to what extent the past spatio-temporal variability of the AMOC can be reconstructed from zonal density gradients; second, to what extent the temporal AMOC variability and its response to external forcing can be reconstructed from a proxy for the FC geostrophic volume transport; and third, whether the FC pseudo-reconstruction and the available proxy-based FC reconstruction (Lund et al. 2006) agree.

This contribution focuses on climate variability above multidecadal time scales due to the relatively coarse temporal resolution of proxies with which our results are compared (Lund et al. 2006). Previous analyses of the simulation used in this study suggested a modest influence of external forcing on the AMOC variability and highlighted the role of surface wind variations from interannual to multidecadal timescales (Ortega et al. 2012). On the latter time scales, the AMOC was found to be largely controlled by convection activity south of Greenland, itself dominated by surface-wind variations associated with a NAO-like pattern. This suggested that the AMOC variability in this simulation is ultimately controlled by atmospheric dynamics, raising the question as to whether it would be captured by procedures based on estimating the thermal-wind component.

This assessment could also have been performed using a control run accounting only for internal climate variability. However, we believe that additionally including external forcing is bound to produce larger and more realistic variations of ocean temperature and salinity, which provides a stricter test of reconstruction methods, and allows a direct comparison of model results against available reconstructions.

The work is structured as follows: Sect. 2 describes the model and the simulation used. Sect. 3 explains the theoretical basis for the decomposition of the AMOC and its application to the FC reconstruction. In Sect. 4, we analyze to what extent this decomposition is able to capture the mean AMOC as well as its variability throughout the last millennium. Furthermore, we validate the FC reconstruction technique in the model and compare it with the available reconstructions. Finally, Sect. 5 discusses the results and summarizes the main conclusions.

2 Model description and experimental design

The ECHO–G model (Legutke et al. 1999) consists of the spectral atmospheric model ECHAM4 (Roeckner et al. 1996) and the ocean model HOPE–G (Wolff et al. 1997). The atmospheric component has a T30 horizontal resolution (ca. \(3.75^{\circ }\times 3.75^{\circ }\)) and a vertical discretization that considers 19 levels extending up to 10 hPa, with five of them being located above 100 hPa. The ocean model has a horizontal resolution of about \(2.8^{\circ }\times 2.8^{\circ }\), with an enhancement of the intertropical meridional resolution toward the Equator, where it reaches a minimum grid point separation of \(0.5^{\circ }\). This refinement in resolution is intended to produce a more realistic representation of tropical and equatorial ocean currents. Its vertical discretization incorporates 20 variably spaced levels, 14 of which are located in the upper 1000 m. To avoid climate drift, the ocean component includes both heat and freshwater flux adjustments. Further details on the methodology for these corrections and their potential impacts can be found in Min et al. (2005).

We use a forced simulation with the ECHO-G AOGCM for the period 1000–1990 (González-Rouco et al. 2006). This simulation was initiated from year 1700 of a previous forced simulation with ECHO–G (González-Rouco et al. 2003a). Both runs were analyzed in depth in previous works exploring the evolution of temperature (Beltrami et al. 2006; González-Rouco et al. 2006; Fernández-Donado et al. 2013), the AMOC (Ortega et al. 2012), and the ocean heat content (Ortega et al. 2013) during the last millennium. Besides, a recent study compared the temperature response to the forcing in this and other simulations and reconstructions of the last millennium (Fernández-Donado et al. 2013). These previous investigations will allow for a better interpretation of the results described below. Further information on this run, the experimental setup, and the model can be found in González-Rouco et al. (2009).

Only a brief description of the external forcings is provided here; for further detailed discussion as well as comparison with forcing configurations used in other model experiments, the reader is referred to Fernández-Donado et al. (2013). The estimates of natural and anthropogenic forcings used to drive the model were total solar irradiance (i.e. solar constant), the radiative effects of stratospheric volcanic aerosols, and concentrations of GHGs (\(\hbox {CO}_2\), \(\hbox {CH}_4\), and \(\hbox {N}_2\hbox {O}\)). Estimations of solar forcing were taken from Bard et al. (2000) based on concentrations of the \({}^{10}\hbox {Be}\) cosmogenic isotope, spliced and scaled to Lean et al. (1995) sunspot solar constant reconstructions, as provided by Crowley (2000). Volcanic forcing is represented in terms of an effective reduction of the solar constant using estimations provided by Crowley (2000). Concentrations of GHGs are reconstructed from the analysis of air bubbles in Antarctica ice cores (Blunier et al. 1995; Etheridge et al. 1996). For the \(\hbox {N}_2\hbox {O}\) concentration, a constant value of 276.7 ppb is used before 1860, while values after 1860 are based on the historic evolution (Battle et al. 1996). Note that the final warming trend might be overestimated due to the absence of both changes in land use and anthropogenic aerosols in the simulation (Fernández-Donado et al. 2013).

The use of more recent forcing reconstructions (e.g., Krivova et al. 2007; Gao et al. 2008; Pongratz et al. 2009) or simulations and models in the Coupled Model Intercomparisons Project phase 5/Paleoclimate Modeling Intercomparison Project Phase 3 (CMIP5/PMIP3) (Taylor et al. 2012) might improve the model-data comparison but would not be expected to have an impact on the validation of the reconstruction technique performed here.

3 Methodology

This section describes the procedure used to reconstruct the variability of the Atlantic Ocean circulation. We first summarize the theoretical basis for the reconstruction of the AMOC and then focus on the technique used to infer past variations of the FC.

3.1 AMOC decomposition

In ocean general circulation models, the AMOC streamfunction is computed from the simulated three dimensional meridional velocity field. In the real ocean, however, the lack of observational data hampers the determination of velocity profiles. In order to infer recent and past changes in the AMOC, Hirschi and Marotzke (2007) proposed two alternative ways to calculate the AMOC streamfunction. Both methods are based on a decomposition of the AMOC streamfunction, \(\psi (z)\), into three dynamical components which allows evaluating the relative contributions of different physical processes (e.g., Lee and Marotzke 1998; Baehr et al. 2004; Sime et al. 2006). These three terms correspond to the meridional components of the barotropic (depth averaged) velocity \(\overline{ v}\), also known as external mode, the geostrophic shear \(v_{sh}\), and the Ekman transport \(v_{ek}\):

$$\begin{aligned} \psi (z)&\,= \, \psi _{ex} (z) \, + \, \psi _{sh} (z) \, + \, \psi _{ek} (z) \, \nonumber \\&\,= \, \int ^{z} _{-H} dz^{\prime} \int ^{x_e} _{x _w} \overline{v} \, dx \, + \int ^{z} _{-H} dz^{\prime} \int ^{x _e} _{x _w} v _{sh} \, dx \,\nonumber \\&\quad + \int ^{z} _{-H} dz^{\prime} \int ^{x _e} _{x _w} (v _{ek} - \overline{v} _{ek} ) \, dx \, , \end{aligned}$$
(1)

where \(x_e\), \(x_w\), and H are respectively the western and eastern limits of the basin and the oceanic bottom depth. The external mode accounts for the impact of zonally nonuniform topography and frictional effects at the ocean boundaries on the AMOC circulation (for further details see Lee and Marotzke 1998). The Ekman term represents the meridional transport forced by the zonal wind stress; it is assumed to take place in the upper ocean and compensated by a barotropic flow. Finally, the geostrophic shear largely consists in thermal wind shear balanced by zonal density gradients, except near the equator.

By assuming thermal-wind balance (Marotzke et al. 1999), the zonally integrated geostrophic shear can be estimated just from seawater densities at the western and the eastern boundaries of the ocean (\(\rho _w\) and \(\rho _e\), respectively), where they can be measured and reconstructed (e.g., Lynch-Stieglitz 2001; Cunningham et al. 2007). One approximation to the zonally averaged velocity component is given by the formula:

$$\tilde{v} (x,z) = - \frac{g}{\rho ^* f} \int ^{z} _{-H} \frac{1}{L(z^{\prime})} \left( \rho _e - \rho _w \right) \, dz^{\prime},$$
(2)

where L(z) is the basin width, g the gravity acceleration at the surface of the Earth, f the Coriolis parameter, and \(\rho ^*\) the reference seawater density. In general, \({\tilde{v}}\) contains both baroclinic and barotropic contributions and is a good approximation to the velocity field if bottom velocities are close to zero (Baehr et al. 2004). Two alternative corrections were proposed by Hirschi and Marotzke (2007) to \({\tilde{{v}}}\). The first one makes \(\tilde{v}\) equivalent to the shear velocity \(v_{sh}\) in Eq. (1) by subtracting the average of \(\tilde{v}\) over all depths:

$$\overline{v}^{\prime} (x) = \frac{1}{H (x)} \int ^{0} _{-H} \tilde{v} (x,z^{\prime}) \, dz^{\prime} .$$
(3)

In this way, however, valuable information on the barotropic mode is lost. The alternative is to calculate a streamfunction directly from Eq. (2), thereby retaining the barotropic mode. To this end, and to ensure mass conservation that is necessary to define the streamfunction, the transport imbalance through a latitude-depth section, further divided by the corresponding section area, must be subtracted:

$$\hat{v} = \tilde{\psi } ( 0 ) / A ,$$
(4)

where \(\tilde{\psi }\)(0) is obtained from the following streamfunction:

$$\tilde{\psi } (0) = \int ^{x _e} _{x _w} dx \int ^{0} _{-H} \tilde{v} \, dz^{\prime} .$$
(5)

By considering these two alternative corrections for \(\tilde{v}\) defined by Eqs. (3) and (4), two streamfunctions can be calculated for the thermal-wind transport, respectively:

$$\psi _{tw _1} (z)\,= \, \int ^{x _e} _{x _w} dx \int ^z _{-H} (\tilde{v} - \hat{v} ) dz^{\prime},$$
(6)
$$\psi _{tw _2} (z)\,= \, \int ^{x _e} _{x _w} dx \int ^z _{-H} (\tilde{v} - \overline{v}^{,} ) dz^{\prime}.$$
(7)

The meridional overturning transport can thus be decomposed in either of the following ways:

$$\psi _{1}\,=\, \psi _{tw _1} + \psi _{ek},$$
(8)
$$\psi _{2}\,= \, \psi _{ex} + \psi _{tw _2} + \psi _{ek}.$$
(9)

While \(\psi _{tw_1}\) still contains a barotropic contribution that is representative of the sum of the first and second terms on the right hand side of Eq. (1), \(\psi _{tw_2}\) corresponds directly to the second term on the right-hand side of Eq. (1). As described by Hirschi and Marotzke (2007), \(\psi _{2}\) is expected to yield a better representation of the AMOC field. However, it requires knowing the barotropic velocity field, which is not easily measurable. By contrast, \(\psi _{1}\) contains only quantities that can be estimated in the real ocean (i.e., surface wind stress and seawater density along the coasts). Note that, in absence of topography, \(\psi _{1} = \psi _{2}\).

We herein follow the same approach to reconstruct the AMOC variability during the last millennium, although for the sake of completeness \(\psi _{2}\) is also evaluated and compared to the fully integrated AMOC streamfunction. This will allow assessing the validity of the approximations and assumptions used to reconstruct the AMOC and its variability, especially at multidecadal and longer timescales. Note that due to the lack of reliable surface wind paleorecords, AMOC reconstructions actually rely on the thermal-wind component only, and thus on the zonal density gradient calculated between the western and the eastern margins in the Atlantic Basin. Density is computed using the 1980 UNESCO International Equation of State (IES80; SCOR 1980) from the simulated salinity, temperature, and pressure at coastal points (Fig. 1). The same equation of state is used in the ocean component of ECHO–G (Legutke et al. 1999). Densities are thus computed across every longitude-depth section at the continental margins along America and Africa/Europe, but also across the Mid-Atlantic Ridge. The latter is nevertheless only present in the last two vertical levels of the model, below 3000 m depth (see Fig. 1b). Taking densities along the intersecting Mid-Atlantic Ridge might somewhat induce some error in the calculated thermal-wind, yet this must be of low relevance given the results shown in Sects. 4.1 and 4.2. The procedure followed here is, in fact, similar to that used for monitoring the AMOC within the RAPID project (e.g., Kanzow et al. 2007). The reconstructed AMOC is only represented from \(15^{\circ }\hbox {N}\) to \(\hbox {90}^{\circ }\hbox {N}\) because the assumption of geostrophy on which the thermal-wind equation is based fails in the Tropics. Since the Strait of Gibraltar is unrealistically broad at T30 resolution in ECHO–G, we mask the Mediterranean Sea from \(0^{\circ }\) to the west, thereby only considering points along the western coasts of Africa and Europe. Moreover, the Arctic Ocean is just considered between \(75^{\circ }\hbox {W}\) and \(30^{\circ }\hbox {E}\) and masked outside these longitudes.

3.2 The Florida Current

The FC forms the beginning of the Gulf Stream/North Atlantic Current system, extending from the Florida Strait to Cape Hatteras (e.g., Leaman et al. 1987). Its strength, of about 31 Sv at 26\(^{\circ }\)–\(27^{\circ }\hbox {N}\) (e.g., Leaman et al. 1987; Smeed et al. 2014) is mainly determined by the western limb of the wind-driven subtropical gyre, which accounts for about 21 Sv of the total transport at \(27^{\circ }\hbox {N}\) (Johns et al. 2002), and the strength of upper branch of the AMOC. In ECHO–G, however, the relatively coarse model resolution only allows for a simplified representation of the Florida Current (Fig. 1a). At \(27^{\circ }\hbox {N}\), the simulated FC is thus stronger (\(\sim\)36 Sv) and broader than in observations (about 600 and 80 km, respectively; Leaman et al. 1987), with a wind-driven meridional Sverdrup transport of about 25 Sv.

The FC is, to first order, in geostrophic balance; therefore, the vertical shear of its flow can be assumed proportional to the horizontal density gradient across its width. To reconstruct the FC variations during the last millennium, Lund et al. (2006) used two vertical profiles of foraminifera-based seawater densities. Samples were taken at two locations at each side of the flow: one at \(24^{\circ }\hbox {N}\)–\(83^{\circ }\hbox {W}\) and depths of 750, 540, 250 and 200 m, and the mixed layer, represented by the upper-50 m average; the other at \(24^{\circ }\hbox {N}\)–\(79^{\circ }\hbox {W}\) and depths of 700, 530 and 260 m, and the mixed layer. Lund et al. (2006) additionally assumed a level of no motion at 850 m depth. Since such profiles are at the same latitude and their maximum depth is well below 850 m (Fig. 1b), the FC meridional transport (FCT) can be approximated as:

$$FCT = - \frac{g}{\rho ^* f} \int ^{0} _{- 850} dz^{\prime} \int ^{z^{\prime}} _{-850} \left( \rho _e - \rho _w \right) \, dz^{\prime\prime}.$$
(10)

Note that Eq. (10) is a particular case of Eq. (2) but applied to the FCT. In the present study, the same approach as that of Lund et al. (2006) is followed; however, and due to the relatively low model resolution, the two locations used here for FC model-based reconstruction (hereafter, FC pseudo-reconstruction) are selected at each side of the simulated northward maximum transport: one at \(27^{\circ }\hbox {N}\)–\(78^{\circ }\hbox {W}\), and the other at \(27^{\circ }\hbox {N}\)–\(70^{\circ }\hbox {W}\) (triangles in Fig. 1a). The performance of the reconstruction technique, described in Sect. 4.3, is nevertheless relatively insensitive to this selection (not shown). The selected simulated density depths are the same as the stratigraphic levels in Lund et al. (2006) in order to make both results directly comparable.

4 Results

4.1 AMOC mean climatology

The simulated meridional overturning streamfunction of the Atlantic Ocean averaged over the last millennium (Fig. 2a) agrees well with the main structure of the AMOC according to observational estimates, as was described by Ortega et al. (2012). It is dominated by a clockwise cell with a maximum of about 18 Sv, comparable to hydrographic data (15 ± 2 Sv; Ganachaud and Wunsch 2000) and recent in-situ measurements at \(26^{\circ }\hbox {N}\) (18.5 ± 1.0; McCarthy et al. 2012). North Atlantic Deep Water (NADW) sinks around \(60^{\circ }\hbox {N}\) down to 5000 m depth, somewhat deeper than in observations (Talley 2003). The overflow over the Greenland–Iceland–Scotland ridge contributes to the AMOC with about 2 Sv, thus below the flux of about 6 Sv reported by Hansen and Østerhus (2000) and Olsen et al. (2008).

In comparison, the main clockwise gyre of the thermal-wind streamfunction based on Eqs. (6) and (7) (Fig. 2b, c) is confined south of the Greenland-Scotland ridge at shallower levels and is considerably stronger. The relatively large streamfunction values near the deep water formation region are partly counterbalanced by the other two circulation terms, namely the Ekman transport and the external mode (Fig. 2d, e).

The meridional Ekman transport (Fig. 2d), which is here assumed to take place in the upper 20 m, is characterized by a dipole circulation forced by the opposing predominant winds at tropical and mid-latitudes, which leads to downwelling around \(30^{\circ }\hbox {N}\). A similar pattern was associated with zonal wind changes by the East Atlantic Pattern in Ortega et al. (2012). In the Tropics, the southern Ekman cell exerts a stronger influence than the geostrophically balanced flow, with values above 12 Sv at the surface. Yet the thermal-wind component clearly dominates the main meridional circulation at mid and high latitudes.

The external mode (Fig. 2e) forms a clockwise overturning cell around \(30^{\circ }\hbox {N}\) with a maximum value of about 20 Sv at 1200 m depth. This term is thus more intense than both the Ekman transport and the thermal-wind component at low latitudes and contributes to a larger fraction of the AMOC there. Note that an approximation of this term is already contained in \(\psi _{tw_1}\).

The sum of the geostrophically balanced flow and the Ekman transport based on Eq. (6) (Fig. 2f) agrees with the target AMOC streamfunction. The resulting circulation is however shallower, and its maximum is shifted to the north, around \(50^{\circ }\hbox {N}\). In addition, its counterclockwise cell between \(55^{\circ }\) and \(60^{\circ }\hbox {N}\) is not present in the target AMOC. The sum of the circulation terms based on Eq. (7) (Fig. 2g) also broadly agrees with the simulated target AMOC, although the addition of the external mode results in a non-realistic higher maximum at subtropical latitudes that is not present in the simulated streamfunction. These discrepancies might result from the assumption that the level of no motion is at the ocean bottom for the thermal wind calculation, where simulated velocities are actually not negligible, together with the relatively coarse resolution of the model.

To summarize, both thermal-wind diagnostics capture the main features of the averaged circulation profile over the last millennium of the simulated AMOC in ECHO–G, with a lesser contribution from the external mode and the Ekman transport at northern latitudes.

4.2 AMOC variability

The variability of the AMOC is analyzed using the meridional overturning index (MOI), commonly defined as the AMOC maximum along the latitude-depth section (Delworth et al. 1993). Above decadal timescales, this index is known to capture well the large-scale AMOC variability (Vellinga and Wu 2004). However, the decomposition of this index in its thermal-wind, Ekman, and barotropic components is complicated by the fact that the location of the AMOC maximum varies in time (Ortega et al. 2012) and will therefore differ from one component to another (Fig. 2). Due to this limitation, the MOI is herein defined instead as the average between \(35^{\circ }\) and \(45^{\circ }\hbox {N}\) at 1500 m depth, where the climatological maximum of the simulated target AMOC is located (Fig. 2a). Analogous indices are defined for the thermal-wind, Ekman and external mode components (MTwI, MEkI, and MExI, respectively). Because our aim is to evaluate reconstruction methods that are applied to low-resolution proxy records (e.g., Lund et al. 2006) we focus on multi-decadal time scales. To this end, an 11-year running mean is previously applied to the simulated AMOC and its components.

The simulated MOI (Fig. 3, blue line), which ranges between 16 and 20 Sv, slightly weakens during the first 800 years at a rate of about 0.13 Sv/century. This decreasing trend accelerates to 1 Sv/century during the last two centuries, coinciding with the increase in GHG concentrations in the industrial period. A similar weakening in the AMOC strength has also been simulated in historical CMIP5 runs (e.g. Drijfhout et al. 2012; Cheng et al. 2013). There is no clear imprint of the solar irradiance or volcanic eruptions (Ortega et al. 2012), which are the dominant natural forcings at decadal and longer time scales (Fernández-Donado et al. 2013).

High agreement is observed between the variability of the MOI and of both MTwIs (black and red lines), whose major excursions and final trends are virtually identical, with correlation coefficients above 0.9 (with a significance level of \(\alpha <0.05\)). By contrast, the correlation coefficient between the MOI and the MEkI (orange line) is only 0.03, while with the MExI (magenta line) it is 0.67 (\(\alpha <0.05\)), thus reflecting no contribution to the variability by the Ekman component and some relatively minor influence from the external mode, with respect to that from the thermal wind. The indices for the sum of the terms given by Eqs. (8) and (9) (Fig. 3, dashed lines) do not completely capture the absolute magnitude of the target, the MOI, presumably because the thermal-wind components are just approximations; nevertheless, their variability is close to the target, as found for the thermal winds, because the latter component is the one that contributes the most to the AMOC variability at these time scales. Similar results to those described above are found for variability above decadal time scales of the overall streamfunctions: whereas both thermal-wind components capture relatively well the target variability, the Ekman transport and the external mode show negligible and minor contributions to the variability at these time scales, respectively (not shown).

Similar correlation coefficients are obtained between the MOI and the indices of the different components for the preindustrial period (1000–1800), when the strong trend due to changing anthropogenic GHGs is not present. We thus rule out that the high significant correlation arises exclusively from the common downward trends since 1800 onwards, and conclude that the thermal-wind term alone is able to capture most of the multi-decadal AMOC variability.

In order to evaluate how the different overturning indices represent the overall AMOC structure, we compute the regression of the simulated AMOC streamfunction on the standardized MOI, MTwI, MEkI, and MExI time series (Fig. 4). The regression of the AMOC against the MOI (Fig. 4a) shows an intensification of the NADW cell when the index increases, and vice versa. The largest regression values are located at \(50^{\circ }\hbox {N}\) at the same depths as the AMOC maximum, at 1500 m depth, within the centre of the shallowest maximum of variability (Fig. 4, contours). The simulated decreasing MOI during the industrial period therefore reflects a deceleration of the AMOC and a weakening of North Atlantic deep water formation.

A similar pattern is found for the regression of the AMOC against both MTwIs (Fig. 4b, c), with the significant coefficients mainly located within the NADW cell. The MTwIs thus capture the main variations of the simulated meridional overturning circulation. Note that the maximum values, as compared to the regression against the MOI, occur at deeper levels, between 2000 and 3000 m depth.

Significant regression coefficients between the AMOC streamfunction and the MEkI (Fig. 4d) are only found at the surface and around \(20^{\circ }\hbox {N}\), \(40^{\circ }\hbox {N}\), and \(60^{\circ }\hbox {N}\), locations where the Ekman cells are centered (Fig. 2d). This reflects the direct imprint of the surface wind-driven transport in the Atlantic gyres.

The regression of the AMOC streamfunction against MExI (Fig. 4e) is also associated with changes in the NADW cell, although with a more notable influence at the bottom of the basin between \(20^{\circ }\) and \(40^{\circ }\hbox {N}\), where the deepest maximum of variability is located (Fig. 4, contours). Thus, the variability of the MExI is mainly associated with circulation changes at tropical mid to deep basin levels.

To summarize, both the high correlations between the MOI and the thermal-wind indices as well as the similarities between the AMOC regression against these indices reveal that the thermal-wind component reproduces very closely the main MOI variations throughout the last millennium at multidecadal and longer timescales. These results support using the thermal wind to estimate the past variability of the Atlantic Ocean circulation on such time scales.

4.3 Florida Current variability

We now proceed to validate the reconstruction technique of the FCT using the simulation as a surrogate for reality. The contributions to the FC variability in the model are also assessed. Additionally, we discuss fluctuations in the FC in the context of changes in external forcing and internal climate variability. Finally, we compare the reconstructed FCT in the model with the proxy-based estimation of Lund et al. (2006).

4.3.1 Validation of the FC reconstruction technique in ECHO-G

The simulated FCT during the last millennium (hereafter \(\hbox {FC}_{S}\); Fig. 5a, green line) is computed in ECHO–G by integrating the meridional velocity field at \(27^{\circ }\hbox {N}\) between \(78^{\circ }\hbox {W}\) and \(70^{\circ }\hbox {W}\) (triangles in Fig. 1a), and from the level of no motion, at 850 m depth, to the surface. The time series is then smoothed with an 11-year running mean. \(\hbox {FC}_{S}\) shows similar decreasing trends to those in the MOI (Fig. 5a, blue line) over both the pre- and the industrial periods. The relatively high correlation between these two indices (0.66; \(\alpha <0.05\)) illustrates the close link between the FC and the AMOC strength, which was also suggested by De Coëtlogon et al. (2006) using AOGCM simulations, albeit of shorter length. The correlation increases up to 0.83 (\(\alpha <0.05\)) after smoothing with a 51-year running mean, due to the similar variability on centennial and longer time scales, and drops down to 0.47 (\(\alpha <0.05\)) if the time series are high-pass filtered.

The pseudo-reconstruction of \(\hbox {FC}_{S}\) proposed in Sect. 3.2 (hereafter \(\hbox {FC}_{\rho }\); Fig. 5a, red line) shows a high correlation of 0.9 (\(\alpha < 0.05\)) with the target (\(\hbox {FC}_{S}\)) and reproduces most of the simulated variability and long-term trends. Both FCTs thus capture a large part of the simulated MOI at multidecadal and longer timescales. Similar conclusions are drawn when only the preindustrial period is considered and, thereby, the simulated negative GHG-driven trend in the AMOC/FC strength excluded.

We now further assess the distribution of the FC with depth. At each level, the transport is calculated from the simulated meridional velocity field and compared with the corresponding pseudo-reconstruction from the simulated zonal density gradient (Fig. 6a, b). Both are filtered with a 51-year running mean to focus on the same frequency range as in Lund et al. (2006). The simulated and pseudo-reconstructed transports show the largest anomalies in the subsurface; these are larger and positive in the first century, remain predominantly positive until 1750, and finish with the most negative values in the whole millennium. The largest contribution to the total simulated transport results from levels below 300 m depth, whereas a minor contribution comes from shallower levels, as shown by the depth profile of its long-term mean (green line in Fig. 6e). This is relatively well captured by the pseudo-reconstruction (red line in Fig. 6e), which shows a maximum at about 400 m depth. The pseudo-reconstruction, on the other hand, tends to overestimate the simulated variability at the surface levels and to underestimate it below 300 m, as shown by their standard deviation depth profiles (red and green lines in Fig. 6d). Nonetheless, the pseudo-reconstruction captures well the simulated maximum contribution to the FC variability at around 400–600 m depth.

As already shown, \(\hbox {FC}_{\rho }\) does not exactly reproduce all the variability of \(\hbox {FC}_{S}\) (Fig. 5a). The temporal evolution of their difference, \(\hbox {FC}_{\rho }\)–\(\hbox {FC}_{S}\), illustrates when this disagreement is largest (Fig. 7a), being, in general, smaller than 1 Sv. This difference is significantly correlated with \(\hbox {FC}_{S}\) but not with the \(\hbox {FC}_{\rho }\) (correlation coefficients of 0.42 and −0.01, respectively; \(\alpha < 0.05\)). A relatively small part (18 % of the total variance) of the simulated FCT variability is therefore not captured by the reconstruction technique. The difference \(\hbox {FC}_{\rho }\)–\(\hbox {FC}_{S}\) shows a relatively low, but significant correlation with the wind-driven Sverdrup transport (0.24; \(\alpha < 0.05\)). Wind-driven currents can thus influence how well the reconstruction technique captures the total variability of the FCT. The difference \(\hbox {FC}_{\rho }\)–\(\hbox {FC}_{S}\) is also found to be related to simulated subsurface temperature changes to the southeast of the FCT maximum at \(25^{\circ }\hbox {N}\)–\(65^{\circ }\hbox {W}\) (not shown). \(\hbox {FC}_{\rho }\) tends to underestimate the strength of the \(\hbox {FC}_{S}\) transport (i.e., negative values in their difference) when temperatures there decrease at 250 m depth and increase simultaneously at 100 m depth, and vice versa. These temperature changes might thus induce an additional ageostrophic term in the FCT (Kontoyiannis and Watts 1990) that the reconstruction is unable to capture.

The difference between the simulated and the pseudo-reconstructed transport is largest at the deepest levels, but always smaller than 0.5 Sv, shifting toward more positive (negative) values above (below) 300 m in the last 100 years (Fig. 7b). This counterbalance between shallower and deeper levels is found at shorter time scales during the entire period as well. The root-mean-square error of \(\hbox {FC}_{\rho }\) with respect to \(\hbox {FC}_{S}\) (Fig. 7c) shows the largest disagreement at the deepest levels, where the underestimation of the variability of \(\hbox {FC}_{S}\) by \(\hbox {FC}_{\rho }\) is also largest (Fig. 6d).

The above-described results illustrate that the pseudo-reconstruction of the FCT in the model reproduces most of the simulated total variability in the area, although local variations can influence its accuracy.

4.3.2 Contributions to the FC variability

To further investigate what drives the FCT variability in ECHO–G, we compare the separate contributions from temperature and salinity to the vertically resolved pseudo-reconstructed transport (Fig. 8). Anomalies of the thermally-driven transport (Fig. 8a) tend to be positive from 1000 to 1400 and negative thereafter, especially after 1650 when negative transport anomalies extend to deeper levels. The haline contribution (Fig. 8b) tends to oppose that from the thermal one, as evidenced by their negative correlation coefficients at all depths (Fig. 8c). Similar phenomena of density compensation have been reported in ocean fronts (Tippins and Tomczak 2003) and within the limits of the mixed-layer depth (Rooth 1993; Rudnick and Ferrari 1999). Thus, the vertical distribution of the pseudo-reconstructed FCT changes in ECHO–G are mostly explained through changes in the zonal temperature gradients, with opposing changes in zonal salinity gradients.

The latter result is confirmed for the total (vertically-integrated) FCT of both contributions (Fig. 9a): the thermally and haline-driven total transports are found to fluctuate in phase opposition at shorter (interdecadal) timescales, with a correlation coefficient of −0.67 (\(\alpha < 0.05\)), as observed in the density-depth profiles (Fig. 8). As anticipated, the thermally-driven variability of the FCT strongly resembles that of \(\hbox {FC}_{\rho }\) (Fig. 9b, red line), with a correlation coefficient of 0.9 (\(\alpha < 0.05\)). Besides, largest excursions in \(\hbox {FC}_{\rho }\) in 1020, 1340, or 1460 are only found in the thermally-driven FCT. Only the thermal contribution (magenta line) weakens throughout the last millennium, including the larger negative trend in the final two centuries. The haline term (Fig. 9a, blue line), by contrast, shows a much weaker correlation with \(\hbox {FC}_{\rho }\) (−0.36; \(\alpha < 0.05\)) since this term predominantly fluctuates on longer time scales (multidecadal and centennial). It does not capture well high frequency (decadal) FC variability either and shows no clear trend at the end of the last millennium. Note that the sum of these two components closely resembles \(\hbox {FC}_{\rho }\) (Fig. 9b); minor discrepancies between them possibly derive from the assumption of a linear equation of state for density, which is estimated by just adding the haline and thermal contributions. We thus conclude that in ECHO-G it is the thermal contribution to the zonal density gradient that mainly drives variations in the FC strength (Fig. 9a). Note that, for both the thermal and haline contributions, the integrated transport variability is primarily determined by changes at mid-levels (50–500 m), where the largest thermal- and haline-driven anomalies are found (Fig. 8). A minor contribution, in contrast, comes from changes at the surface. These results agree with the standard deviation profile of \(\hbox {FC}_{\rho }\) (Fig. 6c).

From the above-described results, it is also possible that the thermal component of the FCT is dominating the \(\hbox {FC}_{\rho }\) anomalies as a result of an enhanced temperature variability range due to direct radiative effect from the external forcings. To reject this hypothesis, we further analyze the thermally and haline-driven contributions of the FC in a control simulation with ECHO–G (Ortega et al. 2012). Here, \(\hbox {FC}_{\rho }\) and the thermal component are also highly correlated (coefficient of 0.8; \(\alpha < 0.05\)), whereas the haline contribution shows no significant correlation (not shown). In the absence of external forcing, density changes that determine the FC variability in ECHO–G are therefore mostly driven by temperature variations, as concluded for the forced run.

4.3.3 Externally-forced influences on the FC variability

We next explore whether the variability of the simulated FCT in ECHO–G is correlated with changes in the external forcing, as well as to what extent the reconstruction procedure is able to capture this relationship. This analysis assumes a direct and linear response of the FC/AMOC to external forcing, which could nevertheless overlook more complex responses; a more detailed study of these was already conducted by Ortega et al. (2012) and is beyond the scope of this paper. We consider all the external forcings used in the simulation: volcanic eruptions and changes in total solar irradiance and GHG concentrations (see Sect. 2 for a detailed description), as well as their linear combination (hereafter, total forcing) following a similar approach as in Fernández-Donado et al. (2013). For a coherent comparison with time series of Fig. 5a, the external forcings are smoothed with an 11-year running mean filter.

\(\hbox {FC}_{S}\), \(\hbox {FC}_{\rho }\) and the MOI are significantly correlated with the total forcing over the entire past millennium (\(-0.63\), \(-0.57\), and \(-0.44\), respectively; \(\alpha < 0.05\)). These correlations are however mainly explained by the trend in the forcing during the last two centuries, because no significant correlations are obtained over the preindustrial period (1000–1800; not shown). This trend is, in turn, mostly explained by the changes in GHG concentration and total solar irradiance, the only individual forcings that correlate significantly with \(\hbox {FC}_{S}\) and \(\hbox {FC}_{\rho }\) over the entire millennium (\(-0.61\), and \(-0.45\), respectively for \(\hbox {FC}_{S}\); \(-0.55\), and \(-0.39\), respectively for \(\hbox {FC}_{\rho }\); \(\alpha < 0.05\)). By contrast, the MOI only shows significant correlations with the GHG concentration changes (\(-0.44\); \(\alpha < 0.05\)). Note that, however, all correlations with the individual forcings become non-significant if only the preindustrial period is considered. Simulated fluctuations in the AMOC, and hence in the FC, can thus be attributed to internal climate variability during this period or potentially non-linear response to the forcing. This result is not obvious, given the strong contribution of temperature to \(\hbox {FC}_{\rho }\), and has important implications for the comparison against proxy data in the following Subsection. Interestingly, correlation coefficients between the \(\hbox {FC}_{S}\) or \(\hbox {FC}_{\rho }\) with the forcing(s) are comparable, albeit larger for the former; this suggests that this reconstruction technique can reproduce without overestimating most of the AMOC/FC variability during the last millennium, regardless of whether this is internally generated or externally forced. Similar coefficients (not shown) to those for the FC are found between the external forcings and its thermal component; by contrast, the haline contribution does not significantly correlate with any forcing. Only the thermal component is thus able to capture their effect on the trends during the industrial period. The largest contribution to the correlations between \(\hbox {FC}_{S}\), or \(\hbox {FC}_{\rho }\), and the external forcing comes from the deepest level (800 m). Correlations decrease toward the surface, where they become negligible (not shown). Ocean-atmosphere flux variability therefore tends to dominate over any externally-forced signal in the upper ocean.

4.3.4 Model-data comparison

We finally compare our pseudo-reconstruction of the FCT in ECHO–G, \(\hbox {FC}_{\rho }\), with the proxy-based reconstruction of Lund et al. (2006). The latter (hereafter \(\hbox {FC}_{R}\); Fig. 5b, black line), which captures only variability on centennial timescales, shows a weakening of about 2 Sv from 1000 to 1800. In this period the temporal evolution of \(\hbox {FC}_{\rho }\) (Fig. 5b, red line) is consistent with that of \(\hbox {FC}_{R}\), lying within its 95 % confidence interval. However, while \(\hbox {FC}_{R}\) increases during the last 200 years, recovering close to its initial values, both the FC and the AMOC strength (represented by the MOI) in the model continue weakening toward their minimum with an even stronger decreasing trend. We anticipate that part of the model-data disagreement can be understood as a result of the dominance of the internal climate variability (see Sect. 5).

The vertically resolved reconstructed FCT (Fig. 6c) agrees with our model results in showing positive anomalies during the first years. However, they show a stark contrast from 1200 onward: while the reconstructed anomalies are predominantly negative between 1200 and 1850 and turn positive in the industrial period, transport anomalies in ECHO-G fluctuate between positive and negative values for the period 1100–1750 and become negative after 1800. The reconstructed long-term mean FCT also differs from the simulated results (Fig. 6e): for the former, the largest contribution comes from the shallowest levels, above 300–400 m, whereas for the latter, the largest contribution is given by mid levels, between 300 and 600 m. The FCT reconstruction shows its largest standard deviation values at the surface (Fig. 6d, black line), with a second relative maximum at mid-levels, similarly to the pseudo-reconstruction. Thus, if we keep in mind the bias in the representation of the FCT pseudo-reconstruction with respect to the simulation, our perfect model assessment suggests that the FC reconstruction might be overestimating the real variability near the surface, and underestimating it at the deepest levels.

Seawater density changes at each location are the basis for the calculation of the zonal density gradient and, therefore, of the FCT. Analyzing their variability in the model thus yields information about how they contribute to the flow, and allows the comparison with the available density reconstructions of Lund et al. (2006). At \(78^{\circ }\hbox {W}\) and \(70^{\circ }\hbox {W}\), the density range is similar in the model (Fig. 10) and in the reconstruction (see their Fig. 2), between 24 and 27 kg/m\(^3\). However, there is disagreement in the simulated and reconstructed long-term trends. According to Lund et al. (2006), the weakening of \(\hbox {FC}_{R}\) for the period 1000–1800 results from a slight decreasing trend in densities at mid-levels (between 200 and 540 m depth) at the western location (Dry Tortugas), while density to the east remains nearly constant. From 1800 onward, \(\hbox {FC}_{R}\) increases in response to lighter densities below the mixed layer at the eastern location (Great Bahama Bank). The simulated seawater densities show an overall increasing trend throughout the entire millennium at both locations, of about 0.01 kg/m\(^3\)/century (Fig. 10). The weakening in \(\hbox {FC}_{\rho }\) between 1000 and 1800 (Fig. 5b, red line) results from a reduction in the zonal density gradient caused by eastern densities increasing faster than the western ones. The larger weakening of \(\hbox {FC}_{\rho }\) during the last 200 years, on the other hand, is related to a slight reversal in the long-term trend of the western densities in the upper levels, controlled by changes in temperature (not shown), which further debilitates the east–west density gradient. This might be related to larger thermal sensitivity to the GHGs in the inner region of the Gulf of Mexico than East of Florida, as previously evidenced in this simulation (Ortega et al. 2013). Interestingly, other short-lived intense \(\hbox {FC}_{\rho }\) excursions, like those in ca. 1100 and 1480, seem to be explained by simulated local density extremes at the eastern location at \(70^{\circ }\hbox {W}\) at the deepest levels, of 530 and 710 m depth. In the model, the differential density trends at both sites are mostly driven by the simulated temperature anomalies throughout the last millennium. However, the vertical temperature and salinity distributions at either site are not available for the reconstructions for a more detailed comparison.

We finally assess whether comparable trends exist in the aforementioned control simulation with ECHO-G that could point to a model drift effect. In this case, densities at the same sites show no trend at the uppermost levels, while densities at mid and low levels increase throughout the simulation, although with a smaller trend, of about 0.005 kg/m\(^3\)/century, than that of the forced run (not shown). Nevertheless, because these trends are constant and relatively similar at the western and eastern locations, the FCT in the control run does not show any long-term trend (not shown). The above-described trend reversal in western densities at mid-levels results from the increasing GHG concentrations over the industrial era in the externally-forced simulation, which reflects the imprint of this forcing on the seawater density field. Additionally, the different trends in the eastern and western densities that lead to the long-term weakening in \(\hbox {FC}_{R}\) could somehow be driven by the external forcing during the preindustrial era. Yet the simulated FC and AMOC do not show significant correlations with any external forcing for this period (see Sect. 4.3.3).

5 Discussion and conclusions

This work aimed at verifying the estimation of the meridional thermal-wind transport calculated from zonal density gradients as an accurate technique for reconstructing the variability of the ocean circulation. In particular, we compared the simulated and the pseudo-reconstructed AMOC and FC above multidecadal timescales during the last millennium in a transient climate simulation with the ECHO-G model that includes the relevant forcings for this period: GHGs, solar irradiance, and volcanism.

Firstly, we decomposed the AMOC into three components, the thermal-wind shear, the Ekman transport, and the external mode, in order to explore agreement between the variability of the thermal-wind term and that of the simulated target AMOC. The mean meridional streamfunction of the thermal-wind component over the last millennium agrees well with the simulated AMOC target, especially at mid and high latitudes of the NH. This subscribes conclusions from previous studies for shorter and/or idealized experiments (e.g., Hirschi and Marotzke 2007; Köhl and Stammer 2008; Josey et al. 2009). Besides, the MOI variability is also well captured by the thermal-wind component, whereas the Ekman transport and the external mode only contribute marginally. The regression of the simulated AMOC against the standardized thermal-wind index is as well found in broad agreement with that against the MOI and suggests that the thermal-wind component also captures the connection between fluctuations in the MOI and the NADW cell. Based on our result we conclude that both the North Atlantic circulation structure and its variability on multidecadal and longer timescales can be well reproduced from reconstructions based on coastal zonal density gradients. This conclusion extends previous results at interannual time scales (Hirschi and Marotzke 2007) with a sparse density data set subject to uncertainties (Hirschi and Lynch-Stieglitz 2006). The good performance of this reconstruction technique can furthermore be explained through a connection between variations in the AMOC strength and those of the transport of temperature and/or salinity anomalies along the Atlantic Ocean, since the latter influences the zonal density gradient and, thereby, the thermal-wind component. It is worth noting that simulated changes in the AMOC strength on these timescales in ECHO-G are ultimately wind-driven by a NAO-like pattern (Ortega et al. 2012). Fluctuations in deep water formation propagate along the oceanic basin, affecting the AMOC and modifying its zonal density gradient. This atmospheric-induced variability is thereby captured by the thermal-wind component of the AMOC.

We further evaluated in ECHO-G the method for reconstructing the FC strength and whether it reproduced the simulated FC and AMOC variability. Both the simulated and pseudo-reconstructed FC show the same long-term trends and variability over the last millennium. They are moreover closely linked with the simulated variability of the AMOC strength. Although a similar connection has been observed for shorter simulations (about 50 years; De Coëtlogon et al. 2006), it is the first time this is illustrated at longer time scales, and, in particular, in a last-millennium transient simulation. This result supports inferences of the AMOC variability through reconstructions of the FC intensity with paleo-estimations of the density field. Variability at depth for both FC representations in ECHO-G was found closely related as well, although the reconstruction procedure tends to underestimate the amplitude of the variability at the deepest levels. When the relative contribution of the thermally and haline-driven FC were compared, we found that the thermal contribution alone can explain most of the FC variability, with salinity changes generally opposing, and thereby moderating, temperature-driven fluctuations.

The analysis of the correlations between the simulated AMOC/FC strength and the external forcings suggested that internal climate variability mainly drives the AMOC/FC variations during the preindustrial era (1000–1800), whereas increasing GHG concentrations forces a weakening of both the FC and the AMOC during the industrial era (1800 onward). These features are well captured by the pseudo-reconstructed FC and, in particular, by its thermal component.

The comparison with the reconstructed FC from Lund et al. (2006) clearly yielded different results for the preindustrial and industrial periods. ECHO-G simulates a FC decrease during the first 800 years due to a weakening in the zonal density gradient that is consistent with a reconstructed FC weakening. However, in ECHO-G such weakening is explained by the relative increase in the eastern densities with respect to the western ones, with both sites showing increasing trends, while in the proxy-based reconstruction it is driven by the decrease in the western densities at mid-levels, with the eastern ones remaining nearly constant. Nonetheless, since the simulated variability of the FC/AMOC during the preindustrial period is not related to any external forcing but to internal climate variability, agreement between the model and the reconstruction is not expected for this period.

As to the industrial period, ECHO-G features a decreasing FC/AMOC trend instead of the recovery shown by the proxy-data. This discrepancy is more difficult to explain and points to different driving mechanisms in the models from those suggested for the real world. From the model side, several possible hypotheses can be proposed to explain such results: first, according to Lund et al. (2006) changes in the FC reconstruction were modulated by salinity anomalies in the tropical North Atlantic linked to a shift in the latitudinal position of the Atlantic Inter-Tropical Convergence Zone (ITCZ) that, ultimately, could have contributed to the recovery of the AMOC after the LIA. ECHO-G does not support this mechanism, however. Whereas an increase in the sea-surface salinity is simulated in the Tropical Atlantic by the end of the past millennium associated with reduced precipitation and enhanced evaporation (Ortega et al. 2012), likely reflecting a southward shift in the ITCZ position, the resultant salinity anomalies are not eventually transported toward regions of oceanic deep mixing in the North Atlantic and, thereby, do not contribute to strengthen the AMOC. Such a particular mechanism, nonetheless, might not be well resolved in this model due to a misrepresentation of potentially important processes influencing the hydrological cycle.

Another possibility is that reconstructed FC changes during and after the LIA were associated with shifts in the wind-driven subtropical gyre (Lund et al. 2006), itself driven by changes in the ITCZ position. ECHO-G does not support this hypothesis either. The simulated wind-driven component of the FC shows no long-term trend similar to that in total FCT. This suggests that the FC weakening is entirely associated with the reduction in the AMOC strength over the industrial period.

In addition, trends of opposite sign in the proxy-based FC during the nineteenth and twentieth centuries could also indicate that some of the response mechanisms to the external forcing are not well represented in ECHO-G. Alternatively, the FC recovery could instead be associated with internal climate variability in the Atlantic Ocean that the simulation is not expected to reproduce. These results thus raise the question of what is the real, ultimate driver of the FCT on decadal and longer time scales and prompt for further analyses with other simulations and data sets.

Finally, due to the relatively low model resolution, sites where densities were collected for the pseudo- and proxy-based FC reconstructions do not exactly coincide. Mechanisms driving density anomalies at both FC locations could thus differ between both setups, and further explain the model-data discrepancy in the past two centuries. On one hand, the reconstructed increase in cross-current density gradient and, hence, in the FC, was driven by decreasing densities at the Great Bahama Bank core sites (i.e. at western location) (Lund et al. 2006), in turn caused by upper-ocean warming and/or freshening. On the other hand, ECHO-G simulates warmer and saltier surface conditions in the Tropical Atlantic during the industrial period than during the period before, although it is the the differential change in temperatures that causes the overall weakening of the cross-current density gradient (Sect. 4.3.2). The origin in recent trends in the pseudo- and proxy-based FCT therefore results from different mechanisms. Model results are nonetheless rather insensitive to the selection of the density points at both sides of the Florida Current (not shown); indeed, a setup closer to the one used in Lund et al. (2006) does not substantially modify our conclusions.

In this line, an important caveat must nevertheless be stressed when interpreting the results of the present study. Due to the relatively low model resolution, ECHO-G presents a simplified bathymetry of the actual one in the area of the Caribbean Sea and Florida Current, unable to resolve the Caribbean Passages, Cuba and the Great Bahama Banks, which form the Florida Strait, nor the complex system of upper-ocean currents within this region. Our validation of the reconstruction technique is therefore only applicable within the limits of such simplified model reality. Caution must then be advised in not over-interpreting our model results, for it is still uncertain whether the reconstruction technique is also valid for a finer, hence more realistic, resolution. It is however worth stressing that previous studies found ECHO-G realistic enough to capture changes on a wide range of spatio-temporal scales (e.g., Beltrami et al. 2006; González-Rouco et al. 2006; Ortega et al. 2012, 2013; Fernández-Donado et al. 2013). It is though unlikely that trends at the spatial scale of the FC depend on the very local character of the bathymetry and spatial resolution in the model. If at all, this would be an argument to render unreliable the reconstructed variability. Earlier modeling studies on the variability of the FC and the AMOC with higher resolution models, albeit for different climate scenarios, supported the idea of covariability between the FC and the AMOC (e.g., Johns et al. 2002; De Coëtlogon et al. 2006; Lynch-Stieglitz et al. 2014).

From the data side, it may also be possible that the differences stem from problems in the proxy data or reconstructions. A first plausible source of error is the relationship between foraminiferal \(\delta ^{18}\hbox {O}\) and seawater density, which is assumed to be constant during the past millennium. Lund et al. (2006) assessed how a cooler high latitude North Atlantic would influence such relationship, thus finding that only very large and unrealistic shifts in the \(\delta ^{18}\hbox {O}_{sw}\) would change the density calibration enough to account for the 2–3 Sv signal in the FC during past millennium. Nevertheless, high resolution reconstructions are nonexistent yet, so these assumptions have not been verified with actual measurements. Another source of error is associated with the stratigraphy and the age models of the cores used for the reconstruction. These cores are characterized for a high sedimentation rates and well-constrained age models. Yet bioturbation can smooth any changes in \(\delta ^{18}\hbox {O}\) and, therefore, in the estimated seawater density time series. Because of this and the uncertainties in the age models, Lund et al. (2006) smoothed their records with 100-year wide windows. Short-term changes, for example during the twentieth century, would thus be difficult to detect with these sediment records. It thus occurs that, although the Bahama Bank cores show a general trend towards lower \(\delta ^{18}\hbox {O}\) values during the last 200 years, they also show a reversal of this trend in the uppermost part (Figure  2d in Lund et al. 2006), which possibly indicates a real trend toward greater densities during the past 50–100 years. Such a reversal trend would reduce the cross-current density gradient and, hence, the reconstructed FC strength, more in agreement with our model results. Last but not least, there might be some error in the assumptions taken in Lund et al. (2006) to interpolate missing data of the Dry Tortuga cores.

Other recent reconstruction and modeling studies can help shed light on the proxy-model discrepancy during the industrial period. In agreement with the FC proxy, a decrease of the surface AMOC over the last millennium with a recovery after 1945 was reconstructed from radiocarbon contents in sea-shells from the North Icelandic shelf (Wanamaker et al. 2012). As pointed out by Rahmstorf et al. (2015), however, this latter recovery contradicts recent instrumental results that describe an AMOC weakening since 1930 (Dima and Lohmann 2010). The index for the AMOC strength proposed by Rahmstorf et al. (2015) from reconstructed surface temperature during the last millennium features an unprecedented AMOC decline over the industrial period, with minimum values at end of the twentieth century. When applied in ECHO-G, this index reproduces most of the variability of the MOI, including the prominent fall during the last two centuries (not shown). Many CMIP5 historical simulations also describe a reduction of the AMOC strength during the last century (Drijfhout et al. 2012; Cheng et al. 2013), similarly to ECHO-G. Jungclaus et al. (2014) explored this modern AMOC reduction in the context of the last millennium using an ensemble of transient Earth-System-model simulations, thus relating the weakening of the AMOC at subtropical and subpolar latitudes to reduced convection activity in the Labrador Sea. These authors also described a strengthening of the subpolar gyre during the twentieth century, which, we conjecture, might explain the recovery found in the surface AMOC reconstruction off North Iceland by Wanamaker et al. (2012); however, the latter could also have been prominently influenced by changes in the Arctic circulation, wind patterns, and/or local currents. The majority of these studies thus would be consistent results shown in ECHO-G for the industrial period. We believe, nevertheless, that new reconstructions and modeling studies are necessary in order to clarify this point, rejecting or confirming conclusions drawn in existing works.

The analysis carried out in the present study followed, as a first step, a perfect-model approach that does not account in the model world for potential uncertainties in the reconstruction. A detailed analysis of the impact of reconstruction uncertainties on the performance of the reconstruction technique was however not performed and is beyond the scope of this work. Recent improvements in climate modeling make also necessary future analyses with different AOGCMs and Earth System models for an assessment of resemblances with present reconstructions and robustness of our results. Higher resolution in climate models will allow for a more direct model-data comparison over the same regions, whereas a deeper understanding of internal climate variability and its representation in climate models (Reichler and Kim 2008) could moreover help to disentangle the variability shown in the reconstructions.

In conclusion, we here provide palpable evidence that, first, the AMOC spatio-temporal variability above multidecadal timescales during the last millennium can be well reproduced by reconstructions based on its thermal-wind transport, and, second, that changes in the AMOC strength on such time scales can be well traced through the reconstructed Florida Current volume transport and, particularly, by its thermally-driven component. In addition, our results do not support the reconstructed recovery of the Florida Current strength suggested by the proxy after the Little Ice Age; the simulated Atlantic Ocean circulation, instead, largely weakens during the industrial period, reaching a minimum strength by the end of the past millennium.

Fig. 1
figure 1

The North Atlantic Ocean in ECHO-G: a margins (dark gray) between \(15^{\circ }\)–\(90^{\circ }\hbox {N}\) at the surface, and b bottom topography at \(45^{\circ }\hbox {N}\) (upper panel) and \(27^{\circ }\hbox {N}\) (lower panel). Grid points where densities are selected to calculate the thermal-wind component of the AMOC are shown in light gray. Those used to calculate the FCT are shown by the two triangles in a (see Sect. 3.2). Also in a average oceanic circulation simulated in ECHO-G for the upper 800 m (arrows, in m/s)

Full size image
Fig. 2
figure 2

a AMOC streamfunction (in Sv) between \(15^{\circ }\) and \(90^{\circ }\hbox {N}\) and averaged for the period 1000–1990 simulated by ECHO-G. b–e AMOC contributions related to: b the thermal-wind from Eq. (6); c and from Eq. (7); d the Ekman transport; e the external mode. Sum of terms defined by f Eq. (8), and g by Eq. (9). Contour interval is 2 Sv

Full size image
Fig. 3
figure 3

Indices for the strength of the simulated AMOC (in Sv) and its components in ECHO-G for the period 1000–1990, defined as the average of their respective streamfunction (Fig. 2) between \(35^{\circ }\) and \(45^{\circ }\hbox {N}\) at 1500 m depth and smoothed with an 11-year running mean: meridional overturning index (MOI, blue line) and contributions for the thermal-wind from Eq. (6) (MTwI\(_1\), black line) and Eq. (7) (MTwI\(_2\), red line), for the Ekman transport (MEkI, orange line), the external mode (MExI, pink line), and the sum of terms defined by Eq. (8) (MI\(_1\), dash black line) and Eq. (9) (MI\(_2\), dash red line)

Full size image
Fig. 4
figure 4

Regression coefficients (in Sv; color shading) for the simulated AMOC (Fig. 2a) against the standardized a MOI, b MTwI\(_1\), c MTwI\(_2\), d MEkI, and e MExI, from Fig. 3. Only values statistically significant below the 0.05 level are shown. Contours account for the simulated AMOC standard deviation in increments of 0.2 Sv. Standardization is done by dividing the long-term anomalies of each index by its standard deviation over the full period. Calculations are performed after smoothing both the indices and the simulated AMOC with an 11-year running mean

Full size image
Fig. 5
figure 5

a FCT anomalies (in Sv) at \(27^{\circ }\hbox {N}\) as simulated by integrating the meridional velocity fields from 850 m depth to the surface (\(\hbox {FC}_{S}\), green line), and pseudo-reconstructed (\(\hbox {FC}_{\rho }\), red line) from Eq. (10) using the density difference between \(70^{\circ }\hbox {W}\) and \(78^{\circ }\hbox {W}\) at the specific levels described in Sect. 3.2. Simulated MOI anomalies (Fig. 3) are also shown for comparison (blue line). Anomalies are calculated with respect to the average between 1000 and 1990 and smoothed with an 11-year running mean. b Estimated FCT anomalies (in Sv) by Lund et al. (2006) (\(\hbox {FC}_{R}\), solid black line), with their 95 % confidence limits (black dashed lines). \(\hbox {FC}_{\rho }\) anomalies are also shown here smoothed with a 51-year running mean for a better representation of the low-frequency (red line)

Full size image
Fig. 6
figure 6

Vertical distribution of the a simulated, and b pseudo-reconstructed FCT (in Sv). Anomalies are calculated with respect to the average between 1000 and 1990 and smoothed with a 51-year running mean for a better representation of the low-frequency. c Reconstructed FCT anomalies (in Sv) versus depth by Lund et al. (2006). d Standard deviation (STDV; in Sv) and e mean (in Sv) versus depth of FCTs depicted in panel a (green line), b (red line), and c (black line)

Full size image
Fig. 7
figure 7

Time evolution of the error in the estimation of \(\hbox {FC}_{S}\) by \(\hbox {FC}_{\rho }\) (in Sv) a from Fig. 5a, and b vertically from Fig. 6, both after smoothing with an 11-year running mean. c Depth profile of the root-mean-square error (in Sv) calculated from panel b

Full size image
Fig. 8
figure 8

Depth profile of the a temperature and b haline contribution anomalies (in kg/m\(^3\)) to the pseudo-reconstructed FCT (Fig. 6b). Anomalies are calculated and smoothed like in Fig. 6. c Correlation coefficient between both profiles, versus depth. Significant values (\(\alpha < 0.05\)) are depicted by solid line

Full size image
Fig. 9
figure 9

a FCT anomalies (in Sv) calculated from the salinity (blue line) and temperature (pink line) contributions to the density gradient between (\(27^{\circ }\hbox {N}\), \(70^{\circ }\hbox {W}\)) and (\(27^{\circ }\hbox {N}\), \(78^{\circ }\hbox {W}\)) using the specific levels described in Sect. 3.2 and integrated for the upper 850 m. b Linear sum of both contributions (black line), compared with the \(\hbox {FC}_{\rho }\) (red line; from Fig. 5). Anomalies are calculated and smoothed like in Fig. 5

Full size image
Fig. 10
figure 10

Values of \(\sigma _t\) (in kg/m\(^3\)) between 1000 and 1990 at a (\(27^{\circ }\hbox {N}\), \(78^{\circ }\hbox {W}\)), and b (27\(^{\circ }\hbox {N}\), \(70^{\circ }\hbox {W}\)), for the specific levels described in Sect. 3.2 and used in Lund et al. (2006). Note that largest densities are at the bottom. Gray/solid-black lines are smoothed values after an 11-/51-year running mean, respectively. For these latter: c vertical distribution of the trends (in \(\hbox {kg/m}^{3}\)/century) between 1000 and 1800 (black line), and between 1800 and 1990 (gray line) at \(70^{\circ }\hbox {W}\) (dashed lines) and \(78^{\circ }\hbox {W}\) (solid lines). Values are also indicated together with the level in a and b, respectively for each period

Full size image