## Abstract

In the Stommel box model, the strength of the overturning circulation is parameterized in terms of the density (and hence the pressure) difference between the two boxes. Straub has pointed out that this parameterization is not consistent with the Stommel–Arons model for the abyssal circulation. In particular, the zonally averaged density field implied by the Stommel–Arons model is unrelated to the strength or the direction of the meridional overturning circulation. Here, the inconsistency is examined using the abyssal circulation model of Kawase and a variant to include the effect of Southern Hemisphere wind forcing. The important parameter is *R,* the ratio of two timescales: the timescale for a perturbation to the density field to propagate, by either wave or advective processes, from a high-latitude source to the equator and the timescale for the dissipation of a perturbation to the density field by diapycnal mixing. If the model is forced only by a deep water source in the northern basin, it is found that the model behaves like the Stommel–Arons model when *R* ≪ 1 (the “weak” damping regime) and like the Stommel box model when *R* ≫ 1 (the “strong” damping regine). Estimates of *R* suggest that coarse-resolution models generally reside in or near the Stommel box model regime (*R* ≫ 1), which is probably why these models generally support the Stommel box model hypothesis and corroborate the momentum-based closure used in zonally averaged models. On the other hand, it is not clear that the real world is also in the strong damping regime. Indeed, it is easy to obtain estimates for *R,* using realistic parameter values, that sit in the weak damping regime. It is shown that, even in the weak damping regime (*R* ≪ 1), adding forcing by the Southern Hemisphere circumpolar westerlies generally moves the model into the Stommel box model regime. It therefore is concluded that, at least in the context of the Kawase model, the inconsistency noted by Straub can be removed by including the effect of Southern Hemisphere wind forcing and that the Stommel box model approach probably has wider applicability than is suggested by estimates of *R* alone.

## 1. Introduction

The meridional overturning circulation is an important mechanism for the transport of heat from low to high latitudes (Hall and Bryden 1982; Roemmich and Wunsch 1985). Our understanding of the overturning circulation has its roots in Stommel's classic papers, notably Stommel and Arons (1960) and Stommel (1961). In his 1961 paper, Stommel proposed a simple two-box model driven by buoyancy fluxes to show that multiple steady states can exist under identical forcing. A key assumption in the Stommel box model is that the strength of the overturning circulation is proportional to the density difference between the two boxes. In the Stommel and Arons (1960) model, a very different approach is taken. Here, the strength and structure of the overturning circulation is specified and the model solves for the horizontal circulation. Straub (1996) has pointed out a fundamental inconsistency between the two approaches. By averaging the Stommel–Arons solution over the regions corresponding to Stommel's boxes, Straub compared the relationship between the meridional overturning circulation and the density field in the Stommel and Arons model. He found that the zonally averaged density field in the Stommel–Arons model is independent of both the sense and structure of the meridional overturning cell or cells, in direct contradiction with the Stommel box model.

The Stommel box model has inspired many studies of the thermohaline overturning circulation and its stability using simplified models. These include more complex box models and zonally averaged models [see Weaver and Hughes (1992) for a review]. A key question in understanding the dynamics of the overturning circulation is how to establish the east–west pressure difference across an ocean basin that is required to geostrophically balance the overturning, especially given the predominently equator/pole structure of the buoyancy forcing (Wajsowicz and Gill 1986; Greatbatch and Peterson 1996). This issue is of great importance for zonally averaged models. Indeed, as noted in the comprehensive review by Wright et al. (1998), zonal averaging does not, of itself, provide information on the east–west pressure difference across an ocean basin. Rather, a parameterization of the east–west pressure difference in terms of zonally averaged quantities is required. Furthermore, it is clear that the closure adopted for the east–west pressure difference is a key aspect of these models, with possible important consequences for the model behavior, including the sensitivity of the modeled overturning to changes in surface forcing, an issue central to the ocean's role in climate. A thorough review of the closures used in zonally averaged models can be found in Wright et al. (1998). These authors note that there are two kinds of closures, one based on the momentum equations (e.g., Marotzke et al. 1988; Wright and Stocker 1991: see section 5 for a detailed discussion) and the other based on the vorticity equation (Wright et al. 1995). The former type of closure is consistent with the Stommel (1961) hypothesis and is a local parameterization that relates the strength of the overturning circulation at a given latitude directly to the north–south gradient of the zonally averaged density field at the same latitude. The vorticity-based closure, on the other hand, is much more akin to the Stommel and Arons approach, including an interior circulation based on Sverdrup dynamics. To complete the closure, it is necessary to apply a boundary condition at the basin extremities, a boundary condition that Wright et al. (1998) describe as “somewhat uncertain.” Wright et al. (1995) made a plausible choice for this boundary condition, and Wright et al. (1998) provide some physical justification for it. In test cases, Wright et al. (1995) find generally good agreement between model results obtained using the two closures, although the role of the uncertain boundary condition remains a concern. Evidence from three-dimensional coarse-resolution models has been used to support the Stommel box model hypothesis. For example, Rahmstorf (1996) finds a linear relationship between the flow of North Atlantic Deep Water (NADW) in his model and the density difference between the NADW formation and outflow regions. In addition, the 3D model calculations of Hughes and Weaver (1994) have been used to calibrate both the momentum and vorticity based closures for zonally averaged models (Wright et al. 1998). Using a model based on the momentum closure, Wright and Stocker (1992) find surprisingly good agreement between model results and observed zonally averaged water mass properties, suggesting that the Stommel box model concept is applicable to the real ocean.

