## Abstract

Simple idealized layered models and primitive equation models show that the meridional gradient of the zonally averaged pressure has no direct relation with the meridional flow. This demonstrates a contradiction in an often-used parameterization in zonally averaged models. The failure of this parameterization reflects the inconsistency between the model of Stommel and Arons and the box model of Stommel, as previously pointed out by Straub.

A new closure is proposed. The ocean is divided in two dynamically different regimes: a narrow western boundary layer and an interior ocean; zonally averaged quantities over these regions are considered. In the averaged equations three unknowns appear: the interior zonal pressure difference Δ*p _{i}*, the zonal pressure difference Δ

*p*of the boundary layer, and the zonal velocity

_{b}*u*at the interface between the two regions. Here Δ

_{δ}*p*is parameterized using a frictionless vorticity balance, Δ

_{i}*p*by the difference of the mean pressure in the interior and western boundary, and

_{b}*u*by the mean zonal velocity of the western boundary layer.

_{δ}Zonally resolved models, a layer model, and a primitive equation model validate the new parameterization by comparing with the respective zonally averaged counterparts. It turns out that the zonally averaged models reproduce well the buoyancy distribution and the meridional flow in the zonally resolved model versions with respect to the mean and time changes.

## 1. Introduction

It is a common assumption in physical oceanography that the magnitude and sign of the zonally integrated meridional transport in the ocean [i.e., the meridional overturning circulation (MOC)] is related to the meridional pressure or density gradient. This assumption originates in the discussion of a two-box model by Stommel (1961), in which the exchange flow between the two boxes is parameterized with the density difference between the boxes. The physical basis of this closure is a hypothetical dynamical balance between the pressure difference induced by the different densities of the boxes and friction in a narrow pipe connecting the two parts of the ocean at depth.

A similar dynamical balance was also assumed by Marotzke et al. (1988) to close the momentum balance of the zonally averaged primitive equations. The Coriolis force is ignored, and a balance between the zonally averaged meridional pressure gradient and some kind of interior friction (Marotzke et al. 1988 choose vertical friction) acting on the meridional velocity *υ* is implemented in the meridional momentum balance, while momentum advection is assumed to be negligible. The last assumption is reasonable for scales larger than the internal Rossby radius. From this regime, a simple diagnostic relation

between the zonally averaged meridional transport and the meridional gradient of the zonally averaged pressure can readily be derived, where the positive parameter *γ* depends on the type for frictional parameterization (we will assume Rayleigh friction for simplicity but other forms are possible). Note that the wind stress forcing in Eq. (1) was ignored. It can be included in all closures discussed in the present study.

This relation for , together with the zonally averaged continuity equation to determine the vertical velocity , allows us to calculate the zonally averaged tracer balances. Here zonal velocity–tracer correlations, which introduce standing-eddy contributions in the tracer balances, are ignored. Wright and Stocker (1991) diagnosed the relation between and in a zonally resolved general circulation model and found indeed a positive value for the constant *γ*, which, however, depends on latitude. However, their particular choice of this relation is void of any dynamical fundament. Wright et al. (1998) give dynamical arguments to motivate a modified version of the closure, which leads to a relation very similar to Eq. (1) (see appendix A for details).

It is one purpose of this study to demonstrate that the closure from Eq. (1) is physically inconsistent. Although this point was already discussed by Straub (1996) and Greatbatch and Lu (2003), it was apparently not well received by the scientific community: there are currently several coupled earth system models of intermediate complexity with zonally averaged ocean model components relying on the closure given by Eq. (1) (Claussen et al. 2002). Because of their low computational costs, such models are often used for paleoclimate simulations and long-term climate projections—several of them are included in the Fourth Assessment Report of the Intergovernmental Panel on Climate Change (IPCC) (Solomon et al. 2007). Ocean-only versions are used for studies discussing the stability of the thermohaline circulation (e.g., Alexander and Monahan 2009). Furthermore, scalings for the global meridional circulation including the Southern Ocean and its impact on the circulation in zonally bounded basins still rely on Eq. (1) (Gnanadesikan 1999; Levermann and Fürst 2010). We would like to point out that the closure by Wright et al. (1995) is an exception; it does not rely on Eq. (1) as we will discuss in appendix A.

It is evident that the closures by Marotzke et al. (1988), Wright and Stocker (1991), and Wright et al. (1998) have in common that they call for a “downgradient” form of the meridional transport similar to what was assumed by Stommel (1961) for the viscous pipe flow in his two-box model, leading to a local relation between and . It was argued by Straub (1996) that this assumption is inconsistent with the model by Stommel and Arons (1960) describing the flow in a two-layer system. In that model, the zonal mean of the interface height between the layers, equivalent to the pressure in primitive equations, becomes independent of the location of the deepwater sources (i.e., independent of the sign and magnitude of the meridional transports), thus proving the closures based on Eq. (1) to be wrong. We call this contradiction between the two models by Stommel and Arons (1960) and Stommel (1961) “Straub’s dilemma” and further detail this point in the following section.

The models by Stommel (1961) and Stommel and Arons (1960) have different conceptual backgrounds and were developed to focus on different aspects of ocean dynamics. Therefore it cannot be a priori expected that both models are consistent with each other. Evidently, both models had success in describing important phenomena of the ocean dynamics. However, applying the strongly simplified assumptions of the Stommel (1961) model to zonally averaged models of Marotzke et al. (1988), Wright and Stocker (1991), and Wright et al. (1998), Straub’s dilemma cannot be ignored anymore because it reveals dynamical inconsistencies of these models.

