Water Mass Transformation in the North Atlantic and Its Impact on the Meridional Circulation: Insights from an Ocean Model Forced by NCEP–NCAR Reanalysis Surface Fluxes

Decadal-scale climate variability in the North Atlantic thermohaline circulation is simulated using a sigma- coordinate primitive equation model, forced by NCEP–NCAR reanalysis surface forcing ﬁelds for the period from 1958 to 1997. Surface heat and freshwater ﬂux are expressed in terms of surface thermal and haline density inputs, diagnosed by the model. Variability in surface density ﬂuxes is closely correlated with the North Atlantic Oscillation and demonstrates differences with the original surface heat and freshwater ﬂuxes. Leading modes of surface water mass transformation are considered in the T – S plane. They identify decadal-scale variability associated with the transformation of the Labrador Sea Waters and Subtropical Mode Waters. Analysis of the model responses to the surface forcing shows an immediate reaction of meridional heat transport to the wind stress curl, resulting in a decrease of meridional heat transport at 48 8 N and an increase in the subtropics. Delayed baroclinic responses to the surface heat forcing are identiﬁed at time lags of 3 and 7 yr. The 3-yr response is represented by an increase in the total meridional heat transport in subpolar latitudes and its simultaneous increase in the Tropics and midlatitudes. The 7-yr delayed response to the surface heat forcing is associated with the strengthening of meridional heat transport at all latitudes. However, 7-yr responses may be inﬂuenced by the self-correlation in the meridional heat transport and forcing function. Meridional overturning is largely responsible for the variability observed, demonstrating high correlation with meridional heat transport.