While the overturning circulation is generally thought to be buoyancy-driven (Stommel 1961), there is growing interest in a possible link, originally suggested by Toggweiler and Samuels (1993), between the wind stress at the latitude of Drake Passage and the strength of the overturning circulation in the North Atlantic. At the latitude of Drake Passage, strong westerly winds drive a northward Ekman drift. The water pushed north in this way has to return southward at depth to maintain mass balance. In the absense of diapycnal mixing to convert surface water to deep water (apart from in the deep convection zones), the northward Ekman drift that enters the South Atlantic Ocean has to flow all the way to the deep convection sites of the North Atlantic Ocean before being returned at depth. Toggweiler and Samuels (1995) suggest that both the deep outflow of NADW from the Atlantic basin and the production of NADW in northern high latitudes vary in direct proportion to the strength of the northward Ekman drift at the latitude of Drake Passage. McDermott (1996) investigated the link between Southern Hemisphere wind forcing and the northern overturning in a study using idealized geometry. His work demonstrates the importance of the boundaries for communicating information about the wind forcing in the Southern Hemisphere of his model domain to the Northern Hemisphere sinking region. Rahmstorf and England (1997) have argued that including feedback on the atmospheric temperature seen by a model [as e.g., in Zhang et al. (1993) or Rahmstorf and Willebrand (1995)] acts to reduce the importance of Southern Hemisphere winds compared to the local production of NADW at northern high latitudes. They argue that, at least in coarse-resolution models, some of the northward Ekman transport at the latitude of Drake Passage can be converted to deep water within the Southern Hemisphere, just north of Drake Passage. Gnanadesikan (1999) discusses the role of Southern Hemisphere wind forcing in driving northern overturning in a simple box model formulation. His study includes the effect of mesoscale eddies in the Southern Ocean. He notes that the eddy-driven mass transport is expected to oppose the northward Ekman transport at the surface and reduce the magnitude of the northward flow out of the Southern Ocean below that predicted by wind-driven Ekman theory. On the other hand, using a combination of observations and simple theory, Speer et al. (2000) argue that the eddies drive a southward, along-isopycnal return flow within the thermocline that feeds the northward Ekman transport at the surface where the isopycnals outcrop (see their Fig. 8). The idea that the wind-driven and buoyancy driven circulations are intimately linked is also highlighted by Nof (2000). In his paper, Nof has shown, subject to certain assumptions, that the strength of the North Atlantic overturning can be expressed entirely in terms of the surface wind stress.

In this paper, we modify Kawase's (1987) abyssal circulation model to mimic the effect of Southern Hemisphere wind forcing. The emphasis here is on simplicity. We seek to illustrate the essential dynamical processes using as simple a model as possible, in this case a shallow-water ocean model. The model is formulated in section 2. In section 3, we use the model to investigate the inconsistency between the Stommel–Arons model and the Stommel box model in the context of closed basin geometry and buoyancy forcing alone. Then in section 4, we investigate the impact of adding a zonally reentrant channel in the Southern Hemisphere and including forcing to mimic the Southern Hemisphere circumpolar westerly winds. We argue that adding forcing by Southern Hemisphere winds can move a model toward the box model regime, even when the purely buoyancy forced problem sits closer to the Stommel–Arons regime. Section 5 discusses our results in the context of the closure used in zonally averaged models that is based on the momentum equations. In section 6 we use the model to provide a simple illustration of how Southern Hemisphere wind forcing can drive northern overturning. Finally, section 7 provides a summary and conclusions.

## 2. Model formulation

We use a linear shallow-water model formulated for the abyssal ocean, following the basic approach of Kawase (1987). The governing equations, written for an equatorial beta plane, are