The central purpose of the present study, however, is to present and validate an alternative closure for zonally averaged models, which generalizes the concept of Wright et al. (1995). Their closure is based on a meridional integration of the vorticity balance in the interior and in the western boundary layer. In their closure, the need emerges for an integration constant that is difficult to determine; but which sets the size and sign of the meridional transports. We also divide the ocean into an interior and a western boundary current, but instead of averaging the vorticity equation over these regions we work with the momentum and buoyancy (layer thickness) equations directly. Although this way we avoid the need to determine an integration constant, we still need parameterizations for the interior pressure difference of the boundary layer and the zonal velocity at the interface between these regions. A detailed comparison of two types of circulation models [a two-layer model (LM) like the one by Stommel and Arons (1960) and a general primitive equation model (PEM) with many levels] with their zonally averaged counterparts demonstrates the feasibility of the closure.

## 2. Straub’s dilemma

### a. Dilemma in a simple layered model

We first consider the model by Stommel and Arons (1960) in a slightly extended form also used by Greatbatch and Lu (2003), which is later referred to as LM. The governing equations for this model are given by

where *H* denotes the mean thickness of the lower layer of a two-layer ocean and its perturbation *h*, with a density difference *δρ*, between the two layers represented by the reduced gravity *g*′ = *gδρ/ρ*_{0}. The velocities *u* and *υ* are the differences between the upper- and lower-layer velocities. A prescribed deepwater source is denoted by *Q* and the interior upwelling is parameterized by the term −*λh* in the thickness balance of Eq. (4). The momentum balance in Eqs. (2) and (3) is taken to be linear, and friction induced by subgrid-scale processes is represented by Rayleigh friction with coefficient *r.* For a detailed derivation of the model equations see, for example, Gill (1982), their section 6.2, or Greatbatch and Lu (2003). There are two equations derived from Eqs. (2)–(4) that we present for later use. The momentum balance yields the vorticity balance

and using this equation to eliminate the divergence from the thickness balance we find

which is the potential vorticity balance. Implementing the geostrophic approximation of Eqs. (2) and (3) to eliminate *u* and *υ* turns this into the familiar form of the quasi-geostrophic vorticity equation

where *R* = *c*/*|f|* is the baroclinic Rossby radius and the Kelvin wave speed. The equation determines the long-term adjustment of the circulation by Rossby waves. It also determines the steady state.

The potential vorticity equation reveals the existence of a western boundary layer of the Stommel type with the familiar width *δ _{W}* =

*r*/

*β*, resulting from the dominant balance between advection of planetary vorticity and the torque by the Rayleigh friction:

*βυ*= −

*r*∂

*or . In the interior the planetary term*

_{x}υ*β*∂

*and the upwelling term −*

_{x}h*λh*/

*R*

^{2}dominate. Approaching the northern (or southern) rim of the domain, with

*h*→ const,

*υ*→ 0, and

*Q*≡ 0, upwelling and the Rayleigh friction term must balance in the steady state. This implies a meridional scale, , where

*δh*/

*h*is the relative variation of

*h.*These considerations can be used to construct an approximate analytical solution of Eq. (7). Here, however, a numerical model will be used.

For the experiments with LM we have used the parameter values *r* = 2 × 10^{−6} s^{−1}, *λ* = 1 × 10^{−9} s^{−1}, *g*′ *=* 0.02 m s^{−2}, *β* = 2.3 × 10^{−11} m^{−1} s^{−1}, and *H* = 400 m, which yields *c* = 2.8 m s^{−1}, *R =* 30 km (at *y =* 4000 km), *δ _{W} =* 100 km, and

*δ*

_{NS }

*=*300 km. For the latter,

*δh*/

*h*~ 0.1 is used. The system is integrated on an equatorial

*β*plane and the horizontal resolution is 20 km in the zonal and meridional directions. The zonal and meridional extent of the model domain is 2500 and 10 000 km, respectively. To demonstrate the influence of the transport and pressure field on the location of the deepwater source, we choose three different locations for

*Q.*The results of the three experiments are shown in Fig. 1. The location of the deepwater source is at the northwestern edge of the model domain for Fig. 1a, at the equator at the western boundary for Fig. 1b, and at the southwestern edge of the model domain for Fig. 1c. The lateral scale of

*Q*is

*δ*in both directions. In each experiment two dynamical different regimes exist: a narrow western boundary layer with a strong meridional flow and a weak interior flow whose meridional component is always poleward. The widths of the boundary layers at the western, northern, and southern rims confirm the above considerations.

_{W}In the interior the velocity field and the thickness contours are almost identical in all three cases. This is because the Sverdrup balance

obtained from Eq. (6) holds to a good approximation for steady conditions and *r*/(*βB*) ≪ 1, where *B* is the zonal width of the basin. We also need to know that *h* is related to *Q* only in an integral sense—that is,

in the integral over the whole model domain; however, the rhs of Eq. (9) has only contributions from the western boundary region. From Eq. (8) it is clear that the meridional interior transport is driven by the interior upwelling *λh*, which is almost identical for the three cases (i.e., almost uniform and of similar magnitude). The major differences occur, therefore, only in the boundary current of the individual experiments, which has to balance the interior flow and the upwelling in the interior and the different inflows of the localized deepwater source *Q*.

It is clear that the location and strength of the deepwater source *Q* in the thickness equation determines the total transport in the lower layer—that is, when *Q* is located at the northwestern corner, the total meridional transport is southward in both hemispheres of the domain (Fig. 1g, solid line) and it is anywhere northward for a deepwater source *Q* located at the southwestern corner of the domain (Fig. 1i), while the total transport is poleward in both hemispheres for an equatorial source (Fig. 1h). The source *Q* drives a total transport of about 5 Sv (1 Sv ≡ 10^{6} m^{3} s^{−1}) in the vicinity of the source in each case, which linearly reduces because of the interior upwelling into the upper layer with increasing distance from the source. Although the western boundary layer is much smaller than the total width of the basin, the transport in the western boundary layer is of similar magnitude to the total transport. It also has the same direction as the total transport, except for the region *y < −*2500, *|y|* > 2500, and *y >* 2500 km for the experiment with northern, equatorial, and southern sources, respectively, where it opposes the total transport.