Introduction
The North Atlantic plays an important role in global climate variations, being largely responsible for the intensity of interocean exchanges and meridional overturning. The leading mode of the atmospheric variability over the North Atlantic is represented by the North Atlantic Oscillation (NAO; Rogers 1984; Barnston and Livezey 1987;Hurrell 1995), which exhibits close correlation with the sea surface temperature (SST) and surface fluxes on interannual to decadal time-scales, as shown in a large number of studies (Bjerknes 1964;Wallace et al. 1990;Cayan 1992a,b;Deser and Blackmon 1993;Kushnir 1994;Iwasaka and Wallace 1995;Halliwell and Mayer 1996;Halliwell 1997;Peng and Fyfe 1996;Hakkinen 2000;Gulev et al. 2002). Thus, the NAO certainly has signatures in the deep ocean circulation driven by the surface fluxes of heat, freshwater, and momentum.
Observational evidence of NAO-related changes in the temperature and salinity of the North Atlantic at intermediate depths has indeed been found in historical hydrographic data (Levitus 1989;Levitus and Antonov 1995). In particular, analysis of long time series (Lazier 1980(Lazier , 1995Levitus et al. 1995;Joyce and Robbins 1996) has shown that decadal variability in the upper and intermediate waters of the North Atlantic shows opposite tendencies in the subpolar and subtropical gyres. Dickson et al. (1988Dickson et al. ( , 1996 and Curry et al. (1998) found that generation and propagation of the Labrador Sea Water (LSW) and Subtropical Mode Waters (STMW) are key indicators of the intensity of meridional overturning in the North Atlantic. LSW properties have changed considerably over the last three de-cades. Its pathways were first described by Talley and McCartney (1982) and have since been the subject of detailed inspections (Pickart 1992;Cunningham and Haine 1995a,b;Sy et al. 1997;Molinary et al. 1998). Joyce et al. (2000) found that the STMW potential vorticity is highly correlated with the NAO and the surface buoyancy forcing. Mode waters (i.e., waters preferably generated by air-sea interactions) represent a large volume of water masses in the ocean and keep the signatures of the year-to-year variations in air-sea exchanges. Since regions of the mode waters formation in the North Atlantic are well known, changes in their properties (in particular, the formation rates) may serve as the indicators of the ocean climate change. One objective of our study is to investigate the projections of the NAO buoyancy forcing onto formation rates of the North Atlantic Mode Waters.
Observational evidence of the changes in meridional heat transport (MHT) is available from the analysis of full-depth oceanic cross sections (Hall and Bryden 1982;Roemmich and Wunsch 1984;Read and Gould 1992;Parrilla et al. 1994;Koltermann et al. 1999). The MHT changes are assumed to be associated with the intensity of the meridional overturning cell. However, such MHT estimates are available for selected years only. Moreover, they are largely influenced by uncertainties in estimation of MHT from hydrographic data and by sampling problems at the oceanic cross sections. Moreover, it is difficult to relate with certitude variability in the transport at a given oceanic section with the basin-scale changes in the atmospheric forcing. Thus, analysis of mechanisms through which the NAO impacts changes in meridional overturning circulation is a challenging task for coupled and oceanic GCMs.
Interdecadal variability of the thermohaline circulation was studied in a number of coupled GCMs. Delworth et al. (1993), using the Geophysical Fluid Dynamics Laboratory (GFDL) model, found interdecadal variability of the thermohaline circulation, associated with density forcing. Delworth and Greatbatch (2000) showed that interdecadal variations in the intensity of the thermohaline circulation can be explained by the ocean response to the atmospheric forcing. By performing several experiments with the oceanic component of the coupled model, forced by surface fluxes, they found that the role of heat fluxes dominates over the fresh water and momentum forcing. Delworth and Dixon (2000) argued that the thermohaline circulation response to the intensifying Arctic and North Atlantic Oscillations can delay the weakening of the former due to increase of greenhouse gases concentrations. Timmermann et al. (1998), analyzing an experiment with the ECHAM3/LSG atmosphere-ocean climate model, stressed that anomalous evaporation and Ekman transport are essential for interdecadal variations in the thermohaline circulation. Groetzner et al. (1998) in the experiment with a coupled ECHO model, found oscillations in the North Atlantic with a period of 17 yr as the result of coupling between ocean gyre circulation and extratropical atmospheric circulation. However, the results of coupled models cannot be directly compared with the observations. Ocean GCMs, on the other hand, while unable to describe the coupled modes of the ocean-atmosphere interaction, can simulate the observed changes in the ocean, forced by the real air-sea fluxes over the last several decades.
These experiments were made possible when long time series of global air-sea fluxes became available from the voluntary observing ship (VOS) data (e.g., da Silva et al. 1994;Josey et al. 1999) and reanalyses (Kalnay et al. 1996;Gibson et al. 1997;Kistler et al. 2001). These products initiated a number of experiments with ocean GCMs driven by forcing functions using real (observed or analyzed) atmospheric variables for the last 20-40 yr (Battisti et al. 1995;Halliwell 1998;Hakkinen 1999;Beismann 1999;Ezer 1999;Seager et al. 2000;Eden and Willebrand 2001). The general concept of these experiments is to use forcing, varying year-to-year, after a model spinup that was driven by monthly climatological forcing fields. Battisti et al. (1995) and Halliwell (1998) showed that surface heat and moisture flux anomalies are the most important factors of the generation of the SST anomalies in the model. Halliwell (1998) in an experiment with 2Њϫ2Њ model found that anomalous air-sea flux patterns, associated with anomalies in atmospheric circulation, controlled by the NAO, are largely responsible for the interannual to decadal-scale SST anomalies generated in the regions of strong westerlies and trade winds. In the Gulf Stream region, however, changes in the heat advection by ocean flow contributed significantly to the generation of SST anomalies. Idealized NAO-like forcing was used by Krahmann et al. (2001) for the simulation of the propagation of SST anomalies, diagnosed by Sutton and Allen (1997). Seager et al. (2000) showing the dominant role of surface heat fluxes in driving SST anomalies, argued that ocean heat transport can be important in the subpolar gyre. Marshall et al. (2001) have shown that changes in the NAO forcing associated with meridional shifts of the zero wind stress curl line, can force an intergyre circulation responsible for the changes in the Gulf Stream location, and suggests a possible feedback mechanism of this intergyre exchanges on the NAO forcing. Shifts of the latitude of the Gulf Stream, lagged by ϳ2 yr with respect to the NAO index, were reported by Taylor and Stephens (1998). Joyce et al. (2000), however, found the highest correlation of the Gulf Stream position with the NAO index at zero time lag. While these studies enlighted our understanding of the response of the SST to the NAO forcing, and suggested possible feedback mechanisms (e.g., Marshall et al. 2001), they paid less attention to the changes of water mass properties associated with the variations in the overturning and MHT. Thus another objective of the present study is to identify the changes in the water mass properties linked to the variability in meridional circulation.
Some studies have, however, analyzed deep ocean circulation changes in response to atmospheric forcing. Hakkinen (1999) in an experiment with an ocean GCM driven by Comprehensive Ocean-Atmosphere Data Set (COADS) surface flux anomalies, found that surface heat forcing led MHT at 25ЊNb yϳ3 yr. Eden and Willebrand (2001), in a series of sensitivity experiments forced by National Centers for Environmental Prediction-National Center for Atmospheric Research (NCEP-NCAR) fluxes found both fast (barotropic) and delayed (baroclinic) responses of the ocean thermohaline circulation to the NAO. Visbeck et al. (1998), who used the NCEP-NCAR winds and heat fluxes derived from the atmospheric boundary layer formulation of Seager et al. (1995), showed, at least in the Tropics and subtropics, a fast response of the ocean circulation to NAO with little contribution from the ocean dynamics. Eden and Jung (2001) studied interdecadal variability in a 130-yr simulation, forced by reconstructed surface fluxes representing a NAO-like signal. They concluded that on interdecadal timescales the role of surface heat fluxes dominates over the wind stress role. This result is not, however, necessarily valid on interannual timescales. Although many general features of these experiments are consistent, less agreement is found in the description of the key mechanisms controlling ocean circulation. Since the variability patterns in different flux products are quite comparable (e.g., Gulev et al. 2003, manuscript submitted to J. Phys. Oceanogr., hereafter GUL), differences in the models performance, algorithms of the surface forcing and coordinate systems have to be responsible for the discrepancies observed. Willebrand et al. (2001) reported strong differences in the circulation patterns in the northern North Atlantic in several GCMs. The DYNAMO project (DYNAMO Group 1997) has identified considerable differences between the circulation characteristics in three ocean GCMs driven by the same forcing functions, but based on different coordinate systems. Thus, more simulations with alternative models are needed to quantify the mechanisms of variability in the North Atlantic thermohaline circulation. A question to address is the robustness of the description by ocean GCMs of the mechanisms generating the variability of the North Atlantic thermohaline circulation.
The present study further answers this question by presenting results from a topography-following coordinate primitive equation ocean GCM driven by 40 yr  of NCEP-NCAR reanalysis surface fluxes. We will try to represent the variability in the surface heat and freshwater forcing in terms of the density fluxes modified by the model (Nurser et al. 1999) and analyze the response of the ocean to the surface energy exchange in order to identify the mechanisms of interannual variability in the North Atlantic. Since the ocean ''feels'' the surface heat and moisture fluxes through the surface water mass transformation, this approach allows us to better link the air-sea exchange of heat and fresh water with the ocean processes such as variations in the intensity of meridional circulation. Section 2 describes the model experiment and the forcing. Analysis of the model surface density flux (effective forcing) is presented in section 3 in terms of transformation rates of mode waters and their variability. In section 4 we analyze the low-frequency changes in meridional overturning and MHT and attempt to identify the links between the variations in the meridional transports and in the water mass transformation rates. In section 5 we discuss our results in the context of the other experiments and observation and address the need for further studies.

a. North Atlantic model configuration
We used the sigma-coordinate primitive equation model (SPEM; Haidvogel et al. 1991). This model was used in an eddy-resolving version (1/3Њ) in the DY-NAMO project, which compared three different North Atlantic models (DYNAMO Group 1997;). The current version was adapted to a coarser resolution and modified on the basis of the DYNAMO results. Details of the model configuration and spinup are described in Knochel (1998) and Barnier et al. (2003, manuscript submitted to Ocean Modell., hereafter BAR). The computational domain covers the North Atlantic from 16.5ЊSto70ЊN (Fig. 1). The horizontal resolution is 1.2Њ cos in latitude () and 1.2Њ in longitude (). The vertical resolution corresponds to 20 nonequidistant levels with a higher discretization near the surface. The model has three buffer zones, at 70ЊN, 16.5ЊS and in the Strait of Gibraltar, where temperature and salinity are restored to climatological values. At 70ЊN we used the restoring data specially processed for the Dynamics of North Atlantic Models (DYNAMO) project (DYNAMO Group 1997). Levitus (1982) climatology is used for the southern boundary and the Strait of Gibraltar. Willebrand et al. (2001) concluded that the restoring allows the model to be prognostic from the Greenland and Norwegian Seas to the equatorial South Atlantic. Bottom topography was smoothed according to the criterion of Barnier et al. (1998). The subgrid-scale parameterizations follow the DYNAMO guidelines . Lateral harmonic mixing has been provided by the rotation of the mixing tensor from to geopotential coordinates (Barnier et al. 1998), with the mixing coefficients (0.8 ϫ 10 4 m 2 s Ϫ1 for momentum and 0.2 ϫ 10 4 m 2 s Ϫ1 for tracers at the equator) varying inversely with the grid size to keep the grid Reynolds number constant with latitude. Other subgrid-scale parameterizations such as mixing along neutral surfaces or eddy induced transport (such as the one proposed by Gent and McWilliams 1990) were not avail- able in the SPEM model at the time of the experiment. Therefore, diapycnal mixing is probably too large in the present experiment, with, however, very minor consequences on the relatively short variability timescales (3-7 yr) discussed here. The vertical mixing coefficient is 10 Ϫ3 m 2 s Ϫ1 for the momentum and depends on the local buoyancy frequency for the tracers (Gargett 1984). A mixed layer model, based on Kraus and Turner (1967), is used to provide the upper-layer mixing. Further details of the model performance, including parameterizations of the bottom friction, are described in BAR.

b. Forcing formulation
Different parameterizations of the surface forcing have been used in numerous model studies focusing on atmospherically driven variability of the ocean. These formulations exhibit a variety of approaches to account for the air-sea coupling or feedback and for possible biases in air-sea flux estimates (Barnier et al. 1995;Seager et al. 1995;Large et al. 1997;Eden and Willebrand 2001). Very few studies have carried out systematic-in-depth comparisons of these formulations [see Paiva and Chassignet (2001) for a recent effort in this direction], and more studies are required to provide objective arguments to prefer one parameterization over another, according to scientific objectives.
In the present study, the surface forcing in the model was applied according to Barnier (1998). The surface heat flux applied at a time t as the surface boundary condition of the temperature equation, Q net (t), is formulated as ‫ץ‬Q is the net surface heat flux from the NCEP-QЈ net NCAR reanalysis at this time (referred to as the original NCEP flux in this study), and ‫ץ‬Q/‫ץ‬T is a feedback term, expressed as the sensitivity (or derivative) of the net heat flux with respect to the sea surface temperature. We use the spatially varying climatological monthly mean values calculated by Barnier et al. (1995). SST(t) is the NCEP-NCAR SST at this time (referred to as the observed SST in this study), and T s (t) is the surface temperature calculated by the model, also at this time.
The salinity equation is forced by a virtual salt flux (Barnier 1998), proportional to a freshwater flux, (E Ϫ P)(t), which reflects the water budget at time t: where E(t) is evaporation and P(t) is precipitation, both taken from the NCEP-NCAR reanalysis at time t, R(t) is the river runoff, S(t) is the surface salinity calculated by the model at time t, is a relaxation time, and ⌬z is the thickness of the first model layer. SSS should be equivalent to an observed sea surface salinity at the same time t. However, the sea surface salinity is poorly observed at weekly or monthly timescales in comparison to the SST. Thus, we used in (2) the climatic seasonal estimates from the climatology of Levitus (1982) interpolated at time t. Since these data do not include interannual variability, the relaxation time scale is set to 365 days, to limit the contribution of the relaxation term in (2). Consequently, interannual variability in the model freshwater forcing is exclusively influenced by the NCEP-NCAR E and P fields. The river runoff R(t) is also difficult to parameterize, since there is no reliable estimates of its variability over the period 1958-97. Therefore, the value of R(t) in (2) has been parameterized by a strong relaxation (time constant of 3 days) of the model surface salinity to the Levitus (1982) climatology in the regions where the SSS is less than 32 psu. This hypothesis assumes that local salinity in these regions is primarily set by the freshwater input from rivers or ice melt. In our model, these regions are the mouth of the Amazon, Congo, Mississippi, and St. Lawrence Rivers, and the shelf regions along the coasts of East Greenland and Canada.
The formulation of the net surface heat flux and freshwater forcing used in (1) and (2) is inspired by the development of surface boundary conditions suggested by Haney (1971), Takano et al. (1973), and Han (1984). However, compared to the simple relaxation to the SST and surface salinity, this formulation accounts for some feedback from the ocean to the atmosphere and allows us to explicitly apply the air-sea fluxes. Thus, it is different from the bulk formulation used by Danabasoglu and McWilliams (1995), or Large et al. (1997), or from the formulation of Seager et al. (1995), which couples a simplified atmospheric advective-diffusive boundary layer to the ocean model to calculate the fluxes. Like these other parameterizations, our formulation allows for prognostic meridional heat and freshwater transports. As a consequence of the formulation of the feedback term used here, the SST and the SSS in our experiments are identified more with forcing functions than with prognostic variables.
The momentum equation is forced by the NCEP-NCAR wind stress as the surface boundary condition. The sea-ice parameterization sets the surface flux to zero and implies a 3-day relaxation of the model SST to the freezing temperature and of the model SSS to the climatic SSS, when the NCEP-NCAR surface temperature falls below Ϫ1.8ЊC.

c. Interannual forcing data: NCEP-NCAR reanalysis
Forcing fields for the model simulation were obtained from NCEP-NCAR reanalysis for the period 1958-97. Net heat surface flux was derived from 6-hourly sensible and latent heat, net short and long wave radiation. Evaporation minus precipitation fields were derived from the latent heat flux and precipitation. Components of the wind stress were used for the momentum forcing and for the computation of the cube of the friction velocity, , which is required by the Kraus and Turner (1967) 3 u * mixed layer model. Values of were computed from 3 u * 6-hourly stress components to account for the strong nonlinearity involved. The NCEP-NCAR SST combines the Reynolds and Smith (1994) SST analyses for the period starting from 1982 and the Global Sea Ice and SST dataset (GISST; Rayner et al. 1996) for the earlier years. Climatological monthly SSS were taken from Levitus (1982). All NCEP-NCAR 6-hourly fields were interpolated from the Gaussian grid (1.875Њϫ ഠ1.92Њ) to the model grid using the modified method of local procedures (Akima 1970). The model was forced by monthly values of forcing components, which were linearly interpolated to the model time steps.
NCEP-NCAR surface fluxes, especially precipitation, include some known biases. Thus, ad hoc corrections were applied to the precipitation fields, particularly for some points in the offshore regions east of Greenland and in the southwestern Caribbean Sea, where NCEP-NCAR shows precipitation anomalies several orders of magnitude higher than in the neighboring locations.
These grid points were excluded from the data and the values were reinterpolated from the neighboring points using the Akima (1970) method. Sensitivity experiments showed that artifacts in precipitation east of Greenland seriously change the properties of waters advected into the Labrador Sea and can affect the convection. In order to evaluate monthly mean SST in the region of possible presence of sea ice, we used NCEP-NCAR surface skin temperatures at 6-hourly intervals to estimate for every model grid point the durations of continuous ''ice'' and ''open water'' periods. Then the total continuity of the durations longer than 72 h was estimated. In the case when this continuity was longer than 50% of the month, the grid point was set to open water, otherwise to sea ice, and the corresponding values were used as the monthly means. In this way we achieved reasonable areas of the ice cover in the northern regions and depicted realistically the conditions leading to convection (or its blocking) in the Labrador Sea. Direct estimation of monthly mean SST from 6-hourly surface skin temperature may result in the drop of SST below Ϫ1.8ЊC even in the case of only 3-4 days of the ice presence in the grid point. Alternatively, the use of only SST values associated with the open water for monthly averaging leads to unrealistic increase of the open water area in the Labrador Sea.

d. Spinup and interannual experiment
Before the 40-yr interannual experiment, the model was spun up for 20 yr. The spinup was forced by the 40-yr climatological monthly NCEP-NCAR fluxes. Thus, the mean forcing is the same for the spinup and variability experiment. A detailed analysis of the spinup experiment is given in BAR. The mean state diagnoses the major features of the North Atlantic circulation (Figs. 1a,b). The maximum overturning streamfunction is slightly over 10 Sv (1 Sv ϵ 10 6 m 3 s Ϫ1 ) and occurs at 43ЊN at 1000-m depth. The advective MHT reaches the maximum of 0.72 PW at 15ЊN. It decreases to 0.3 PW in the latitudinal band 40Њ-43ЊN, where the Gulf Stream runs zonally and the MHT is done by the eddy diffusivity (Fig. 2). Then MHT increases to 0.44 PW at 54ЊN and decreases to zero at the northern boundary. This is typical behavior of these models run with similar parameterizations for diffusivities (Marsh et al. 1996). More recent parameterizations of mixing along neutral surfaces tend to produce greater overturning (Böning and Bryan 1996).
The present model solution represents the location of the Gulf Stream well after its separation at Cape Hatteras (Fig. 1a). It also shows a cyclonic circulation cell of cold water between the north wall of the Gulf Stream and Newfoundland, which is a remarkable advantage of the present SPEM model simulation. This feature, which seems to be a common feature of sigma-coordinate models (Dengg et al. 1996), has been related to an improved representation of the deep western boundary current around the Flemish Cap (BAR). Temperature and salinity drifted away from the initial conditions during the spinup, getting warmer and saltier, which is also common for this type of simulation. The drift is below 0.002ЊCy r Ϫ1 in temperature, and 0.0005 psu yr Ϫ1 in salinity. BAR analyzed the water masses produced by the model at the end of the spinup and concluded that the major mode waters (subtropical, subpolar, and Labrador Sea Mode Waters) are represented. We describe below the main characteristics of mode waters in the interannual experiment.
STMW are found in the western subtropical gyre to the south of the Gulf Stream front, between 25Њ and 35ЊN ( Fig. 3a). They are characterized by a thermostat (i.e., an area of minimum vertical temperature gradients) extending from 50 and 300 m in the temperature range of 19Њ-21ЊC. In the observations (Levitus 1982) the STMW thermostat is found slightly deeper and colder (17Њ-20ЊC). The warm to cold conversion of subtropical waters occurs along the Gulf Stream and the North Atlantic Current (NAC). Waters are progressively cooled as they flow northeastward. The Subpolar Mode Waters (SPMW) are the product of this conversion. They are found in the northeast Atlantic between 40Њ and 55ЊN, and are characterized by a thermostat between 50 to 400 m in the temperature range of 13Њ-11ЊC (Fig. 3b). The characteristics of the model SPMW compare well with those described by McCartney and Talley (1982). The third type of mode waters represented in the model is the LSW which is found at a temperature range of 4Њ to 4.5ЊC. Mean properties in T, S of the model mode waters are summarized in Table 1.
For the diagnosis of the mixed layer depth (MLD) we required that 0 (surface) Ϫ 0 (MLD) ϭ 0.01. The MLD according to this criterion has been exactly computed by a vertical interpolation, and not approximated to a discrete model level. March values of the MLD in the Labrador Sea, slightly over 1200 m, are realistic. During the interannual experiment the deepest mixed layer in the Labrador Sea was reached in 1973 (2000 m) and in 1993 (1800 m). Climatological values of the MLD in the Labrador Sea and its interannual variability simulated in our model agree well with the observations of Clarke and Gascard (1983), Clarke (1983), Lasier (1980), and Curry et al. (1998), showing the major extremes of MLD, such as negative anomalies in the early 1980s and high values in the mid 1970s and 1993-94. Eden and Willebrand (2001) in a model experiment forced by the surface heat flux only got approximately the same tendencies in the MLD anomalies in the Labrador Sea as in our simulation. After the spinup phase, the experiment has been continued, being forced from 1958 to 1997 by the NCEP-NCAR fluxes applied as specified in section 2b. The drift of the model properties (T, S) observed during the 20-yr spinup is nearly zero during the 40 years of the interannual experiment.

a. Surface density flux
The actual heat and freshwater fluxes applied to the model significantly differ from the original NCEP-NCAR flux estimates, as the forcing is formulated as (1), (2). These differences result from the correction (or feedback) terms in (1), (2). Since the model is forced by the changes in buoyancy at the surface, we have calculated the density flux, f , from the net heat and freshwater fluxes diagnosed by the model, using a formula proposed in a number of works (see Large and Nurser 2001, e.g.). According to these works, the density flux (in fact a virtual mass flux since it has the unit of kg m Ϫ2 s Ϫ1 ) at the sea surface has been derived as where C P is specific heat of seawater at constant pressure, 0 is a reference density of seawater, s is salinity in portions of unity, and ␣ and ␤ are the thermal expansion and haline contraction coefficients: Analytical expressions for ␣ and ␤ are obtained by differentiation of the equation of state used in the SPEM model (Jackett and McDougall 1995) with respect to T and S, and surface values are obtained by using surface temperature and surface salinity in these analytical expressions. Figure 4a shows the 40-yr mean surface density flux derived from the NCEP-NCAR data (the original fluxes hereafter). Figure 4b shows the correction of this flux diagnosed by the model. The algebraic sum of these is the actual density flux applied to the model.

1) CLIMATOLOGICAL 40-YR MEAN
The original density flux ( Fig. 4a) indicates large positive values (generation of denser water) over the Gulf Stream (where the maximum is 12 ϫ 10 Ϫ6 kg m Ϫ2 s Ϫ1 ), in the North Equatorial Current (the southern limit of the subtropical gyre), and in the Labrador Sea where the convection region at 60ЊN is well identified. The most intense negative density flux (generation of lighter water) is found in the equatorial Atlantic (Ϫ5toϪ6 ϫ 10 Ϫ6 kg m Ϫ2 s Ϫ1 ), off Newfoundland and in coastal regions near Portugal and Senegal. Figure 4a reflects major geographical features of the distribution of the heat and freshwater fluxes and agrees with that obtained by Garnier et al. (2000) from the European Centre for Medium-Range Weather Forecasts (ECMWF) reanalysis.
The flux correction term (Fig. 4b) allows the model MHT to be prognostic. Locally large corrections (with either sign) indicate that inconsistencies exist between the model dynamics and the forcing data (SST, the original NECP fluxes, and to a lesser degree the climatological SSS). Small correction values indicate that forcing fields are consistent with the model dynamics. The spatial distribution of the correction term (Fig. 4b) is qualitatively very similar to the original fluxes (Fig. 4a). It generally reduces the strength of the original flux (signs are opposite) almost everywhere except in the Labrador Sea, where the correction increases the density loss, and along the equator. It can be hypothesized that underestimation of surface velocities (first of all in the western boundary currents) in a coarse resolution model leads to a longer time for the surface water to be subjected to the real surface flux than in nature implying a negative flux correction. In the western equatorial Atlantic (off Brazil), the correction amplifies the buoyancy increase of the surface water. The pattern is, however, peculiar since the maximum buoyancy gain occurs in the central Atlantic in the original flux and in the western Atlantic in the corrected flux. Thus, for these areas the correction changes the magnitude of fluxes but does not displace the forcing centers of action. The reason for that is the correct location of the main surface currents in the model. The correction is smaller than the original flux, contributing from 5% to 50% to the resulting buoyancy forcing. Along the Gulf Stream and in the North Equatorial Current the correction accounts for 50% to 40% of the total forcing. In regions of deep convection the correction is quite small; that is, the model dynamics seems to be consistent with the original NCEP-NCAR fluxes. The strong relaxation applied in the regions of low salinity (S Ͻ 32 psu), which provides a parameterization of river runoff and sea-ice melt, generates light waters along the coast of Greenland and in the western Labrador Sea. The spatial pattern of our flux correction term is qualitatively comparable to that of Nurser et al. (1999). Our correction is somewhat stronger in the Gulf Stream area and is considerably smaller at low latitudes and in the Labrador Sea.

2) INTERANNUAL VARIABILITY
The interannual variability patterns for the original NCEP fluxes and the corrected fluxes are quite comparable in the midlatitude and subtropics. The first three empirical orthogonal functions (EOFs) of the original NCEP-NCAR and total model fluxes are nearly identical and the corresponding principal components (PCs) are highly correlated. Correlation coefficients between the wintertime flux anomalies of the original and the corrected flux (not shown) are higher than 0.9 over most areas north of 20ЊN. Somewhat smaller correlation (0.75-0.9) is observed over the Gulf Stream region. Equatorial areas are characterized by the drop of correlation to 0.3-0.5. In general this is in agreement with Eden and Willebrand (2001), who, reported however, slightly negative correlations in the Equatorial area. If the linear trends are removed, correlation increases by approximately 10%. Consequently, the feedback term has not changed the nature of the variability of the forcing, which remains the one contained in the original NCEP-NCAR fluxes. Analysis of correlation between the model SST and the observed SST ( Fig. 4c) also shows very high correlation, ranging from 0.7 to 0.98 and increasing by 8%-13% after the trends were removed. This is higher than in the experiment of Hakkinen (1999) and is in agreement with our forcing formulation, where the model SST, being a part of the forcing function, is strongly constrained.
The variability of the surface density flux is described by its first and second EOFs, which explain 36% and 15% of the total variance, respectively (Figs. 5a,b). The first EOF identifies subpolar, subtropical and tropical centers of action. The second EOF exhibits a dipole structure with centers of action in the Florida Current and the Gulf of Mexico, and in the region south of Newfoundland along the path of the North Atlantic Current. A similar mode of variability has been identified in a number of studies of surface heat fluxes (Cayan 1992a,b;GUL). The first normalized PC of the total density flux shows a correlation of 0.73 with the NAO index derived from the NCEP-NCAR sea level pressure (SLP) data (Fig. 6). GUL reported higher correlation (0.85) between the first normalized PC of the NCEP-NCAR net heat flux and the NAO index. Visible disagreements between the two time series occur between 1962 and 1966 and during the 1990s. We have investigated separately the contribution of the surface thermal and haline density fluxes to the variability of the total density flux (no figure shown). The first EOF of the thermal component accounts for 37%, and that of the haline component for 19% of the total variance, respectively. The magnitudes of the leading mode for the thermal density flux are higher than for the haline flux by a factor of 5 to 10 in the mid-and subpolar latitudes and subtropics and by a factor of 1.5-3 in the equatorial regions. The first EOF of the thermal flux largely resembles the corresponding pattern of the total density flux (shown in Fig. 5a) with 10% to 20% smaller magnitudes. The first EOF of the haline density flux exhibits an east-west dipole in the midlatitudinal North Atlantic. The first normalized PCs of the thermal and haline density fluxes (not shown) are highly cor-related with each other (correlation of 0.85) and with the first PC of the total density flux (correlation with the first PC of the thermal flux is 0.97, and with the first PC of the haline flux is 0.81).

b. Transformation of water masses
Now we investigate the changes in the water masses, produced by the density forcing. Following Walin (1982), Speer and Tziperman (1992), Marsh (2000), and others, we integrate in space and time the density flux (3) to obtain the transformation rate of waters at given density, F( ): ͵͵ ͵ ⌺ where the ␦ function samples the density flux f over the area where waters of density are outcropping within the integration area ⌺. This quantity represents the time average over a period (generally 1 yr) of the density change of a unit water volume of density due to the action of the atmospheric forcing. It has a unit of kg s Ϫ1 and can be scaled with density (kg m Ϫ3 ) (Nurser et al. 1999;Speer et al. 1995) to be expressed in m 3 s Ϫ1 or Sv. In this form the transformation rate shows the volume of water of density , which is transformed during a given period into higher densities [F( ) Ͼ 0] or lower densities [F( ) Ͻ 0]. Figure 7a compares the climatological mean transformation rate derived from the original NCEP-NCAR surface fluxes with that diagnosed by our experiment (the difference between the curves comes from the flux correction). The reference climatological winter map of surface density (Fig. 7b) helps to locate geographically the density scale in Fig. 7a. The curves in Fig. 7a are qualitatively comparable with those of Speer and Tziperman (1992), who used the Isemer and Hasse (1987) flux climatology and Levitus (1982) surface temperature and salinity, and found slightly larger rates than those obtained from the NCEP-NCAR fluxes. Relationships between the transformation rate derived from the original surface fluxes and the fluxes diagnosed by the model are somewhat smaller than those reported by Nurser et al. (1999). In agreement with the pattern of the flux correction described before (Fig. 4b), the model transformation rates are weaker in the tropics and subtropics (density range from 24 to 27, the flux correction reduces the net flux), and are stronger at the highest subpolar latitudes indicating more intensive transformation of the densest waters ( 0 ϭ 27.7).
To better identify the waters whose transformation can be associated with significant ocean-atmosphere exchanges, we have calculated the transformation rate only for the area where the winter MLD is greater than 100 m (shaded area in Fig. 7b). This allows us to consider only the regions where the impact of the atmospheric buoyancy forcing reaches a significant depth. The corresponding curve (Fig. 7a) shows two transformation peaks. Part of a broad peak (25.5 Ͻ 0 Ͻ 26.5) at densities lower than 26.5 corresponds to the formation of the STMW (18ЊCat 0 ഠ 26-26.5). They are formed by winter densification of waters south of the 0 ഠ 26 outcrop density line inside the western subtropical gyre to the south of the Gulf Stream front. At densities greater than 26.5, this peak corresponds to the formation of the SPMW (10Њ-14ЊCa t 0 ഠ 26.8 to 27.1) by the densification of waters with 26.5 Ͻ 0 Ͻ 27.0 located in the path of the NAC and the eastern North Atlantic between 40Њ and 55ЊN. The narrow peak centered at 0 ഠ 27.7 (densities that are only found in the Labrador Sea and Irminger Sea) corresponds to the formation of the dense LSW (4ЊCa t 0 ഠ 27.8). Figure 8 shows anomalies of the winter transformation of surface water for the surface density classes higher than 0 ഠ 25.5 over the area where the winter MLD is deeper than 100 m, from the early 1960s to the mid 1990s, smoothed with a 3-yr running mean. Pronounced decadal scale variability is observed in the transformation of very dense waters ( 0 ഠ 27.75-27.85) and of waters in the density classes corresponding to the STMW (26.25 Ͻ 0 Ͻ 26.5). Periods of primarily outof-phase behavior of the transformation of these waters, pronounced in the late 1960s-70s, and in the 1990s are superimposed with periods of coordinated changes in the 1980s. The overall correlation coefficient between the winter averages of water mass transformation in the classes mentioned above (Ϫ0.52) is significant at the 99% level.
The next step is to consider the behavior of the transformation rate F in the space of surface temperature and salinity, because the T-S properties of mode waters are clearly defined from observations. Figure 9a shows climatological winter surface water mass transformation in the T-S plane (instead of density), and mode waters are identified according to Table 1. It shows fairly continuous, positive transformation rates in the temperature range of 24Њ-12ЊC, associated with the formation of the STMW (17Њ-20ЊC), and with the densification of the central North Atlantic waters flowing to the northeast in the NAC to give birth to the SPMW (10Њ-14ЊC). Large positive transformation (density flux to the ocean) is also associated with the LSW (4ЊC). The leading EOF (36% of the total variance) of the winter transformation rate derived from the model experiment is shown in Fig.  9b. It presents a pattern with the highest anomalies of opposite signs associated with the LSW and STMW, and a significantly weaker signal for SPMW. Note that the transformation of STMW is associated with weak transformation of the opposite sign at the same densities, but slightly smaller temperatures and salinities. This local extremum, however, is characterized by very small locally explained variance. The corresponding normalized PC (Fig. 10a) exhibits a pronounced decadal variability with the strongest surface diapycnal density flux in the subtropics in the late 1960s and strong positive anomalies, corresponding to the T-S classes of LSW in the mid 1970s and 1990s with the extremes in 1969 and 1989.
The first normalized PC of the water mass transformation exhibits a very close relationship with the NAO index (Fig. 10a). The correlation between the two time series is fairly high (0.79). Figure 9c shows the associated correlation between the NAO index and the water mass transformation anomalies in the T-S plane. An extreme positive correlation higher than 0.6 is identified over the T-S classes corresponding to the transformation of LSW and over the other edge of the T-S diagram corresponding to the warm (T Ͼ 20ЊC) and salty (S Ͼ 37 psu) waters of the eastern tropics. A significantly negative correlation of Ϫ0.5 is observed for the T-S classes associated with water mass transformation in the western Atlantic subtropics (STMW). We therefore expect, on the one hand, a high correlation between the NAO and changes in formation rates of STMW (out of phase) and LSW (in phase, and in agreement with the positive correlation found between the NAO and the MLD in the Labrador Sea noticed in Fig. 10b). On the other hand, SPMW, which shows a weak negative correlation, is expected to respond differently to changes in the atmospheric forcing.

c. Linking patterns of variability of the buoyancy flux and transformation rates
In order to identify the optimally correlated patterns of the transformation anomalies in the T-S plane and surface density fluxes, we applied a canonical correlation analysis (von Storch and Zwiers 1999) to the anomalies in surface density fluxes and transformation rates. The first five EOFs, which account altogether for 83% of the total variance of surface density and 90% of the total variance of the water mass transformation, were taken for the computations. The first canonical pair (Figs. 11a,b) (correlation coefficient is 0.92) largely resembles the first EOFs of the surface density flux (shown in Fig. 5a) and of the water mass transformation (shown in Fig. 9b), and exhibits negative anomalies of the transformation of STMW associated with the enhanced production of LSW. The second canonical pair (Figs. 11c,d) (correlation is 0.81) is represented by the second EOF of the surface density flux (shown in Fig. 5b) and is optimally correlated with the transformation of the SPMW and the North Atlantic Central Waters (NACW) by the anomalous surface density flux over the central midlatitudinal Atlantic.

d. Wind forcing
An essential part of the forcing function is represented by the wind stress curl, which is largely driving the barotropic changes in the ocean circulation, in particular contributing to the changes in the Gulf Stream location (Marshall et al. 2001). The leading EOF of the wind stress curl (not shown) accounts for 29% of variance and shows the intensification of the large-scale cyclonic vorticity in subpolar latitudes and of anticyclonic vorticity in the midlatitudes and subtropics in agreement with Eden and Willebrand (2001). The corresponding normalized PC (Fig. 10a), representing the NAO signal (correlation with the NAO index is 0.89), shows coordinated variability with the leading PC of the water mass transformation (correlation is 0.80).
Thus, the most important signals in the surface forcing and water mass transformation in the North Atlantic clearly demonstrate NAO signatures. NAO-like surface thermal forcing results in the out-of-phase anomalies of the water mass transformation in the subtropics (STMW) and subpolar latitudes (LSW). The transformation of SPMW, which is well-correlated to the second mode of variability of the atmospheric buoyancy forcing, appears not to be linked to the NAO. Enhancement of the midlatitude atmospheric flow is associated with the dipole pattern in the wind stress curl.
The NAO signature is also clearly pronounced in the MLD diagnosed by the model, representing the immediate reaction of the ocean to surface fluxes. The first EOF (41% of the total variance) of the March MLD (not shown) is presented by the main center of action in the Labrador Sea in agreement with the results of Eden and Willebrand (2001). The corresponding normalized PC (Fig. 10b) shows clear decadal scale variability and a correlation with the NAO index of 0.72. Note that the correlation with the NAO during the period 1958-77 is somewhat higher than during the last two decades (0.77 and 0.66, respectively).

Variability in the meridional ocean circulation associated with the surface forcing a. Temporal variability of MHT
We now analyze the variability in the North Atlantic, associated with the surface forcing, of the MHT and the intensity of the meridional overturning cell. Figure 2a shows the climatological (1958-97) estimates of MHT from the equator to 66ЊN, implied by Q net (T s ) (thus including the correction) and by (the original NCEP-QЈ net NCAR flux) according to the formula ͵͵ net 0 0 where is latitude, is longitude, MHT 0 is a boundary condition at the north (66ЊN), which is set to the advective MHT diagnosed by the model. As we mentioned above, flux correction considerably modifies the original flux, imposing primarily positive correction in mid-and high latitudes and negative correction in the subtropics and Tropics (with differences over 0.2 PW in the 10Њ-20ЊN latitude band). The model advective MHT can be decomposed as (e.g., Kamenkovich et al. 2000) MHT ϭ Cd l ͗T͘ dz where ͗͘ is the zonal averaging operator and prime stands for the variations of temperature T and meridional component of velocity around the zonal mean, L is longitude, z is depth. The first term in the rhs of (7) is an overturning component, which corresponds to the transport of the zonally averaged circulation, and the second term in the rhs of (7) represents a gyre component, which corresponds to the transport by horizontal gyres.
The results of the MHT decomposition is shown back in Fig. 2. The gyre MHT has the usual shape and amplitude shown by other model solutions at similar resolution (Böning and Bryan 1996). In the Tropics the gyre component opposes the overturning component, but the latter dominates and reaches a maximum of 0.8 PW at 14ЊN. North of 20ЊN, both components provide the northward transport of heat, the gyre component reaching its maximum at 27ЊN. This maximum is usually reached at a higher latitude (33ЊN) in other models (Böning and Bryan 1996). A reasonable explanation for that is the possible northward shift of the center of the gyre in these models caused by the overshoot of the Gulf Stream. The overturning usually dominates, except at high latitudes (55ЊN and northward) where the gyre component carries most of the heat. The gyre component is insignificant in the 30Њ-40ЊN band, where currents are zonal and a considerable part of the MHT is provided by the diffusive component (Fig. 2a). The latter has been defined from the model potential temperature gradients and horizontal diffusivity. Ideally, the total advective MHT diagnosed by the model together with the diffusive component should converge to that derived from the actual surface fluxes applied to the model. Figure 2a shows that these two estimates are balanced on a large scale. Local imbalances are identified for different latitudes, which partly represent the model heat storage, but also the contribution of the buffer zones, which are not accounted for in this calculation, and thus appear in the residual. Thus, between 30Њ-40ЊN, the treatment of the Mediterranean outflow is certainly responsible for the difference between the curves. The impact of the northern buffer zone on the diffusive transport is drastic beyond 60ЊN. Figure 12 shows anomalies of different MHT components from the early 1960s to the late 1990s smoothed with a 3-yr running mean. The advective MHT (Fig.  12a) exhibits pronounced decadal variability, superimposed on the secular changes. There is a tendency of increasing MHT in the mid-and subpolar latitudes (30Њ-60ЊN) and weakening MHT in the Tropics (10ЊS-20ЊN). High positive anomalies in the Tropics in the 1970s correspond to negative anomalies north of 40ЊN. Starting in the late 1980s, strong negative anomalies in the Tropics are associated with primarily positive anomalies in the midlatitudes and subtropics (20Њ-40ЊN). The variability of the overturning MHT (Fig. 12b) is very similar to that for the total advective MHT from the equator to the mid latitudes. The gyre component (Fig. 12c) contributes mainly to the variability of the total advective MHT north of 40ЊN, where it increases from the 1960s to the 1990s, and is responsible for the strong positive anomaly of the total MHT at 45Њ-50ЊNi nt h e mid 1990s. In the equatorial band and in the midlatitudes gyre and overturning components demonstrate out-ofphase behavior. Anomalies of the diffusive component (not shown) are much smaller than for the other components. They are correlated with the anomalies of the gyre component between 35Њ and 50ЊN, contributing to the secular increase of the midlatitude MHT.
MHT in the low latitudes is largely influenced by the Ekman MHT (EMHT), primarily associated with the transport in the upper ocean. We derived the EMHT from the model temperature of the surface layer e and zonal wind stress x , as (e.g., Levitus 1989) where L is longitude, z is the vertically averaged potential temperature, is the density of water in the surface layer, and C p is the specific heat at constant pressure. Recently, Klinger and Marotzke (2000) and Jayne and Marotzke (2000) have accounted for the isopycnal return flow to achieve more accurate estimates of the Ekman transport. However, this primarily improves the time-mean estimates, while for the variations of EMHT, on which we are focused, (8) gives acceptable accuracy (Jayne and Marotzke 2000). The interannual variability of the EMHT (Fig. 12d) exhibits strong positive trends in the low latitudes and pronounced decadal variability in the subtropics. The anomalies of the maximum of the overturning stream function (no figure shown) are in clear agreement with the anomalies of the total advective MHT, as in Willebrand et al. (2001). Disagreements are observed only during the 1960s and mid-1990s in the Tropics, when the contribution of the gyre component to the MHT variability is significant. We computed the correlation at all latitudes between the anomalies of MHT and the maximum overturning stream function computed from the original and detrended time series. For the original time series the correlation is fairly high from 10Њ to 43ЊN and decreases sharply in the subpolar latitudes. After the detrending, the correlation between 45Њ and 60ЊN grows considerably, indicating that secular changes influence the relationships between overturning and MHT. In order to investigate whether an important part of the MHT variability originates from changes in temperature, we produced a diagram showing time evolution of the zonally averaged temperature anomalies in the upper 200-m ocean layer (Fig. 13). The most important heat storage occurred in the 50Њ-60ЊN latitude band during the period 1989-94, and with a smaller amplitude between 1972 and 1976, in relation with maximum deep convection in the Labrador Sea (Fig. 10b). The period 1976-92 was characterized by heat storage in the subtropical and the subpolar gyre (10Њ-50ЊN latitude band). No significant correlation was found between anomalies of the upper-layer heat content and anomalies of the model heat transport divergence, except for the 40Њ-50ЊN latitude band. At these latitudes, where the MHT is primarily associated with the NAC, we can speculate about a possible link between the advection of temperature anomalies and the heat storage.

b. Linking MHT changes with surface forcing characteristics
MHT responds to the surface thermal forcing, with some delay needed for baroclinic adjustment. Barotropic responses, on the other hand, may be quite fast. In order to quantify the dominant timescales at which the various components of the MHT respond to the variability in the forcing fields, we computed cross-correlation functions for relevant quantities at every latitude for time lags from 0 to 10 yr. Figure 14a shows the lagged correlation of the first PC of the surface water mass transformation (5) (leading) with the overturning MHT. The immediate reaction (lag 0) is quite weak (Ϫ0.30 to Ϫ0.34 between 45Њ and 54ЊN). At a 3-yr lag the overturning MHT shows a significantly positive correlation of 0.54 at 46Њ-53ЊN and a significantly negative correlation of Ϫ0.6 in the latitudinal band 20Њ-33ЊN. Significantly positive correlations of 0.5 to 0.58 are also observed at lag 7 between 40Њ and 50ЊN. The negative correlation at lag 7 in the Tropics is less pronounced.
Cross correlations of the gyre MHT (Fig. 14b) show that in the midlatitudes (20Њ-60ЊN), this component is positively correlated (0.72) with the leading PC of the surface water mass transformation at the 3-yr lag, and shows a negative correlation of Ϫ0.54 at lag 8 yr. The immediate reaction (lag 0) is characterized by a weak negative correlation of about Ϫ0.3 in the midlatitudes and south of 20ЊN. The total MHT shows a strong correlation with the leading PC of the wind stress curl at lag 0 (Fig. 14c). The strongest negative correlation is observed at 52ЊN( Ϫ0.71) and the highest positive correlation of 0.6 is identified at 32ЊN with a clear change of sign at 42ЊN.
Most of the responses of the MHT to the surface forcing are, thus, associated with changes in the intensity of the meridional overturning. Indeed, the lagged correlation for the maximum overturning (no figure shown) is basically identical to the results derived for the overturning MHT (Fig. 14a). Thus, the changes in MHT can be dynamically explained in terms of processes controlling the overturning. Note, that relationships for the lags of 7-8 yr can be influenced by the autocorrelation in both MHT components and forcing functions. We investigated the autocorrelation of the total and gyre MHT (no figure shown), and found significantly negative autocorrelation at lags 7-8. Conversely, the 3-yr lag is clearly not influenced by any significant autocorrelation.
In order to depict the dominant patterns of variability in the MHT, and to relate them to the variability patterns of the surface forcing, we decomposed MHT anomalies into EOFs. Figure 15a shows the first two EOFs of the total MHT, which account for 55% and 21% of the total variance respectively. They are well separated from each other and from the higher order EOFs. Corresponding normalized PCs are shown in Figs. 16a and 16b. The first EOF resembles the overturning contribution to MHT (Fig. 2) with the sign changed. Its PC (Fig. 16a) is lagged by a few years with respect to the first mode of the surface water mass transformation and the NAO index (Fig. 10). The second EOF shows a maximum in the latitudinal band 20Њ-40ЊN and demonstrates pronounced decadal-scale variability (Fig. 16b), reaching a maxima in the late 1960s, early 1980s, and early 1990s. This pattern and its variability resemble, in many respects, the MHT changes simulated by Hakkinen (1999). Figures 15b and 15c show the first two EOFs (46% and 27% of the total variance, respectively) of the overturning streamfunction. The first mode describes variations in downwelling at 55Њ-60ЊN (the region of formation of LSW) and in the strength of the northward flow north of 42ЊN, and is clearly compatible with the first MHT mode of Fig. 15a. The pattern of the second mode shows an intensified downwelling between 45Њ-55ЊN (the region of formation of SPMW), and increased upwelling between 15Њ and 35ЊN, particularly intense near the surface at 25ЊN (the region of formation of Fairly good correlation between the corresponding PCs of MHT and meridional overturning streamfunction (0.87 and 0.56) indicates that meridional overturning is largely responsible for the first two leading modes of the MHT. However, during the 1990s second modes demonstrate a lesser agreement, primarily due to the trends of opposite sign in the two time series. When the trends are removed, correlation between the second PCs increases up to 0.63. In order to identify the patterns of variability in surface forcing and meridional circulation which are associated with each other, we applied canonical correlation analysis to the surface forcing fields and MHT for time lags from 0 to 10 yr. Table 2 shows first three canonical correlation coefficients for different time lags between surface water mass transformation F and MHT and between wind stress curl and MHT. For the wind stress curl the highest correlation for the first pattern is at lag 0. Surface water mass transformation shows high correlation for the first two patterns with maxima at lags 3 (0.91 and 0.74) and 7 (0.87 and 0.76) and for the first canonical pattern at lag 0 (0.84), which are consistent with the cross-correlation analysis (Fig. 13). The canonical patterns for the corresponding parameters for the time lags of 0, 3, and 7 yr are shown in Figs mally correlated at lag 0 with the MHT pattern displayed in Fig. 17c. The corresponding leading canonical pattern of the barotropic streamfunction (Fig. 17b) is associated with high (low) NAO phases and describes a simultaneous weakening (strengthening) of the subtropical and subpolar gyres. This canonical pattern is slightly disrupted by a variation of the flow around the Azores plateau, which varies in phase with the subpolar gyre, being a part of the subtropical gyre. Thus, the MHT curve in Fig. 17c, which represents the first EOF of MHT with a reversed sign, corresponds to wind-driven changes in the barotropic streamfunction associated with the deintensification (intensification) of both subpolar and subtropical gyres. According to Fig. 17c, this imposes a decrease (or increase) of MHT north of 42ЊN. When NAO index is high, a strong negative anomaly in wind stress curl north of 40ЊN reduces the strength of the western part of the subpolar gyre (the first EOF of the barotropic flow indicates a weaker Gulf Stream and NAC), which imposes a decrease of MHT north of 42ЊN. This is in agreement with the results of Eden and Willebrand (2001), who found primarily negative correlation between MHT at 48ЊN and the NAO index in their experiments forced with wind anomalies. The strong anticyclonic intergyre gyre (emphasized by reddashed contours in Fig. 17b), induced by a positive NAO wind anomaly, is indeed very similar to the schematic proposed by Marshall et al. (2001) in their Fig. 6. It is interesting to note that this circulation pattern mostly brings subpolar waters across the boundary of the gyres (the blue line in Fig. 17b), and therefore is consistent with a cooling of the SST in the subtropics as suggested by the scheme of Marshall et al. (2001). However, since the SST is not prognostic in our experiment, we cannot demonstrate further a warming of the northern latitudes and a cooling of the southern latitudes, which would feed back on the position of the jet stream, as suggested by Marshall et al. (2001). Nevertheless, our model unambiguously confirms the pattern of the wind-driven ocean response suggested by the idealized study of Marshall et al. (2001).

2) WATER MASS TRANSFORMATION AND MHT AT 3-YR LAG
Figures 18a-c show the first two canonical patterns of the surface water mass transformation and MHT at lag 3 (water mass transformation is leading). The first canonical pair shows the highest correlation among all lags studied (0.91 in Table 2). It is represented by the first EOF of the winter surface water mass transformation and a MHT curve, which resembles the first EOF of MHT at least north of 20ЊN (solid line in Fig. 18c). As discussed in section 3b, the pattern of this first EOF presents strong anomalies of opposite signs associated with the LSW and STMW. Note that the negative anomaly associated with the STMW extends to the waters with characteristic temperature 16Њ-14ЊC (26.5 Ͻ 0 Ͻ 27). Thus STMW, formed in the Gulf Stream, cool as they move to the northeast within the NAC. This is a preconditioning phase for their later transformation into SPMW; which varies in phase with the formation of STMW. Nevertheless, the locally explained variance corresponding to SPMW ( 0 ϭ 27.1 at 12ЊC) is small and their contribution to the variability described by this canonical pair is not significant, which is consistent with the insignificant correlation between the NAO and the formation rate of these mode waters (Fig. 9c). The 3yr delayed response of the MHT to this buoyancy forc-  ing is associated with negative anomalies south of 40ЊN and large positive anomalies at 48ЊN (Fig. 18c), and is likely linked to the variability of LSW. Periods of the maximum formation of LSW corresponding to high NAO index and deep MLD in the Labrador Sea (1972-76, 1982-84, and 1989-95), are followed nearly three years later (1976-80, 1986-87, 1991-97) by an increase in downwelling at 50Њ-60ЊN, as suggested by the first EOF of the overturning (Fig. 15b) and its PC (Fig. 16a), a stronger surface northward flow north of 42ЊN and a significant increase of the northward MHT.
The second canonical pattern of water mass transformation at lag 3 (Fig. 18b) is presented by positive anomalies of transformation of LSW and STMW and negative anomalies of the surface transformation of the SPMW. Compared to the first pair, the peculiarity of the second pair is to have the LSW and STMW in phase, and the SPMW in opposition. Again the peak of explained variance relative to the SPMW is small. This pattern is optimally correlated with the MHT curve, resembling the second EOF of MHT (Fig. 18c), which appears as the dominant mode in the simulation of Hakkinen (1999). This mode has no impact north of 50ЊN, and is only relevant to the variability in the subtropical gyre. It represents an increase (or decrease) in the MHT over the subtropical gyre (10Њ-50ЊN), related to a positive (or negative) anomaly in overturning at the same latitudes. The corresponding pattern of the overturning streamfunction shows an intensified downwelling between 45Њ-55ЊN (region of formation of SPMW), and increased upwelling between 15Њ-35ЊN. The corresponding changes in the buoyancy forcing, (positive anomalies in the formation of LSW and STMW, and negative anomalies in SPMW), occur three years ahead of the anomalies in the overturning and MHT. Although we have no clear physical explanation for the moment, we can hypothetically assume that the 3-yr delay is the time needed for the baroclinic adjustment of the western subtropical gyre to the changes in STMW properties.

3) WATER MASS TRANSFORMATION AND MHT AT 7-YR LAG
The first canonical pair at lag 7 shows the second highest correlation among the lags studied (0.87 in Table  2). The corresponding canonical pattern of the water mass transformation (Fig. 18d) is characterized by a general weakening of the transformation rate of all mode waters, associated with the decreased density flux over the northern North Atlantic. The associated MHT pattern (Fig. 18f) represents a general intensification of MHT at all latitudes. Thus the warming of the ocean surface at all latitudes results in a general increase of the northward MHT. The second canonical pair at lag 7 (Fig. 18e) is largely the same as the first one at lag 3 (Fig. 18a)-that is, a subpolar-subtropical dipole in the surface buoyancy forcing and positive anomalous MHT between 30Њ and 60ЊN. However, both 7-yr responses should be considered with caution because they may be influenced by autocorrelation in the MHT and the forcing function.