where *g*′ is reduced gravity, *h* is the interface displacement (positive upward), and *β* = 2Ω/*R*_{e} is the gradient of the Coriolis parameter at the equator. The continuity equation allows for a deep-water source *Q* and, following Kawase (1987), parameterizes the cross-interfacial flow (other than that specified in the source region) using a Newtonian damping term, –*λh.* It should be noted that in steady state, positive *λh* corresponds to upward vertical motion, negative *λh* to downward vertical motion. In some of the model runs, *λ* = *λ*(*x, **y*) is a specified function of the spatial coordinate, allowing for spatial inhomogeneity in the diapycnal mixing. The deep-water source *Q,* on the other hand, parameterizes the addition of mass to the abyssal layer by means of deep convection. Another important parameter is the Rayleigh friction *K.* Unless otherwise stated, this parameter has value *K* = 1.3 × 10^{–6} s^{–1}, corresponding to a Stommel boundary layer of width 130 km for *β* = 1. × 10^{–11} m^{–1} s^{–1}. The stress term, *τ*^{x}/(*ρ*_{0}*H*) mimics the effect of Southern Hemisphere wind forcing (see below). For simplicity, the northward wind stress, *τ*^{y}, has been set to zero. The stress is modeled using a Gaussian profile in the meridional coordinate, centered about a specified latitude, *y* = *y*_{0}, with an *e*-folding scale of 15° latitude:

Since our model is for the abyssal circulation, the stress in the model acts westward (*τ*_{0} < 0), as described below, even though the actual surface wind stress acts toward the east.

A detailed derivation of Eqs. (1)–(3), including the stress and Rayleigh friction terms, is given in appendix A for the case of a two-density layer ocean. Here, we summarize the argument. We begin by noting, following Rintoul et al. (2001), that in the Southern Ocean, the zonal wind stress at the surface is thought to be balanced by topographic form stress acting in the opposite direction at the bottom. This balance is expressed, in the zonal average, by Eq. (A10) in appendix A. In the two-density layer model, the eastward wind stress acts as a body force over the top layer, and the topographic form stress acts in the opposite direction (i.e., westward), on the lower, abyssal layer. The vertical flux of momentum is represented by an interfacial form drag, resulting in the Rayleigh friction term in (1) and (2). Since there are two density layers, there are two vertical modes: the barotropic mode and the baroclinic mode. In the former, the horizontal velocity corresponds to the vertically averaged velocity, and the projection of the wind stress forcing on this mode acts eastward with the topographic form stress term acting westward (in equilibrium, these almost balance). For the baroclinic mode, the horizontal velocity corresponds to the difference in velocity between the two layers. The projection of both the wind forcing and the topographic stress terms on this mode acts westward (eastward) if the velocity is interpretated as the baroclinic contribution to the velocity in the abyssal (surface) layer. The model equations (1)–(3) being considered here correspond to the baroclinic mode, and the stress term acts westward in the model since the velocity in the model is interpreted as that in the deep, abyssal layer. It should be noted that the stress term is a linear combination of the wind stress and topographic form stress terms. The relative importance of the two contributions depends on the vertical structure of the thermocline. Since, however, both the surface wind stress and bottom, topographic form stress are expected to roughly balance each other, their relative contribution to *τ*^{x}/(*ρ*_{0}*H*) is not important here.

The argument given in appendix A can be generalized to continuous stratification following the argument of McCreary (1981). McCreary described a generalization of the concept of vertical normal modes (Gill 1982) to include vertical mixing of momentum and density. When the vertical mixing coefficients are chosen to be inversely proportional to the square of the buoyancy frequency (e.g., as for the form stress term in Greatbatch and Lamb 1990), then, using a slightly modified form for the vertical mixing of density, the separation into vertical normal modes proceeds as in the standard case. The vertical mixing terms appear as Rayleigh friction and Newtonian damping terms in the resulting shallow-water equations describing each vertical mode, exactly as in Eqs. (1)–(3), and a forcing term appears in the equation for each vertical mode, analogous to the stress term in Eq. (1). McCreary's approach therefore provides justification for using a Newtonian damping term to represent vertical (i.e., diapycnal) mixing of density and a Rayleigh friction term to represent vertical mixing of momentum. The association between eddy form stress and vertical mixing of momentum has been discussed by Rhines and Holland (1979), Greatbatch and Lamb (1990), Gent et al. (1995), and Greatbatch (1998).

Unless otherwise stated, the baroclinic layer depth *H* is set at 400 m and the reduced gravity *g*′ = 0.02 m s^{–2}. The corresponding gravity wave speed is 2.5 m s^{–1}, which can be thought of as corresponding to the first baroclinic mode in a continuously stratified ocean. The model actually integrates the spherical geometry form of Eqs. (1)–(3). We use finite differencing on an Arakawa B grid (Mesinger and Arakawa 1976). The numerical resolution is 1° × 1°.

In the discussion that follows, zonally averaged quantities are represented by an *h*; for example, *h* is the zonal average of the interface height. Also, since we are working with a reduced-gravity model, the velocity is interpreted as that of the first baroclinic mode so that, for the purpose of volume conservation, the volume flux in the deep layer, given by *V* = *W**υ* (*W* is the east–west width of the ocean), is compensated by an implicit, opposite volume flux in the overlying layer: *V* is, therefore, a measure of the strength of the meridional overturning circulation in the model. Similarly, *h* can signify either the pressure in the water column in the lower dense layer or the density anomaly. Therefore, the zonal average of the interface height *h* corresponds to either the zonally averaged density or pressure (Straub 1996).

## 3. Experiments in a rectangular basin using buoyancy forcing only

We begin by describing several experiments in a rectangular, closed basin that is symmetric about the equator. In this section, the model is driven by a source, *Q,* in the continuity equation and the stress term in (1) is set to zero (the influence of wind forcing is discussed in section 4); *Q* has a Gaussian structure, centered in the northwest corner of the domain with an *e*-folding scale of 5° in both latitude and longitude. Our objective here is to investigate the effect of linear damping (–*λh*) on the relationship between the meridional overturning and the north–south gradient of the zonally averaged density field. The damping coefficients *λ* are spatially uniform with values of 1 × 10^{–9}, 3 × 10^{–8}, 6 × 10^{–7}, and 1.2 × 10^{–6} s^{–1}, corresponding to damping timescales of roughly 30 yr, 1 yr, 20 days, and 10 days, respectively (the experiments are summarized in Table 1). Figure 1 shows plan views of the interface height in the final equilibrium state of experiments I–III. The flow pattern in the weak damping case (expt I) (which can be deduced to a good approximation from geostrophy) is very like the Stommel–Arons solution (Kawase 1987) with strong boundary currents and interior poleward and eastward flow, in both hemispheres, and with the western boundary current reaching deep into the Southern Hemisphere (as indicated by the gradient of the interface displacement across the western boundary layer, implying southward flow to near 30°S). With increasing *λ,* the flow is increasingly confined to the Northern Hemisphere western boundary and the equatorial regions (Fig. 1c). Kawase (1987) has given a thorough explanation, both physically and theoretically, for the dependence of the solution on *λ.* Of interest here is how the zonally averaged overturning, *V,* is related to the zonally averaged interface displacement *h*. The profiles of *h* and *V* for each experiment are shown in Fig. 2. The weak damping case (expt I) is closest to the Stommel and Arons solution. The structure of *h* in this case is almost symmetric about the equator, apart from the source region and near the southern boundary. The symmetry arises first from the Kelvin wave adjustment around the model boundaries, and then from the Rossby wave adjustment expanding the contours of *h* into the interior on both sides of the equator. The structure and even the sense of *V* are almost completely independent of *h*, just as described by Straub (1996). In particular, there is no local relationship between *V* and *h* as implied by the Stommel box model approach or the momentum closure used in zonally averaged models (Wright et al. 1998).

With increasing *λ,* the symmetry in *h* between the hemispheres is increasing lost, as first the Rossby wave adjustment and then the Kelvin wave adjustment is increasingly arrested. In Fig. 1c, the Kelvin wave excited by the source in the northwest corner of the domain is strongly damped and does not succeed in making it far into the Southern Hemisphere along the eastern boundary. In this case, the zonally averaged north–south flow is down the zonally averaged pressure gradient (Fig. 2c) over most of the domain, consistent with the Stommel (1961) box model, but there is an exception near the equator where the relationship breaks down. Experiments using even stronger damping (e.g., Fig. 2d) show that the region of discrepancy is removed only when the damping is so strong that the Kelvin wave is finally prevented from reaching the equator, emphasizing the importance of the Kelvin wave in our model experiments for determining whether the box model relationship is valid or not. We conclude that the Kawase model agrees with the Stommel–Arons model in the limit of “weak” Newtonian damping (implying weak diapycnal mixing), but agrees with the Stommel box model when the Newtonian damping is “strong,” corresponding to relatively strong diapycnal mixing, and that the Stommel box model regime starts to break down as soon as the damping is weak enough to allow the Kelvin wave to reach the equator.

Let us now examine more closely what we mean by “strong” and “weak” in the above. In the simple model we have used here, information about the northern source, *Q,* is communicated toward the equator by Kelvin waves propagating along the western boundary. Döscher et al. (1994) found similar behavior in both 1° and a ⅓° (eddy permitting) models of the North Atlantic Ocean in response to a change in the high latitude boundary condition. Hsieh and Bryan (1996) have suggested that sea level changes associated with global warming may be communicated to different parts of the global ocean by means of baroclinic Kelvin waves propagating around the ocean boundaries, followed by Rossby wave propagation into the ocean interior. Huang et al. (2000) have also argued that wave processes are important for communicating information about changes in the rate of deep-water formation in the North Atlantic to other ocean basins, and the model of Johnson and Marshall (2002) also relies on wave propagation along the ocean boundaries to communicate changes in deep-water formation from high to low latitudes and across the equator. Further, McDermott (1996) has argued that in his model, waves propagating along the model boundaries communicate the information about Southern Hemisphere wind forcing to the northern high-latitude deep convection regions. Several authors (e.g., Winton 1996; Greatbatch and Peterson 1996) have also noted that boundary waves play an important role in the mechanism of self-sustained interdecadal variability in flat-bottomed ocean models. However, there is evidence from 3D primitive equation models that advective as well as wave processes can be important for communicating information from high-latitude source regions to lower latitudes, for example, Gerdes and Köberle (1995), Goodman (2001), Marotzke and Klinger (2000), and Eden and Greatbatch (2003). For an advective process, the southward propagation speed is considerably reduced compared to the Kelvin wave speed typically used in the Kawase model. For example, in the model of Marotzke and Klinger (2000), perturbations in the vertical velocity take roughly 45 years to propagate to the equator, implying a propagation speed roughly two orders of magnitude smaller than the Kelvin wave speed in our model experiments. However, the essential physics is not changed by the reduced propagation speed. Perturbations to the density field that are generated in high latitudes are subject to dissipation by diapycnal mixing as they propagate equatorward, whether this be by advective or wave processes, exactly as the perturbations in interface height *h* are subject to dissipation by the Newtonian damping term in our model.

As shown in appendix B [see Eq. (B6)], for our model experiments (in which the Rayleigh friction *K* is held fixed), the relevant parameter is

where *L* is the distance along the western boundary from the high-latitude source to the equator and *c* = is the propagation speed for Kelvin waves; *R* is the ratio of the propagation timescale (from the high-latitude source to the equator) to the dissipation timescale associated with the Newtonian damping. When *R* ≫ 1, we are in the strong damping regime of the Stommel box model limit, as in experiment IV. When *R* ≪ 1, we are in the weak damping regime of the Stommel–Arons solution (e.g., expt I). Further, the model results show that the box model regime breaks down as soon as the Kelvin wave can reach the equator (that is *R* ≈ 1), as in experiment III. It is straightforward to extend the definition of *R* so that it is applicable to more complex models or to provide an estimate of the damping regime appropriate to the real world. For example, when advective processes dominate, the wave speed *c* is replaced by the advective speed *U.* In a similar vein, a more realistic estimate for the damping timescale 1/*λ* is the diffusive timescale *H*^{2}_{T}/*κ,* where *κ* is the diapycnal diffusivity and *H*_{T} is the thickness of the pulse of dense water that is propagating southward. So, as an alternative to Eq. (5), we have

Using Eq. (5), the values of *R* for experiments I, II, III, and IV are 0.003, 0.08, 1.4, and 2.8 respectively, corresponding to a range of values from the weak through to the strong damping regimes. For the experiment described by Marotzke and Klinger (2000), we use (6) to estimate *R* by taking the advective timescale *L*/*U* to be 45 yr, as implied by Fig. 7 in their paper, *κ* = 5 × 10^{–4} m^{2} s^{–1} (the value used in the boundary region of their model) and *H*_{T} = 100 m, at the lower end of possible estimates. This gives *R* = 71, which puts the model squarely in the strong damping regime of the Stommel box model. On the other hand, if we take *H*_{T} = 1000 m, at the large end of possible values, then *R* = 0.7, a value between that of our experiments II and III. Webb and Suginohara (2001) have noted there is considerable uncertainty in the value of the diapycnal diffusivity appropriate to the ocean. A value typical for the deep ocean, based on turbulence observations and dye diffusion, is 1 × 10^{–5} m^{2} s^{–1} (Gregg 1989; Ledwell et al. 1998) [although much larger values are found locally near rough topography (Polzin et al. 1997)]. If we replace *κ* = 5 × 10^{–4} m^{2} s^{–1} by *κ* = 1 × 10^{–5} m^{2} s^{–1} in the above estimate, we obtain *R* = 0.01, a value that is now more typical of the weak damping regime of the Stommel–Arons model. There is also evidence from observations (and also other models, e.g., Gerdes and Köberle 1995) that the advective timescale appropriate to the North Atlantic may be closer to 10 yr than 45 yr. For example, Molinari et al. (1998) found newly formed Labrador Sea Water appearing at 24°N roughly 10 years after its formation in the Labrador Sea. If we use 10 yr, rather than 45 yr, *R* is further reduced to 0.003, a value similar to that in our experiment I. If we use *H*_{T} = 100 m, rather than 1000 m, we can increase *R* to 0.3, a value between that of the weak and strong damping regimes. Nevertheless, a more rapid advection timescale, or a greater importance for wave processes, as in the model of Döscher et al. (1994), acts to reduce the effective value of *R,* and put the estimate firmly in the weak damping regime more typical of the Stommel–Arons model.