Since the zonal integral of *h* is dominated by the interior, the zonally integrated *h* becomes independent of the location of the deepwater source *Q.* Consequently, *h* and in particular the meridional gradient of *h* become independent of the location of *Q*, and thus the sign and strength of the meridional transport is neither related to the zonal average of *h* nor its meridional gradient. This statement is in contrast to the box model by Stommel (1961) where the flow between the two boxes is parameterized by the meridional density difference between the boxes, and also in contrast to the closures by Marotzke et al. (1988), Wright and Stocker (1991), and Wright et al. (1998), which all depend on Eq. (1). To summarize, the parameterization in the box model by Stommel (1961) and the closures based on Eq. (1) are not consistent with the dynamics of the model by Stommel and Arons (1960). This inconsistency (Straub’s dilemma) was first noted by Straub (1996).

### b. The dilemma in primitive equations

The independence of the meridional gradient of the zonally averaged thickness from the meridional transport is not specific to layered models but is also found in a primitive equation model (Viebahn and Eden 2010; https://wiki.zmaw.de/ifm/TO/cpflame), referred to as PEM, with a similar configuration as the LM. In PEM we have neglected momentum advection (as before), and, for simplicity, the only tracer is temperature. The model domain is identical to the layered model, but there are 20 vertical levels of 50-m thickness, such that the domain is 1000 m deep. PEM is forced by relaxation of temperature in the uppermost grid box toward a target temperature, which is zonally and meridionally uniform except for a small region of meridional width *r*/*β* (equivalent to the western boundary-layer width) at the northern or equatorial region with a 3-K-smaller target temperature. This way, a northern or equatorial deepwater formation region is introduced as in the layered model. A case with a southern source is just a mirror of the one with a northern source and therefore not further discussed. The time scale of relaxation at the surface is 20 days. Convection in the case of unstable stratification is parameterized by setting the vertical diffusivity to very large values. As in LM, there is no wind forcing—that is, we focus here on the thermohaline circulation. Friction is identical to the LM, except that we introduce in addition lateral and vertical friction with viscosities of 3.2 × 10^{4} m^{2} s^{−1} and 1 × 10^{−3} m^{2} s^{−1}, respectively, since otherwise unphysical oscillations on a short time scale develop. We use the Quicker advection scheme (Leonard 1979) for tracers and vertical diffusivity of 1 × 10^{−4} m^{2} s^{−1} in addition.

The steady solution of PEM, shown in Fig. 2, indeed has much resemblance to LM. In the experiment with a northern source, there is a deep temperature minimum at the equator, and isopycnals below about 500-m depth are symmetric with respect to the equator, bending toward the bottom and toward the poles (Fig. 2b). A similar “hill,” symmetric around the equator, can be seen in the experiment with the equatorial source (Fig. 2d), although it is located more to the surface than at depth. Figures 2a,c also show the meridional transports for both experiments with PEM by the meridional streamfunction *ψ* with . The surface forcing drives a volume transport of a couple of Sv in both cases. In the case of the northern source, there is southward flow at depth, almost-uniform upwelling in the interior, and northward return flow at the surface. Sign, magnitude, and meridional structure of the meridional transport is also very similar to LM in the experiment with an equatorial source (see Figs. 2c,d).

As for LM, the meridional gradient of the zonally averaged density (i.e., temperature) or pressure (not shown) is similar in both experiments with PEM (cf. Fig. 1d with Fig. 2b and Fig. 1e with Fig. 2d for the case with northern and equatorial deep water sources, respectively) and is of opposite sign in both hemispheres. Their depth dependence differs. The meridional transport, on the other hand, does not show any direct dependency on , with respect to the individual hemispheres or experiments, proving the downgradient closures based on Eq. (1) to be wrong in primitive equation models as well. The reason is of course the same as in LM—notably the zonally averaged pressure is dominated by the interior zonal mean of *p*, which is in turn governed by the frictionless Sverdrup relation in the interior.

Greatbatch and Lu (2003) increased the vertical mixing to unrealistic values in LM and found that the zonally averaged thickness becomes dominated by the values in the western boundary layer, such that resembles more and more the meridional transports. We note here that introducing mesoscale eddy mixing results in a similar effect. However, again, unrealistically large values of the isopycnal thickness diffusivity are needed.

## 3. A consistent closure

We regard the closure for a zonally averaged model proposed by Wright et al. (1995) as dynamically consistent. It is based on a meridionally integrated, zonally averaged balance of vorticity (see appendix B). However, integrating the vorticity balance requires integration constants to be specified in each hemisphere. Wright et al. (1995) locate these constants at the northern and southern boundaries of the domain and relate them to the interior flow at the respective boundaries. The result is a nonlocal relation between the zonally averaged pressure (or thickness ) and the average meridional flow, as detailed in Eq. (B7) of appendix A. We have shown that the interior flow, and thus (or ), are nearly independent of the deepwater convection region and the direction and magnitude of the MOC. Since the choice of the integration constant determines the direction and magnitude of the MOC, it appears therefore problematic to use the interior flow for the choice. We therefore propose and evaluate a closure for zonally averaged models that avoids unknown integration constants. The closure is developed for the shallow water and the primitive equations.

### a. Closure for the layer model

Instead of considering the vorticity balance in the interior and the western boundary layer we simply use separate thickness and momentum balances averaged over these domains and keep also all time derivatives. Zonal averages of variables over the whole basin are indicated with overbars without an index, and zonal averages over the boundary layer or the interior carry an additional index *b* or *i*, respectively. The total meridional transport can be obtained by

and the zonally averaged height by

Here *B*, *B _{b}*, and

*B*denote the total basin width from the western to the eastern boundary, the width of the western boundary layer, and the width of the interior, respectively. Note that the boundary-layer width might be defined as a multiple of

_{i}*r*/

*β.*Obviously,

*B*≪

_{b}*B*

_{i}.#### 1) Motivation and evaluation of the closure

Averaging the system of Eqs. (2)–(4) over the western boundary at *x _{W}* to the offshore edge of the western boundary layer at

*x*, yields

_{δ}with Δ*h _{b}* =

*h*(

*x*=

*x*) −

_{δ}*h*(

*x*=

*x*)

_{W}*= h*In the thickness balance the zonal velocity

_{δ}− h_{W}.*u*=

_{δ}*u*(

*x*=

*x*) at the interface between the interior and the boundary layer appears. Furthermore, it was assumed in Eq. (14) that

_{δ}*u*(

*x*=

*x*) = 0 and the source

_{W}*Q*was located entirely in the western boundary layer. Likewise, the respectively averaged equations for the interior regime, extending from

*x*to

_{δ}*x*, are

_{E}with Δ*h _{i}* =

*h*(

*x*=

*x*) −

_{E}*h*(

*x*=

*x*)

_{δ}*= h*−

_{E}*h*Here

_{δ}.*u*(

*x*=

*x*) = 0 is used. To allow for the northern and southern boundary layers, described in section 3a, the friction terms have been retained though they are negligible in the actual interior. The pressure differences over the respective domains, Δ

_{E}*h*and Δ

_{b}*h*, as well as the zonal velocity

_{i}*u*, have to be parameterized. The resulting model will be referred to as the zonally averaged layer model (ZALM).

_{δ}We start by assuming that *u _{δ}* must be close to , hence we put

^{1}

in the thickness balances from Eqs. (14) and (17). Note that a linear increase of *u* within the western boundary layer would yield *γ*_{1} = 2. However, *u* is not increasing linearly over the western boundary layer (not shown) and we found that the best fit is obtained for *γ*_{1} = 1.7 [see also Figs. 3c,f where lhs and rhs of Eq. (18) are shown]. Next we demand that the thickness balance for the interior regime yields the averaged form of the Sverdrup balance from Eq. (8), assuming steady state and vanishing friction. To insert from the momentum balance of Eq. (15) into the thickness balance given by Eq. (17), we compute the meridional divergence of from the zonal momentum balance in Eq. (15), which becomes

We note that the choice

leads to an interior thickness budget in which only the *βυ* contribution, frictional, and tendency terms remain—that is, to an interior Sverdrup balance analogous to Eq. (8) [assuming that *βυ* dominates and in Eq. (19)]. Using Eq. (20) together with Eq. (18) as parameterization of Δ*h _{i}*, the integrated form of this relation reads

The integration constant at *y* = 0 follows from the steady state zonal balance at the equator:

Note that the Rayleigh model has the deficit that the equatorial field is completely decoupled from the rest. For this reason in Eq. (22) is replaced by the mean over three grid points across the equator. Note that the lhs and rhs of Eq. (21) are shown in Figs. 3b,e.

It remains to specify a parameterization for the pressure difference Δ*h _{b} = h_{δ}* −

*h*across the boundary layer. We have experimented with a variety of closures for Δ

_{W}*h*analogous to Eq. (20)—that is, motivated by the potential vorticity budget in the western boundary layer. However, many possible forms for the closure, which often yield an excellent fit to the respective variable in the zonally resolved model, turned out to lead to unstable numerical integrations. The simple ansatz

_{b}with another tuning parameter *γ*_{2} of order one, on the other hand, yields stable integrations in all cases, which we have considered, and is also a reasonable fit to Δ*h _{b}* from the zonally resolved model as discussed next [lhs and rhs of Eq. (23) are shown in Figs. 3a,d]. Further, the results of the integrations with the resulting zonally averaged model compare well with the zonally resolved counterparts, as discussed below, giving confidence to the parameterization from Eq. (23).

Figure 3 shows the western boundary thickness difference Δ*h _{b}* and its parameterization , both diagnosed from the experiments with the layered model LM, which is shown in Fig. 1. For the experiment with a northern (southern) source, the parameterization fits well with Δ

*h*except for the southernmost (northernmost) part of the domain in the experiment with a northern (southern) source, where Δ

_{b}*h*becomes negative, while stays positive. A similar deviation in the sign of the parameterization can be seen in the experiment with the equatorial source for large distances from the source. However, the structure of the meridional changes in are in all cases similar to Δ

_{b}*h*We also note that the quality of the parameterization depends on the exact definition of the width of the western boundary layer. Here, we have used

_{b}.*B*=

_{b}*r/β*with values for

*r*and

*β*as in the numerical experiments. The middle panels of Fig. 3 display Δ

*h*from the zonally resolved model and its parameterization from Eq. (21) which in fact agree very well. Figure 3 also shows the outflow from the western boundary region,

_{i}*u*, from the zonally resolved model together with its parameterization . Our choice of Eq. (18) fits

_{δ}*u*well with respect to the meridional structure and sign, while the magnitude is underestimated, which might be resolved by tuning the parameter

_{δ}*γ*

_{1}to values greater than one.

#### 2) Performance of the zonally averaged layer model

The complete ZALM consists of Eqs. (12)–(14) for the boundary layer and Eqs. (15)–(17) for the interior domain together with the parameterizations expressed in Eqs. (18), (21), and (23). ZALM was programmed in FORTRAN 90 and the source code together with a detailed documentation of all numerical details can be downloaded from https://wiki.zmaw.de/ifm/TO/zom. The steady state of a numerical integration of ZALM is shown as dashed lines in Fig. 4 and compared with the correspondingly averaged quantities diagnosed from LM. Note that the configuration of ZALM is identical to LM in all respects (except for the zonal extent and the closure). Transports and thickness height are reasonably well reproduced by ZALM for the northern and equatorial sinking case using *γ*_{1} = *γ*_{2} = 1 and a boundary-layer width of *B _{b}* =

*r*/

*β.*However, by changing the tuning parameters to

*γ*

_{1}= 1.7,

*γ*

_{2}= 1.2, and

*B*= 2

_{b}*r*/

*β*, the broad central hill in and the structure at the northern and southern boundaries are even better reproduced (dotted lines in Fig. 4). This improvement was found after some educated trials. Probably an even better improvement could be reached by using a parameter optimization procedure but this is not focus of this study.

We next take a closer look at the physical processes that establish the circulation in ZALM. Kawase (1987) showed that the establishment of the deepwater circulation involves basin-wide propagating Kelvin and Rossby waves: a thickness anomaly generated at the northern boundary of the basin propagates along the western boundary southward in the form of a Kelvin wave; at the equator, the Kelvin wave turns into an equatorial Kelvin wave and crosses the basin toward the east where it is again reflected and propagates at the eastern boundary north- and southward; westward propagating long Rossby waves, emanating from the eastern boundary, then transfer the signal into the interior of the ocean. These processes are of course realized in LM and a very similar adjustment process can be found in ZALM.

It is easily confirmed that the dynamics in the western boundary layer of ZALM allow for a Kelvin wave. With and vanishing friction, diffusion, and forcing, the equations become

which yields the familiar wave speed of . The zonal velocity is in geostrophic balance. In a corresponding way, Kelvin waves exist for the interior regime of ZALM, which are, however, attached to the eastern boundary in the zonally resolved model. Both regimes also support meridionally propagating gravity waves, which are coupled via the pressure terms and the *u _{δ}* term in the thickness balances.

Because of the zonal averaging, equatorial waves and midlatitude Rossby waves appear in a quite hidden way in ZALM. The Rossby wave response in midlatitudes is governed by the potential vorticity equations for the boundary and interior regime derived from Eqs. (12)–(14) and Eqs. (15)–(17) together with the parameterizations expressed in Eqs. (18), (21), and (23). We find, omitting again friction, diffusion, and forcing,

with the Rossby radius . Here all tuning parameters of the parameterizations are set to one. Note that *τ _{α}* =

*B*/(

_{α}*βR*

^{2}) is the time that a baroclinic Rossby wave needs to cross the respective region

*α = i*,

*b.*The corresponding time scale for the interior is roughly 10 yr in the northern part; it decreases toward the equator to several days (7 days at

*y*= 300 km). For the boundary layer, the time scale is considerably smaller (by the factor

*B*/

_{b}*B*≃ 25). The Rossby wave communication between the two regimes is thus represented by an oscillation of the mean layer thicknesses and with a period proportional to . In addition there is a meridional propagation of the perturbation with the correct Rossby wave speed, expressed by the meridional derivative terms in the tendency terms.

_{i}To assess the temporal behavior, Fig. 5 shows and for ZALM in comparison to LM from the start of both integrations. In the initial phase of both simulations, the anomaly in the interface height produced by the deepwater source is distributed via a Kelvin wave propagating from the northern edge of the model domain along the western boundary toward the south. During this stage the interior is still quiet. The Kelvin wave reaches the equator after approximately 20 days (see Figs. 5a,c) in both ZALM and LM. We note that in ZALM the propagation speed of the wave response at the western boundary depends to some extent on the western boundary width *B _{b}.* Increasing (decreasing)

*B*from

_{b}*r*/

*β*—the value used in Fig. 5—to larger (smaller) values, the southward propagation speed decreases (increases) slightly. The reason is the increasing (decreasing) importance of for the dynamics, which should be zero for a pure Kelvin wave, but which is present in both numerical integrations, reducing the southward propagation (Kelvin wave) speed.

The signal is then transported from the boundary layer into the interior (see Figs. 5b,d). This is first achieved by the increasing imbalance of the Kelvin wave dynamics: approaching the equator, the geostrophic relation for cannot be sustained, and a zonal velocity must increasingly develop. This disturbance thus couples into the interior thickness balance and the resulting thickness perturbation spreads to the north and to the south, involving the northward propagating Kelvin wave response at the eastern boundary and zonally and meridionally propagating Rossby waves in the interior. The time scale of this subsequent adjustment is also very similar in LM and ZALM.

### b. Application to primitive equations

The application of the closure, discussed so far for the layered model, to primitive equations is straightforward. Averaged separately over the western boundary layer and over the interior, the equations become

with *α = b*, *i* indicating the boundary or interior part, respectively, and *ε _{b}* = 1 and

*ε*= −1. The (scaled) pressure is related to the buoyancy by the hydrostatic relation . Note that standing-eddy fluxes are neglected in Eq. (29). Friction is contained in and and is specified below. Momentum advection has been neglected as before for the layered model. Convection is parameterized by using large values of the vertical diffusivity

_{i}*K*in case of unstable stratification.

_{α}The pressure differences over the western boundary layer and the interior, Δ*p _{b}* =

*p*(

*x*=

*x*) −

_{δ}*p*(

*x*=

*x*) and Δ

_{W}*p*=

_{i}*p*(

*x*=

*x*) −

_{E}*p*(

*x*=

*x*), respectively, and the zonal velocity

_{δ}*u*at the offshore edge of the western boundary need parameterizations. Analogous to the closure in the layered model, we use

_{δ}The interior pressure difference at the equator, Δ*p _{i}*(

*y*= 0), is again set by the steady zonal momentum balance at the equator. The model includes a rigid lid surface boundary condition and a diagnostic relation to find the surface pressure as usual in ocean general circulation models. This will be called the zonally averaged primitive equation model (ZAPEM).

Figure 6 shows Δ*p _{b}*, Δ

*p*, and

_{i}*u*diagnosed in PEM and their parameterizations given by Eq. (31). There is a good agreement concerning sign and structure of the variables and their parameterizations. Only in the southernmost part does the parameterization for Δ

_{δ}*p*not show the correct sign, which is similar to what we have seen for LM (see Fig. 3a). It turns out that an important parameter is the boundary-layer width

_{b}*B*, which we have chosen here as

_{b}*B*= 3.7

_{b}*r*/

*β*since this value seems to match best the boundary-layer width in PEM. Note that the boundary layer is broader than expected from the Rayleigh friction term because we also have included harmonic friction in PEM, which leads to a wider boundary layer.

ZAPEM was programmed in FORTRAN 90 and the source code as well as a documentation of all important numerical details can be downloaded from https://wiki.zmaw.de/ifm/TO/zom. The results of ZAPEM after 200 yr of integration for two different surface boundary conditions are shown in Fig. 7.

The configuration and relevant parameters of ZAPEM are identical to the zonally resolved model version (PEM) shown in Fig. 2, except for the zonal extent, the closure, and that we have omitted the harmonic zonal friction terms (meridional friction is kept) in ZAPEM. Two different surface boundary conditions are implemented using the two different target surface temperatures, as described in section 2b. Without parameter optimization (i.e., for *γ*_{1} = *γ*_{2} = 1 and *B _{b}* = 3.7

*r*/

*β*) we found already good agreement between ZAPEM and PEM; we therefore made no further attempt of parameter tuning. However, for the polar sinking case (see Figs. 7a,b) the overturning rate in ZAPEM is slightly too strong while the vertical stratification is slightly too weak. For the equatorial sinking (see Figs. 7c,d) the reverse statement holds.

For comparison we also present a primitive equation simulation with the inconsistent closure of the form from Eq. (1). We use the zonally averaged meridional momentum equation

which is the time-dependent case of Eq. (1) and similar to the closure proposed by Marotzke et al. (1988) and Wright and Stocker (1991). Note that the Coriolis term is omitted in the meridional momentum balance in Eq. (32) and replaced by a large Rayleigh damping term. Note also that horizontal and vertical friction is included here only for a consistent comparison with the other simulations; there is no qualitative difference in the results with and without these terms (not shown). The zonally averaged meridional velocity from Eq. (32) is used in the budget for for which no further closure is needed; is calculated from the continuity equation.

Figure 8 shows results of an integration using Eq. (32) as closure with *γ*_{WS} = 40. While the overturning circulation is similar to ZAPEM and PEM, the structure of the buoyancy field reveals a major disagreement for the case of a northern deepwater source (Fig. 8b). According to Eq. (1), the sign of the meridional buoyancy (pressure) gradient cannot change if the streamfunction consists of one single overturning cell. This clearly contradicts the results of PEM and ZAPEM. Only for the case of the equatorial deepwater source, the buoyancy distribution (Fig. 8d) conforms better with that of ZAPEM and PEM, although the northern and southern boundary layers as observed in ZAPEM and PEM are not reproduced.

Figure 9 illustrates the effect of wind forcing and a Southern Ocean part in simulations with ZAPEM and PEM. It is straightforward to include wind forcing and/or zonally periodic boundary conditions in the zonally averaged models ZALM and ZAPEM: the zonally averaged wind stress is used as upper-boundary condition in the vertical stress divergences contained in and . For the zonally unbounded periodic part of the domain, as found in the Southern Ocean, the zonal pressure differences are simply set to zero. Their dynamical role is replaced by the effect of mesoscale eddies; see, for example, Olbers and Visbeck (2005). This process can be included by interpreting the momentum balance as a balance for the residual velocity—that is, the sum of the Eulerian mean velocity and the eddy-driven (bolus) velocity (Andrews et al. 1987; Ferreira and Marshall 2006; Viebahn and Eden 2010). The effect of mesoscale eddy density mixing is then represented by vertical friction with viscosity *K*_{gm}*f*^{2}/*N*^{2} where *K*_{gm} denotes the isopycnal thickness diffusivity according to the Gent and McWilliams (1990) parameterization. We simply take a constant value of *K*_{gm} = 1000 m^{2} s^{−1}.

For the simulations shown in Fig. 9, we have chosen the same model domain as before, but for the southern quarter of the basin we apply zonally periodic boundary conditions to represent the Southern Ocean. Note that the setup is similar to that used in Viebahn and Eden (2010): the wind stress over the Southern Ocean region is zonally constant and sinusoidal in the meridional coordinate with a maximum of 2 × 10^{−4} m^{2} s^{−2} located at the center of the periodic domain. The wind stress in the zonally bounded part of the domain is set to zero. The surface boundary condition for buoyancy is a relaxation toward a target buoyancy restoring function with a linear increase (with a rate of 10^{−9} s^{−2}) in the Southern Ocean region, a constant value from *y* = −2560 to *y* = 2560 km, and a linear decrease of the target buoyancy at the northern part (with the same rate), which generates an equivalent temperature difference of about 20 K between the equator and the polar boundaries. As before, the zonally averaged model (ZAPEM) with our new closure is compared to a simulation with PEM in an identical (but zonally resolved) setup (Fig. 9). ZAPEM again reproduces well PEM, although ZAPEM again slightly underestimates the overturning and overestimates the stratification in comparison with the corresponding PEM experiment. We further note that temperature and salinity and further passive tracers can be added as variables to the zonally averaged model (not shown). Isopycnal mixing is also implemented as an additional mixing term in the tracer balances. It is also straightforward to include variations in the ocean depth.

## 4. Summary and discussion

The box model by Stommel (1961) and the model by Stommel and Arons (1960) both aim to describe the meridional overturning circulation of the ocean. Much of our knowledge about this important aspect of the ocean’s circulation is based on these models. However, Straub (1996) pointed out an inconsistency between the box model and the Stommel and Arons (1960) model, which proves the assumption in Eq. (1) to be inconsistent, which is an inherent assumption in the box model of Stommel (1961) and also in many zonally averaged ocean models (Claussen et al. 2002). We call this inconsistency Straub’s dilemma, representing the fact that it appears not possible to infer the meridional transport from the meridional gradient of the zonally averaged pressure. This is because the zonally averaged pressure is dominated by the interior pressure, which, on the other hand, is governed by frictionless and linear dynamics expressed by the Sverdrup relation of Eq. (8). Since this vorticity balance is driven only by the interior upwelling, it is unrelated to the sign of the meridional flow.

In this study, we present and evaluate a new and consistent closure for zonally averaged models to replace the inconsistent closure given by Eq. (1), illustrated by numerical integrations with a layered model version and a version based on the full primitive equations. Following Wright et al. (1995), the model domain is divided into an interior part—governed by the Sverdrup relation of Eq. (8)—and a boundary layer part, where friction plays an important role in the vorticity balance. However in contrast to Wright et al. (1995), we do not use the vorticity balances of the interior and the boundary layer directly, but instead use the zonally averaged, interior and boundary layer, momentum and thickness (or buoyancy) budgets. The reason for doing so is that using the meridionally integrated vorticity balances, as suggested by Wright et al. (1995), introduces the need to specify an unknown integration constant. We find the choice of this integration constant to be problematic, since it sets the sign of the meridional transport.

Therefore, we use the vorticity balance only to motivate the parameterization of the zonal pressure difference over the interior, which is needed for the zonally averaged interior zonal momentum balance. The zonal pressure difference across the boundary layer, on the other hand, is parameterized by the difference of the zonally averaged pressure in the interior and the pressure averaged over the boundary layer. The advective exchange between the boundary layer and interior is parameterized using the mean zonal velocity in the boundary layer. The standing-eddy fluxes in the nonlinear buoyancy budgets are simply neglected. Both in the layered model and the primitive equation model we find good agreement with respect to the evaluation of the parameterizations and model results in terms of the mean simulation of the transports and the thickness (buoyancy) and its time changes.

We advocate replacing the inconsistent closure of Eq. (1) with the new closure discussed in this study in zonally averaged ocean models. However, we do not imply that the box model by Stommel (1961) is inconsistent as well. On the one hand, the interpretation of the meridional flow in the ocean as driven by the pressure difference between two boxes and controlled by friction in a hypothetical pipe connecting the two boxes by Stommel (1961) is certainly an incorrect oversimplification of the real dynamics. On the other hand, many results and predictions of the box models can be reproduced by models, including the correct dynamics. This agreement might give some confidence in the box model, although we know that its dynamics are incomplete and only a very rough analog to the real dynamics. We hope that the new zonally averaged model presented here can contribute to further confirm and extend knowledge from the box model about the meridional flow in the ocean.

## Acknowledgments

This work was supported Grant by BMBF-SOPRAN FKZ 3F0611A.

### APPENDIX A

#### Some Frequently Used Inconsistent Closures

In this section and the next, we discuss some closures analogous to Eq. (1) using the layer equations for simplicity, but the results easily transfer to primitive equations. We also neglect wind forcing, which can, however, easily be incorporated. Marotzke et al. (1988) proposed a closure by abandoning the Coriolis force and implementing (unrealistically) large friction into the meridional momentum balance:

which leads directly to Eq. (1) with *γ* = 1/*r*. Here *g*′ is the reduced gravity, *h* the layer thickness, *υ* the meridional velocity component, and and their zonal averages. Note that we use here Rayleigh friction with friction coefficient *r* to connect to the model by Stommel and Arons (1960), while Marotzke et al. (1988) originally used vertical diffusion of momentum. However, the specific choice of the friction does not change the fundamental relation in Eq. (1).

A similar relation was proposed by Wright and Stocker (1991). They consider the zonal momentum balance in the zonally averaged form, where the east–west pressure difference *h _{E}* −

*h*over the basin width

_{W}*B*needs a parameterization. They choose

where *φ* denotes latitude and γ_{WS} a constant of order 1. The zonal pressure difference is thus expressed in terms of the local meridional pressure gradient, which is also of the form in Eq. (1), assuming that the meridional flow is in geostrophic balance, with a parameter *γ* proportional to cos*φ.* This setting is supported by numerical experiments with a three-dimensional (but highly simplified) circulation model, but as argued by Greatbatch and Lu (2003), the support is due to the highly diffusive nature of the ocean model. Note that the closure is not based on any dynamical concepts.

Wright et al. (1998) avoid a direct closure for the pressure difference *h _{E} − h_{W}.* The zonal momentum balance is entirely abandoned. In fact, the zonally averaged meridional momentum balance is written as

with Coriolis parameter *f*, the zonally averaged zonal velocity , and its geostrophic component . To determine the meridional velocity from Eq. (A3), the ageostrophic zonal velocity must be known and thus needs to be parameterized. For this reason, Wright et al. (1998) divide the zonal extent *B* of the ocean again into a western frictional boundary layer part of width *B _{b}* and an interior part of width

*B*≫

_{i}= B − B_{b}*B*. They write

_{b}with the zonal velocities and averaged over the interior and western boundary layer, respectively, and where the superscript (*g*) denotes the geostrophic component of the velocity. The interior flow is largely geostrophic and, thus, Wright et al. (1998) assumed the product to be small. In the boundary layer, the flow has both a geostrophic and an ageostrophic component, but *u* vanishes on the continental side of the layer and should be largely governed by the geostrophic balance on the offshore edge of the western boundary layer. The interior geostrophic component continues only moderately changed into the boundary layer and to the actual boundary. Hence, the magnitudes of should be similar but of opposite signs in the boundary layer, and

should hold. Inserting the parameterized ageostrophic velocity into the meridional momentum balance then yields

which is identical to Eq. (1) with a suitable parameter *γ* = *B _{b}*/(

*rB*). Note that this closure is entirely of geometric nature: it uses the observed structure of a basin-wide circulation with a narrow western boundary current but not any further dynamics.

### APPENDIX B

#### The Consistent Closure by Wright et al.

Wright et al. (1995) propose a dynamically consistent closure by splitting the ocean basin into a western boundary layer and an interior and considered the vorticity budgets averaged separately over both regions. Assuming a frictionless interior and using specific parameterizations for friction in the western boundary layer, they derive a nonlocal relation between the meridional transport and zonally averaged pressure. Because the original closure of Wright et al. (1995) needs the specification of integration constants that are difficult to determine, we have presented in section 3a generalization of the concept by Wright et al. (1995), which does not need the specification of integration constants.

We found the derivation of the closure in Wright et al. (1995) unnecessarily complicated. Here we give an alternative simplified derivation with fewer assumptions to arrive at a similar equation. The analysis is again performed for the layer model and starts with the zonal momentum balance in the zonally averaged form

which are identical to Eqs. (12) and (15), neglecting the time tendency terms. Again, indices *W*, *E*, and *δ* denote that the values are taken at the western or eastern boundary or at the interface between interior ocean and boundary layer, respectively. The overbars denote zonal averages over the interior with additional index *i* or boundary layer with index *b,* and *B _{b}* and

*B*denote the width of the boundary layer and the interior, respectively. Only a few approximations now lead to the closure by Wright et al. (1995). First, the friction term in Eq. (B2) will be neglected. Because of the kinematic boundary condition at the eastern boundary,

_{i}*u*= 0, it follows from Eq. (3) that ∂

_{E}*= 0 or*

_{y}h_{E}*h*= const. Second, the thickness perturbation

_{E}*h*along the western boundary in Eq. (B2) is eliminated by the meridional velocity

_{W}*υ*using the steady meridional momentum balance at

_{W}*x*in the form

_{W}and *υ _{W}* is parameterized by . Note that

*u*= 0 was assumed. Next the friction coefficient

_{W}*r*in Eq. (B3) is replaced by

*βB*using

_{b}*B*=

_{b}*r*/

*β*as the boundary-layer width according to Stommel (1948). The meridional integral of Eq. (B3) with starting point at

*y*

_{0}can be used to eliminate

*h*from Eq. (B1) to end up with

_{W}where integration limit *y*_{0} is arbitrary. The meridional velocities and then follow from

and are seen to be both determined by *h _{δ}. *Wright et al. (1995) propose the closure where denotes the zonally averaged thickness, and neglect the last term in Eq. (B5) related to friction. Both are quite good assumptions for the layered model outside the northern and southern boundary layers (not shown). Note, however, that Eq. (B5) only determines the derivative of and thus it it is necessary to set an integration constant for . One may take , which by Eq. (B5) is obviously related to the unknown

*h*(

_{W}*y*

_{0}). Note that the frictionless interior balance leads to .

To arrive at the central equation of the Wright et al. (1995) model, Eq. (B5) is divided by *f* and integrated from *y*_{0} to *y* (in the same hemisphere to avoid the singularity at *y* = 0) to give , involving now the unknown . The total meridional flow is then governed by

Wright et al. (1995) use as integration constant the boundary transport at the northern and southern boundary, which they relate to the interior flow at the respective boundary [Eq. (B7) is used twice to circumvent the singularity at the equator]. It becomes clear that this closure implies that the information about the placement of the deepwater source—which is invisible to the interior flow—must be contained in at the northern and southern boundary layers.

Figures 1d–f show indeed that at the northern (southern) boundary layer for the experiment with northern (southern) sinking is slightly higher and reaches a larger value at the northern (southern) end of the domain than in the experiment with equatorial and southern (northern) sinking. It is this small difference that has to determine the sign of the flow in the nonlocal relation between and of Wright et al. (1995). Consequently, an evaluation (not shown) of the closure based on Eq. (B7) in the layered model shows that it is not able to predict using only from the model. The reason is that the assumption leading to the closure (i.e., a frictionless interior flow) breaks down in the northern and southern boundary layer. The unknown integration constant then determines the meridional flow. We propose in section 3a a more robust way to determine the meridional flow, which avoids Eq. (B3) and thus the meridional integration and appearance of unknown integration constants.

## REFERENCES

## Footnotes

^{1}

In this paper we denote by the sign that we introduce a parameterization.