Summary and discussion
We have described the results of the numerical simulation of the large-scale interannual variability of the North Atlantic circulation in a numerical simulation carried out with the SPEM model driven by NCEP-NCAR reanalysis surface fluxes for the period 1958-97. Contrary to the experiments of Battisti et al. (1995), Hakkinen (1999), Eden and Willebrand (2001), Eden and Jung (2001), our formulation of the forcing includes a relaxation to the observed monthly SST and a weak relaxation to the seasonal SSS climatology. Since it is part of the forcing, the model SST is a strongly constrained parameter. Nevertheless, the formulation of the forcing used here allows the model to be prognostic in most other respects, and our approach, which focuses on the relationship between water mass transformation rates and meridional transports, is consistent with the way the forcing was applied.
In support of our choice of forcing, it is not clear whether applying a forcing function without any feedback to the atmosphere, as in the experiments cited above, is a better choice for the temperature. Indeed, in nature the SST is not determined by externally prescribed fluxes. The fluxes arise from coupled interactions between the ocean and the atmosphere. Ferry (2001) showed that fixed heat fluxes (with no parameterization of any feedback) tend to overestimate the magnitude of the SST changes due to the lack of the mechanism for SST moderation. In nature, the SST changes are moderated by the direct feedback action of the SST on the fluxes. On the other hand, since the SSS has no direct effect on the surface fluxes, it seems reasonable to use as weak a relaxation to the observed surface salinity as possible. We were able to get a quite reasonable simulation with a weak relaxation to the surface salinity.
Although the surface density input diagnosed by the model is closely correlated with the original NCEP-NCAR surface fluxes, it exhibits some differences in the 1960s and 1990s and demonstrates closer links with the NAO in comparison to the original surface heat flux. Interannual changes in the buoyancy fluxes considered in the T-S plane identify leading modes of variability associated with the low frequency coordinated changes in the transformation of surface waters into the three major types of mode waters observed in the North Atlantic (STMW, SPMW, and LSW). The dominant role of these mode waters is realistic and consistent with the known features of the North Atlantic circulation. It appears very clearly in the numerical experiment because the modeled large scale circulation patterns (e.g., the location of the Gulf Stream and the NAC) are consistent with the large scale patterns of the forcing components.
Joint statistical analyzes of the surface forcing characteristics and parameters of meridional circulation allowed us to link the intensity of MHT with surface buoyancy and wind forcing. To the extent that it is possible to make physical inferences from the statistical approach applied here, we draw the following conclusions.
• The intensity of the ocean meridional circulation has an immediate response of a barotropic nature to the wind stress curl, dominated by coordinated changes induced in the barotropic circulation of the western subpolar and subtropical gyres. This is consistent with the results of Eden and Willebrand (2001), who found primarily negative correlation between MHT at 48ЊN and the NAO index in their experiments forced with wind anomalies. The process involved, according to Eden and Willebrand (2001), is a topographic Sverdrup balance. Our model responds according to the same principle, but due to differences in the position of the major currents (Gulf Stream and NAC) with respect to the forcing fields, the pattern of the barotropic streamfunction is slightly different, explaining the quantitative discrepancies. In the tropics and subtropics MHT changes are largely dominated by the Ekman transport, whose variability is comparable to the total MHT variability (Eden and Willebrand 2001;Hakkinen 1999). The immediate response of the barotropic streamfunction to the changes in the wind stress curl shows some peculiarities in low latitude regions [e.g., a stronger deep western boundary current, and a stronger North Brazil Current (NBC)], indicating an increase of the meridional overturning. This is consistent with positive anomalies in MHT having the largest amplitude in the NBC region. • The dominant baroclinic response seems to be associated with the baroclinic adjustment of the subpolar gyre to the changes in stratification related to the formation of LSW induced by anomalous density fluxes, delayed with a 3-yr lag, likely through a joint effect of baroclinity and relief (JEBAR) effect or boundary wave propagation. This response involves extremes in the formation of LSW and STMW, associated with the MHT enhancement at 45ЊN (Fig. 15a). The associated overturning variability pattern (Fig. 15b) shows a surface convergence at 45ЊN (the boundary separating the subtropical and the subpolar gyre), and an increase of the southward deep flow, which is felt only to the south of 35ЊN (within the subtropical gyre). This suggests that the mechanism involved is not an acceleration of the overturning cell by entrainment of dense water in the subpolar gyre (which should drive an increase of the deep southward flow from high latitudes), but a geostrophic adjustment within both gyres to the changes in stratification driven by the STMW and LSW formation process. This is slightly different from Eden and Willebrand (2001) who also explained the model response at 3-yr lag by a delayed baroclinic adjustment of the subpolar gyre, mainly in response to the wind forcing, the contribution of the buoyancy forcing appearing of secondary importance. It is likely that our model gives more importance to the buoyancy forcing due to a stronger JEBAR effect in sigma models (Penduff et al. 2001). Nevertheless, it is a robust feature that the baroclinic adjustment of the gyres induces a large positive correlation at 3-yr lag between the NAO index and the MHT at 48ЊN. Note that Groetzner et al. (1998) with a coupled-ECHO model also showed that heat flux leads the oceanic SST changes by 2-3 yr. • A 3-yr delayed response of the MHT in the subtropical gyre to the buoyancy forcing is the second-order mode in our simulation. This agrees in terms of the time lag with Hakkinen (1999). A comparison of our MHT (Fig. 12) with those of Hakkinen (1999) shows that positive MHT anomalies in the subtropics in the mid-1970s and in the midlatitudes in 1990s were observed in both simulations. On the other hand, there is little agreement for the period from the late 1970s to the late 1980s. The structure of Hakkinen's (1999) MHT anomalies indicates coordinated changes at all latitudes with the strongest anomalies at 25Њ-35ЊN, whereas in our simulation we observe an out-of-phase relationship between subtropical and subpolar latitudes. This suggests that the major mode of variability simulated by Hakkinen (1999) is largely our second leading mode of MHT. It accounts in our simulation for only 21% of the total variance, but dominates in the Hakkinen (1999) simulation. Without a specific intercomparison of Hakkinen's and our models we can only speculate about the reasons for these differences. Reasons could lie in many model parameters or details of the model configuration or forcing. Our intuition is that the formation of LSW at high latitudes, driving the baroclinic response to the buoyancy forcing in our simulation, differs significantly between the two experiments.
General intensification of MHT over all latitudes occurs 7 yr after the decrease of the transformation rates of all mode waters, at a timescale needed for the baroclinic adjustment of the subtropical gyre. The first mode of the 7-yr delayed response to the surface heat forcing, associated with a strengthening or weakening of the northward transport is consistent with the temporal scale of the baroclinic adjustment of the great subtropical gyre. A secondary mode of variability is similar to the first mode found at the 3-yr lag. However, we suggest that one consider 7-yr delayed responses with care, since they can be influenced by the self-correlation of the MHT and forcing functions. Thus, our model responds in a way consistent with the findings of Eden and Willebrand (2001), and with some of the findings of Hakkinen (1999). Therefore, the mechanism of the MHT variability shows some robustness to model characteristics and forcing formulation. Using a statistical analysis, it is difficult to speculate about delayed response to the wind stress, obtained by Eden and Willebrand (2001), as well as about the immediate response to surface buoyancy forcing. Our results (Figs. 14a,c) show an increase of correlation for both cases. For the wind stress curl the sign is reversed as in the experiment of Eden and Willebrand (2001). However, these responses may be modulated by the strong correlation of the wind stress with the heat flux (and correlation of both with NAO). To isolate clearly the different factors, it would be necessary to make simulations with idealized forcings (e.g., Eden and Willebrand 2001). In particular in the Groetzner et al. (1998) experiment wind stress curl also leads the SST changes by approximately 3 yr, although, as it was pointed out above the direct comparisons of our results, based on the strongly constant SST, with those obtained using coupled models, is difficult.
Comparison of our results (as well as the results of the other simulations) with experimental data may be difficult to interpret, because most quantities being compared are, in the experimental works, not measured but computed from observations. These computations also have some limitations and their accuracy still remains quite low, first of all because of irregularity of sampling at the ocean cross sections. This leads to inconsistency of estimates of baroclinic and barotropic MHT with the Ekman MHT, which is usually taken as climatological in many works of practical oceanographers. Recently Koltermann et al. (1999) and Lorbacher (2000) derived MHT estimates at 48ЊN using pre-WOCE (World Ocean Circulation Experiment) data for 1957 and 1982 and WOCE sections for the period 1992-98. Quantitative comparison of their MHT estimate at 48ЊN with our model estimates shown in Fig. 19a reveals a first inconsistency: the magnitude of interannual changes is much smaller in the model than that reported by Lorbacher (2000). Koltermann et al. (1999) and Lorbacher (2000) report quite large error bars, ranging from 1.1 to 1.5 ϫ 10 14 W, for their estimates. In order to compare general tendencies in the simulated and measured MHT at 48ЊN, we normalized both time series with respect to the standard deviations (Fig. 19b). We can conclude that some of the very general tendencies are qualitatively comparable. Further development of this study can go in a number of directions. The respective role of the thermal and haline forcing could be investigated from sensitivity experiments with only one component of the forcing varying. But beside new simulations with different models, forcing functions, and formulations of forcing, it will be of great importance to understand the respective role of surface forcing and internal model variability in the circulation changes simulated. Of special interest is the behavior of the Gulf Stream, and particularly its location in response to the wind and buoyancy forcing (Marshall et al. 2001;Taylor and Stevens 1998). If our model qualitatively produces a wind-driven ocean response whose pattern is very consistent with the interaction scenario between the NAO and the ocean circulation proposed in the idealized study of Marshall et al. (2001), its robustness still has to be confirmed due to small variability of the limit between the major ocean gyres noticed in coarse-resolution experiments. In the ocean models these phenomena can be effectively studied with high-resolution simulations, which allow us to identify variability in the Gulf Stream location. Simulations for longer time periods can be very helpful in both understanding the mechanisms and quantifying the variability on the timescales of several decades. In this context the experiment of Eden and Jung (2001) with the reconstructed NAO-like forcing for 130 yr may set a reference for further experiments with other models. Due to the chosen forcing formulation we did not study in this work the response of the SST variability to the changes in the North Atlantic large-scale circulation. However, positive NAO should result in a heat convergence around FIG. 19. (a) Variability of the actual values of MHT at 48ЊN in the model simulation (thin line, annual means; bold line, 3-yr running mean) and from experimental works of Koltermann et al. (1999) and Lorbacher (2000) (squares) and (b) anomalies of MHT at 48ЊNi n the model experiment normalized with respect to their standard deviations (thin line, annual means; bold line, 3-yr running mean) and derived by Lorbacher (2000) from the hydrographic sections (squares). 40ЊN (Fig. 15), that leads to the SST increase in the region of convergence. Besides the coupled models, SST responses (including those that can affect the forcing fields) can be partly studied in simulations with a forcing formulation, that accounts for the feedbacks with the atmosphere (e.g., Rahmstorf and Willebrand 1995;Seager et al. 1995).