These estimates for *R* suggest that coarse resolution modeling studies, such as those used to fit the closures for zonally averaged models (Wright et al. 1998), are in the strong damping regime of the Stommel box model, which is why these models are able to provide credence to the momentum closure, in particular, used in zonally averaged models. On the other hand, it is not clear from the above analysis that the real world is also in the strong damping regime. Indeed, it is easy to obtain estimates for *R,* using realistic parameter values, that clearly sit in the weak damping regime, especially given that in the model experiments, the box model regime starts to break down for values of *R* ≤ 1.

## 4. Basin plus reentrant channel: Adding wind forcing

In this section the model domain (see Fig. 4) consists of a enclosed rectangular domain 60° wide from 30°S to 60°N plus a reentrant channel from 60° to 30°S that is also 2 times as wide as the enclosed part, in accordance with having a circumpolar region in the Southern Hemisphere. Two types of forcing are used to drive the model; a local deep water source, *Q,* and the stress term, *τ*^{x}/(*ρ*_{o}*H*) mimicking forcing by Southern Hemisphere winds (hereafter referred to as the wind forcing term; see the discussion in section 2). As noted in section 2, the meridional stress is set to zero and the zonal stress is given by Eq. (4) and is centered 20° of latitude to the north of the southern boundary. All the model experiments use the same parameter values as for experiment II in section 3, an experiment in which the model behavior is more akin to the Stommel–Arons model than the Stommel box model.

For the first run (expt QI), the model is forced by a deep-water source in the northwest corner of the domain, as in experiment II, and the wind forcing term is set to zero. In the second run (expt WI), the source term is set to zero and the model is driven by the wind forcing term: *τ*_{0} is chosen so that the Ekman transport at the northern edge of the channel is representative of that associated with the zonal wind stress at the latitude of the southern tip of South America. Integrated around the globe, this transport is about 34 Sv (Sv ≡ 10^{6} m^{3} s^{–1}), whereas integrated only over the longitudinal range corresponding to the Atlantic basin, it is about 9 Sv [based on the National Centers for Environmental Prediction (NCEP) climatology; L. Talley 2003, personal communication]. For experiment WI, we choose *τ*_{0} so that the Ekman transport in the model is 25 Sv when integrated around the latitude at the northern boundary of the channel. This choice lies between 9 Sv on the one hand and 34 Sv on the other. The sensitivity to the choice of *τ*_{0} is explored in later experiments.

Both cases are run to steady state and then zonally averaged. The profiles of *h* and *V* are shown in Fig. 3. In experiment QI, forced by the source in the northwest corner of the basin, there is no qualitative difference in the structure of *V* and *h* compared with experiment II (cf. Fig. 3a with Fig. 2b). Experiment QI therefore belongs to the weak damping regime. On the other hand, in the wind-forcing only case (expt WI), *V* flows down the gradient of *h* almost everywhere (except at the latitude of the channel, as we expect). In this experiment, information about the wind forcing is communicated into the northern basin through the propagation of Kelvin and Rossby waves, as described by McDermott (1996), with an anticyclonic gyre eventually being set up in the northern part of the basin (see Fig. 4b). Experiment QWI includes both the source and wind forcing terms and, since the model is linear, the solution is the sum of the solutions for experiments QI and WI. Here, the wind forcing dominates, and *V* again flows down the *h* gradient (Fig. 3d). We have tested the sensitivity of the solution to the relative strengths of the Northern Hemisphere source and the wind forcing. The results are shown in Fig. 5. Only in the case with the wind forcing reduced by one-half in comparison with experiment WI, and the strength of the source doubled, does the source finally become dominant over the wind stress, and *V* no longer flows consistently down the *h* gradient. These results show that by including forcing due to the Southern Hemisphere winds, the model can be pushed into the regime described by the Stommel box model, even though when the model is driven only by a northern source, it is in the weak damping regime, and hence closer to the Stommel–Arons regime. This is the main result of our paper.

Last, in this section, we note that for our results to hold up, it is important that the stress profile in our model intersect the northern boundary of the channel, as is easily demonstrated by moving the location of the maximum stress in Eq. (4) to the south. This is important in order to excite the Kelvin wave that communicates the information about the Southern Hemisphere wind forcing to the northern deep convection region, as in McDermott (1996).

## 5. Understanding our results and their relationship to the momentum-based closure used in zonally averaged models

To further understand our results and their relationship to closures used in zonally averaged models, we first consider the balance of terms in the steady, zonally averaged, meridional momentum equation

where *f* = 2Ω sin*θ* is the Coriolis parameter (*θ* is latitude). We can rewrite (7) as

where

and *u*_{g} is the zonal component of the geostrophic velocity. It is clear from (8) that for *υ* to flow down the gradient of *h*, *ε* must be positive. Indeed, the closure used in the zonally averaged model of Wright and Stocker (1991) (a typical example of the momentum closure discussed by Wright et al. 1998) is to assume that *ε* is a positive constant.

Another way to write (8) is

It is interesting to consider (10) in the context of our experiments. Since in all our experiments *υ* < 0, corresponding to a net southward transport at all latitudes *θ,* we shall restrict to *υ* < 0 throughout the following discussion. It follows from (10) that, for *ε* to be positive, *u*_{g} must change sign either side of the equator. Since for our model *f**u*_{g} = –*g*′*h*_{y}, we can deduce *u*_{g} from the plots of *h* in Figs. 2, 3, and 5. In the Stommel–Arons solution, and also our experiments I and II, *u*_{g} is dominated by the average of *u*_{g} over the interior of the basin (Straub 1996) and, as is clear from Figs. 1a,b and 2a,b, *u*_{g} is eastward in both hemispheres, meaning, from (10), that *ε* cannot be positive everywhere. Correspondingly, *υ* does not flow consistently down the gradient of *h* in these experiments. On the other hand, in experiment WI, *u*_{g} does indeed change sign either side of the equator, and *υ* flows down the gradient of *h* (Figs. 3b, 4b). The height contours in Fig. 4b show that the Southern Hemisphere basin (north of the channel) in experiment WI is dominated by the cyclonic (clockwise) gyre associated with the Sverdrup transport, with predominantly eastward interior geostrophic flow and *u*_{g} > 0. On the other hand, the northern basin is filled by an anticyclonic (also clockwise) gyre with predominantly westward interior geostrophic flow and *u*_{g} < 0. Both gyres are associated with the large-scale descent required to feed the Ekman suction in the southern part of the channel associated with the wind forcing. However, in the northern basin, the downwelling is stronger around the periphery of the anticyclonic gyre, as expected since the downwelling signal has been propagated along the equator and the model boundaries from the Southern Hemisphere and is then spread into the basin interior by Rossby waves [note that in Figs. 1 and 4, positive (negative) *λh* corresponds to upward (downward) vertical motion]. In the Southern Hemisphere, on the other hand, the downwelling is associated with downward Ekman pumping on the northern periphery of the band of zonal wind stress. In contrast to the northern basin, long Rossby waves radiated from the eastern boundary act to reduce the downward Ekman pumping, setting up the Sverdrup balance in their wake (Gill 1982, chapter 12), with the result that the strongest downward motion is now in the west, in the interior of the gyre, not around its periphery. It is this contrast in behavior between the northern and southern basins that is responsible for *u*_{g} having a different sign either side of the equator, and hence ensuring that *ε* is positive almost everywhere. It should be noted that the situation is different in experiments I and II. Here the influence of the forcing must be propagated into the interior by long Rossby waves in both hemispheres, with the result that *u*_{g} has the same sign in both hemispheres. It is also of interest to note that downwelling in both the northern basin and in the Southern Hemisphere north of the channel are important features of our model solution, suggesting a link to the work of Rahmstorf and England (1997) who noted that some of the northward Ekman transport at the latitude of Drake Passage is converted to deep water within the Southern Hemisphere, just north of Drake Passage, in their model. In those of our model experiments that include both a northern source as well as wind forcing (Fig. 5), then as long as the wind forcing dominates the effect of the source, an anticyclonic gyre is found in the northern basin, as in Fig. 4c, and the dynamics is essentially the same as in the wind-only case, with the result that *υ* flows down the gradient of *h*.

It is of interest to ask why the closure used by Wright and Stocker (1991) breaks down in experiments I and II. In deriving the Wright and Stocker (1991) closure, Wright et al. (1998) note that it is assumed that the zonal geostrophic velocity *u*_{g} does not change significantly over the width of the western boundary layer. Since *u* = 0 on the western boundary, it is easy to show that this assumption implies that *ε* in (9) is always positive, and thereby ensuring from (8), that *υ* flows down the gradient of *h*. In our experiments I and II, however, the assumption is not valid. As can be deduced from Figs. 1a and 1b outside the western boundary layer, *u*_{g} is eastward, as it is in the Stommel–Arons solution, and the associated north–south pressure gradient is equatorward. On the other hand, inside the boundary layer, the western boundary currents flows down the pressure gradient, southward from the high-latitude source. It follows that the north–south pressure gradient is reversed in the boundary layer and that *u*_{g} must also reverse in the boundary layer, in contradiction to the assumption that *u*_{g} remains essentially unchanged across the boundary layer. This is another way to understand why our experiments I and II, and also the Stommel–Arons solution, are not consistent with the Stommel box model or the momentum closure used in zonally averaged models.

## 6. A simple illustration of how Southern Hemisphere wind forcing can drive northern overturning

In this section, we use our model to illustrate how wind forcing over the Southern Ocean can influence northern overturning. Throughout the discussion, the northern source *Q* is set to zero. We begin by zonally averaging (3) to obtain

where *W* is the east–west width of the ocean and *V* is the volume flux associated with the overturning; (11) makes clear the close relationship between the meridional overturning and diapycnal processes (represented by *λ*). In particular, if *λ* = 0 everywhere, then *V*_{y} = 0 everywhere, and hence *V* = 0 everywhere, since *V* = 0 at the northern and southern boundaries, and there is no overturning. Nevertheless, a solution with wind forcing can be found when *λ* = 0 everywhere. Since *V* = 0, zonally averaging (1) shows that in steady state

Using the interpretation given in section 2, (12) implies that the combined stress forcing, from the zonal wind stress at the surface and the topographic form stress at the bottom, is exactly balanced by the form stress due to the eddies. It follows that (12) corresponds (for the baroclinic modes) to the dominant balance thought to operate in the Antarctic Circumpolar Current (Rintoul et al. 2001; Olbers and Ivchenko 2001).

Now consider experiment WII. In this experiment, the Newtonian damping coefficient *λ* is set to zero everywhere except in the far northern and southern parts of the model domain. The diapycnal mixing is therefore assumed to be zero everywhere except in two regions, one chosen to mimic the deep convection regions in the northern North Atlantic and the other the region of strong diapycnal mixing where isopycnals outcrop in the surface Ekman layer of the Southern Ocean. Since *λ* is no longer zero everywhere, a meridional overturning circulation is now possible, from (11). However, over the latitude band where *λ* = 0, *V*_{y} = 0. It follows that any northward transport where the basin joins the channel must extend northward all the way to the far northern region where *λ* ≠ 0. Indeed, as can be seen in Fig. 3c, the meridional circulation in the model extends all the way from the far southern to the far northern part of the model domain, connecting the two regions where *λ* ≠ 0. It is important to realize that to have a meridional circulation, we must have both regions where *λ* ≠ 0. Putting *λ* = 0 in the northern basin suppresses the overturning circulation there since, from (11), *V*_{y} = 0 and the boundary condition that *V* = 0 at the northern boundary then requires *V* = 0 everywhere in the northern basin. Similarly, putting *λ* = 0 in the circumpolar region requires *V* = 0 there, and hence at the southern end of the northern basin. Integrating (11) northward then requires that *V* = 0 everywhere in the northern basin.

We also carried out an experiment consisting of a two-basin system, in one of which *λ* ≠ 0 in the northern part of the basin but in the other *λ* = 0 everywhere. In this case, the meridional overturning connects the Southern Hemisphere to the basin in which *λ* ≠ 0. The implication is that the Southern Hemisphere wind forcing may be important for driving a deep overturning circulation in the Atlantic, rather than the Pacific, precisely because the Atlantic basin supports deep convection in northern high latitudes, whereas the Pacific basin does not. (Note that here, *λ* is mimicking the diapycnal mixing required to convert the northward upper thermocline flow that is driven by the Southern Hemisphere winds, to the southward deep return flow. To include a local conversion of water from the upper to the abyssal layer, such as can happen with deep convection, would require the addition of a source *Q,* as in the experiments discussed in section 4.) Of course, we are not saying that diapycnal mixing is zero in the North Pacific. It undoubtedly is not. But the presence of deep convection in high northern latitudes in the North Atlantic provides a mechanism to convert northward flowing near-surface waters to much denser North Atlantic Deep Water, a process that is not possible in the North Pacific.

As noted at the end of section 4, for the Southern Hemisphere winds to drive northern overturning, the wind stress profile in our model must intersect the northern boundary of the channel. This is important in order to excite the Kelvin wave that communicates the information about the Southern Hemisphere wind forcing to the northern deep convection region, as in McDermott (1996). The role of the Kelvin wave is clearly seen in the adjustment of the model to the final state, following the switch-on of the wind forcing.

Next, we note that, when *V* ≠ 0, (12) must be modified to

showing that the meridional circulation, which in view of (11) arises from diabatic processes (*λ* ≠ 0), implies an imbalance between the stress term, *τ*^{x}/(*ρ*_{o}*H*), and the eddy form stress. As in the argument given by Gnanadesikan (1999), this implies a net northward transport in the surface layers, where the surface wind stress wins out over the eddy form stress and a net southward transport in the deep layers where the topographic form stress wins out over the eddy form stress.

Last, it is clear from Fig. 3c that the relationship between *V* and *h*_{y} is weak in experiment WII, mainly because *h* also shows very little variation with latitude in the region where *λ* = 0. This can be understood, following the argument of Yamagata and Philander (1983) who note that in the presence of Newtonian damping with coefficient *λ* and Rayleigh friction with coefficient K, the equatorial radius of deformation scales as (*K*/*λ*)^{1/4}, and so is very large when *λ* = 0, but *K* ≠ 0 (see also the scaling given by Cane 1989).

## 7. Summary and discussion

We have used the Kawase (1987) abyssal circulation model, together with a variant to mimic forcing by Southern Hemisphere winds, to investigate the inconsistency noted by Straub (1996) between the Stommel (1961) box model and the Stommel and Arons (1960) model of the abyssal circulation. In the Stommel box model approach (e.g., the momentum closure used in zonally averaged models; Wright et al. 1998), the strength of the meridional overturning at a given latitude is directly related to the north–south gradient of the zonally averaged density at the same latitude. The important parameter is defined by Eqs. (5) and (6) and is the ratio of the timescale for a perturbation to the density field to propagate from a high-latitude source to the equator, to the timescale for the dissipation of a perturbation to the density field by diapycnal mixing. The propagation may be by wave processes (as in the Kawase model) or by advective processes, and we noted that results from coarse-resolution ocean models suggest advective processes can be important (e.g., Gerdes and Köberle 1995; Marotzke and Klinger 2000; Goodman 2001; Eden and Greatbatch 2003).

When the model is forced only by a deep-water source in the northern basin (e.g., associated with high-latitude deep convection), we showed that the model behaves like the Stommel–Arons solution when the diapycnal mixing is “weak” (*R* ≪ 1) and like the Stommel box model when the diapycnal mixing is “strong” (*R* ≫ 1). Furthermore, in the Kawase model, the Stommel box model assumption starts to break down as soon as the Kelvin wave excited by the northern source can reach the equator (i.e., for *R* ≤ 1). Estimates of *R* suggest that coarse-resolution models generally reside in or near the Stommel box model regime (*R* ≫ 1), which is probably why these models have been successful in corroborating the momentum closure used in zonally averaged models (see Wright et al. 1998). On the other hand, it is not clear that the real world is also in the strong damping regime. Indeed, it is easy to obtain estimates for *R,* using realistic parameter values, that clearly sit in the weak damping regime. However, adding forcing by the eastward wind stress associated with the Southern Hemisphere circumpolar westerlies puts the model into the Stommel box model regime, even when *R* ≪ 1 (although it should be noted that there is some dependence on the relative strength of the northern source and the wind forcing, as illustrated in Fig. 5). Wright and Stocker (1992) found remarkable agreement between a zonally averaged model based on the momentum closure and the zonally averaged density field found in observations. Our results suggest that this may be because of the role played by Southern Hemisphere wind forcing in the dynamics of the overturning circulation.

In the introduction, mention was also made of the vorticity-based closure for zonally averaged models introduced by Wright et al. (1995). This closure is much more akin to the Stommel–Arons approach than the Stommel box model, and is also in much better agreement with coarse-resolution ocean models than the momentum-based closure (see Fig. 3 in Wright et al. 1998). Although the vorticity based closure, like the momentum approach, is strictly inconsistent with the Stommel–Arons solution, the agreement with coarse-resolution ocean models is impressive, and it seems likely to be a valid closure for zonally averaged models over a much broader range of values of *R* than the momentum-based closure. The weakness of the vorticity based closure is the boundary condition that must be imposed at the boundary extremities. Wright et al. (1998) provide some justification for the particular choice made by Wright et al. (1995), but even they admit that the form of the boundary condition is “somewhat uncertain.” Nevertheless, the results presented here, and the work of Wright et al. (1998), suggest that the vorticity-based closure has wider applicability than the momentum-based closure and is to be preferred. Furthermore, the vorticity approach (or a variant thereof) may be able to distinguish the direction of the overturning circulation from the zonally averaged density field, even when quite close to the Stommel–Arons limit, since it is only in the Stommel–Arons limit itself (i.e., *R* → 0) that the relationship between the meridional overturning and the zonally averaged density field breaks down completely. Indeed, Wright et al. (1995) argue that the vorticity-based approach can reproduce the Stommel–Arons solution in an example where the location of the deep-water source is known.

Last, we air a note of caution regarding the Stommel box model approach. If the robustness of the box model approach relies on forcing by the Southern Hemisphere winds, then the applicability to different climate scenarios may be suspect, especially climate scenarios in which the effect of Southern Hemisphere winds is either eliminated or radically altered as compared with that in the present climate. We also acknowledge that the results here are not conclusive regarding the role of Southern Hemisphere wind forcing. We have used the simplest possible model [a shallow-water ocean model like that used by Kawase (1987)] to make our case, but clearly further work is required using more complete models. Of particular importance in future work on this topic is the need to control the diapycnal mixing in models in order to fully explore both the “strong” and “weak” damping regimes.

## Acknowledgments

This work has been supported by funding from NSERC, CICS, the Meteorological Service of Canada (MSC), and MARTEC, a Halifax-based company. We are grateful to Carsten Eden and Youyu Lu for their encouragement and for detailed comments on earlier versions of this work, and Lynne Talley for her comments and encouragement throughout the review process. Comments from Paola Cessi and two anonymous reviewers also led to considerable improvement in the manuscript. We are also grateful to Dirk Olbers and Richard Karsten for stimulating conversations on the dynamics of the Antarctic Circumpolar Current.

## REFERENCES

### APPENDIX A A Two-Density Layer Model

Consider a two-density layer ocean, as described in section 6.2 of Gill (1982). The upper (lower) layer has density *ρ*_{1} (*ρ*_{2}) and mean depth *H*_{1} (*H*_{2}), and the bottom is assume flat apart from unresolved irregularities that contribute to the topographic form stress. The equations governing the upper (with suffix 1) and lower (with suffix 2) layers and linearized about a state of rest are

Here *η* is the upward displacement of the free surface; *h* is the upward displacement of the interface between the two layers; *ν* is a vertical eddy viscosity coefficient associated with eddy form stress; *λh* represents diapycnal transport; *Q* is a source, for example, deep convection that transfers fluid from the upper to the lower layer (*Q* > 0); and *g*′ is reduced gravity: *τ*^{x}_{s} is the zonal component of the surface wind stress (the meridional component is set to zero), measured positive eastward, and *τ*^{x}_{t} is the zonal component of the topographic form stress, measured positive westward, and which acts on the lower layer (again we put the meridional component to zero).

Vertically averaging these equations, and neglecting the terms involving *g*′ (see Gill), leads to the equations for the barotropic mode; that is,

where (*u*_{b}, *υ*_{b}) is the vertically averaged horizontal velocity. Let us consider steady state dynamics (i.e., ∂/∂*t* = 0) for the reentrant channel geometry shown in Fig. 4. Denoting zonally averaged quantities by an overbar, it follows from (A9), using *υ*_{b} = 0 at the southern boundary, that *υ*_{b} = 0 everywhere; that is, there is no net meridional transport. Zonal averaging (A7) then shows that at the latitude of the channel

that is, the zonal wind stress is balanced by the topographic form stress, as we expect (Rintoul et al. 2001; Olbers and Ivchenko 2002).

We now follow section 6.3 in Gill and form the baroclinic mode. Denoting (*û, **υ̂*) = (*u*_{2} – *u*_{1}, *υ*_{2} – *υ*_{1}) and subtracting the momentum equations governing layer 1 from those governing layer 2, we obtain

where *K* = 2*ν*/(*H*_{1}*H*_{2}) is the Rayleigh friction coefficient appearing in (1) and (2). Equations (A3) and (A6) lead to the continuity equation for the baroclinic mode (neglecting the contribution from *η _{t}*)

It is equations (A11)–(A13) that correspond to equations (1)–(3) in section 2, with (*u, **υ*) in (1)–(3) replacing (*û, **υ̂*), *H* replacing *H*_{2}*H*_{1}/(*H*_{1} + *H*_{2}) and *τ*^{x}/(*ρ*_{o}*H*) replacing

Let us examine more closely the stress term

We note that the surface wind stress, *τ*^{x}_{s} and the topographic form stress *τ*^{x}_{t} both contribute to *τ*^{x} in (1). Further, since the surface wind stress is eastward (*τ*^{x}_{s} > 0) and the opposing topographic form stress is westward (*τ*^{x}_{t} > 0), with both roughly balancing each other in view of (A10), it follows that *τ*^{x} in (1) is negative, as we have assumed. Further, if *H*_{2} ≫ *H*_{1}, then the wind stress contribution dominates whereas, if *H*_{1} ≫ *H*_{2}, the topographic form stress contribution dominates. We note, however, that, since both the surface wind stress and the topographic form stress are expected to roughly balance each other, their relative contribution to *τ*^{x}/(*ρ*_{o}*H*) is not important here. The importance of the topographic form stress, emphasized by (A10), is a consequence of the rather special dynamics in the Southern Ocean where the zonal flow associated with the Antarctic Circumpolar Current extends to the bottom.

Last, we note that the treatment given here for a two density layer ocean can be extended to continuous stratification, following McCreary (1981), resulting in an equivalent expression for the stress term, provided the vertical eddy viscosity and the vertical diffusivity are inversely proportional to the buoyancy frequency, and using a slightly unconventional form for the vertical diffusivity term.

### APPENDIX B Scaling of the Western Boundary Layer

In steady state, it is straightforward to use (1) and (2) to express *u* and *υ* in terms of *h* to give (with *f* = *βy*)

Substituting these expressions for *u* and *υ* in (3) then gives the vorticity equation

The appearance of both the wind stress term and buoyancy forcing in (B3) emphasizes their equivalence as forcing mechanisms in the shallow-water context (cf. Gill 1982).

We note that for *K*^{2} ≪ *f*^{2} (in particular, away from the equator) and, since the scale of variations across the western boundary layer is much smaller than the scale of variations along the western boundary itself, (B3) can be reduced, near the western boundary and in the absence of local forcing, to

Balancing the *β* term with the second term on the right-hand side leads to the usual Stommel scaling for the boundary layer width, *L*_{w}; that is, *L*_{w} ≈ *K*/*β.* On the other hand, balancing the first two terms on the right-hand side gives a boundary layer width

where *c* = . In this case, the boundary layer width is the Rossby radius of deformation (as for inviscid Kelvin waves on an *f* plane) scaled by (*K*/*λ*)^{1/2}, as in Cane (1989), and as implied by the extension to midlatitudes of the analysis in Yamagata and Philander (1983).

What is of more interest here is the decay scale along the western boundary. We note that in Figs. 1c and 1d, there is no flow outside the boundary layer and that, in particular, *u* = 0 there. Since *u* = 0 at the western boundary itself, we can therefore put *u* = 0 inside the boundary layer. Equations (2) and (3) then give (in steady state with the forcing set to zero)

where *y* measures distance along the western boundary. We now nondimensionalize (B5) using the length scale *L,* which is the distance along the western boundary from the northern source to the equator, as in section 3. The nondimensional parameter to emerge is the product

where *R* = *λL*/*c* is the parameter introduced in Eq. (5). Hence, for fixed *K,* as in this paper, the important parameter is *R.* When *R* is large, the decay scale *c*/*λ* is less than *L,* and the influence from the northern source fails to reach the equator. On the other hand, when *R* ≈ 1, the Kelvin wave can reach the equator and the Stommel box model approximation starts to break down. For smaller values of *R* we eventually reach the Stommel–Arons regime.

## Footnotes

*Corresponding author address:* Dr. Richard J. Greatbatch, Department of Oceanography, Dalhousie University, Halifax, NS B3H 4J1, Canada. richard.greatbatch@dal.ca