A simple and self-consistent geostrophic-force-balance model of the thermohaline circulation with boundary mixing

A simple model of the thermohaline circulation (THC) is formulated, with the objective to represent explicitly the geostrophic force balance of the basinwide THC. The model comprises advective-diffusive density balances in two meridional-vertical planes located at the eastern and the western walls of a hemispheric sector basin. Boundary mixing constrains vertical motion to lateral boundary layers along these walls. Interior, along-boundary, and zonally integrated meridional flows are in thermal-wind balance. Rossby waves and the absence of interior mixing render isopycnals zonally flat except near the western boundary, constraining meridional flow to the western boundary layer. The model is forced by a prescribed meridional surface density profile. This two-plane model reproduces both steady-state density and steady-state THC structures of a primitive-equation model. The solution shows narrow deep sinking at the eastern high latitudes, distributed upwelling at both boundaries, and a western boundary current with poleward surface and equatorward deep flow. The overturning strength has a 2/3-power-law dependence on vertical diffusivity and a 1/3power-law dependence on the imposed meridional surface density difference. Convective mixing plays an essential role in the two-plane model, ensuring that deep sinking is located at high latitudes. This role of convective mixing is consistent with that in three-dimensional models and marks a sharp contrast with previous two-dimensional models. Overall, the two-plane model reproduces crucial features of the THC as simulated in simple-geometry threedimensional models. At the same time, the model selfconsistently makes quantitative a conceptual picture of the three-dimensional THC that hitherto has been expressed either purely qualitatively or not self-consistently.


Introduction
The oceanic thermohaline circulation (THC), a key agent in Earth's climate, has as yet denied researchers important insights into its dynamics.Attempts to capture THC dynamics in a reduced-complexity model (Marotzke et al., 1988;Wright and Stocker, 1991;Sakai and Peltier, 1995;Wright et al., 1995) have long disregarded at least one of two fundamental observations: that diapycnal mixing, one of the driving mechanism of the meridional overturning, is essentially confined to lateral boundaries (Munk, 1966;Wunsch and Ferrari, 2004) and that the western boundary current and, by implication, the zonally integrated thermohaline overturning are in geostrophic balance (Johns et al., 2005).Marotzke (1997) presented an analytical THC theory that adheres to these principles but must make several ad-hoc assumptions.The current paper extends the Marotzke (1997) theory such that THC dynamics are represented in reduced yet self-consistent form.
Our starting point is the assumption that upwelling from the abyss and, by implication, total THC strength are limited by the rate of diapycnal mixing (e.g.Munk, 1966;Wunsch and Ferrari, 2004).Munk (1966) already conjectured that diapycnal mixing might be concentrated near the ocean boundaries.Enhanced boundary mixing would make possible a realistic THC strength while accounting for the observed low diapycnal mixing rates in the ocean interior (e.g.Armi, 1978;Ledwell et al., 1993;Wunsch and Ferrari, 2004).However, boundary mixing was not implemented into a primitiveequation model until Marotzke (1997), who found flow patterns very similar to the uniform-mixing case, except that vertical motion was entirely confined to lateral boundary regions and that interior flow was more zonal (see also Scott, 2000, who confirmed the results using improved resolution and a more realistic isopycnal mixing scheme).
An alternative explanation of the meridional overturning strength that does not require vigorous diapycnal mixing is that wind forcing over the Southern Ocean drives the overturning (Toggweiler and Samuels, 1995; for theoretical models implementing this suggestion, see Gnanadesikan, 1999;Johnson et al., 2007).Here, however, we follow Munk's (1966) conjecture and work within the boundary-mixing framework.
The quantitative dependence of the THC strength on diapycnal mixing has been confirmed repeatedly to show a 2/3-power-law (e.g.Marotzke, 1997;Zhang et al., 1999;Park and Bryan, 2000).This power law complies with a simple scaling law derived by assuming a unique vertical scale for thermal-wind balance, the continuity equation, and advective-diffusive balance (Welander, 1971;Bryan, 1987;Colin de Verdière, 1988).This scaling requires the zonal density difference to scale with the imposed surface meridional density difference, which Park and Bryan (2000) showed to hold true in the parameter range considered, but which hitherto has remained unexplained.
Another attempt to gain insight into THC dynamics was the formulation of two-dimensional models.Some of these (Marotzke et al., 1988;Wright and Stocker, 1991;Sakai and Peltier, 1995) are based on the zonally averaged momentum equations and a closure assumption relating the meridional flow to the meridional density gradient (Wright et al., 1998).The model formulation given by Wright et al. (1995) draws on vorticity dynamics; the closure is based on vorticity dissipation in a western boundary layer.All these two-dimensional models yield density and velocity fields in qualitative agreement with the zonal averages of primitiveequation model solutions.The parameters, however, must be chosen by calibration to these primitive-equation model solutions.
The closure assumptions relating the meridional flow to the meridional density gradient (Marotzke et al., 1988;Wright and Stocker, 1991;Sakai and Peltier, 1995) rely on a substantial viscous deviation from geostrophy in the western boundary layer (Wright et al., 1998), although western boundary currents have been observed to be geostrophically balanced (Johns et al., 2005).The vorticity-based closure (Wright et al., 1995) is in principle consistent with a geostrophically balanced western boundary current.How this balance is affected by their approximations, however, remains unclear.
Another deficiency of these two-dimensional models is that convective mixing appears to be unimportant for THC dynamics (Marotzke et al., 1988).When convective mixing is switched off, statically unstable stratification occurs, but the gross density and overturning structures only change marginally.This marks a sharp contrast with threedimensional models that show a reverse cell at high latitudes and a substantially lower abyssal density when convective mixing is switched off (Zhang et al., 1992;Marotzke and Scott, 1999).
Due to their zonally averaged nature, the two-dimensional models are also incompatible with boundary mixing.When the temperature and salinity equations are averaged zonally, vertical eddy diffusivities are assumed to be constant in all these studies (Marotzke et al., 1988;Wright and Stocker, 1991;Sakai and Peltier, 1995;Wright et al., 1995).Marotzke (1997) proposed an analytical THC theory that is based on boundary mixing.Motivated by his primitiveequation model results, he prescribed density structures at both eastern and western boundaries, allowing him to deduce the overturning from vertical advective-diffusive and thermal-wind balances.Bulk features like the mass and heat transports' latitudinal dependence are well captured, and the THC strength's scaling can be derived.But the conceptual disadvantage of his theory is that density structures are prescribed at both boundaries -only a depth scale, the pycnocline depth, can be determined as part of the solution.Because density structures are crucial in setting both shape and strength of the THC, this is a major drawback.Furthermore, Marotzke (1997) assumed level isopycnals at the eastern boundary, which deviates from those found in primitiveequation model solutions (Marotzke, 1997).
In the present study, we overcome these deficiencies and extend the Marotzke (1997) theory to a two-plane model such that self-consistent density and velocity fields can be deduced from a specified meridional surface density profile.We hence cast THC dynamics based on boundary mixing and geostrophy into a conceptual model capable of explaining salient features such as narrow high-latitude downwelling and distinct density structures at eastern and western boundaries.By formulating and solving the equations of our model, we make quantitative the purely qualitative conceptual picture put forward by Colin de Verdière (1993) and Zhang et al. (1992), which, in addition to boundary mixing, had been the conceptual starting point of Marotzke (1997).
In turn, we here make the Marotzke (1997) This paper is organized as follows: in Sect.2, the model formulation is given.Section 3 presents model solutions plus an analysis of parameter dependencies and an assessment of the role of convective mixing.Section 4 draws conclusions and summarizes the findings in a conceptual picture of the THC.In the Appendix, we elaborate on the computational procedure applied to solve the model equations.

The governing equations
The starting point of the model formulation is the geostrophically and hydrostatically balanced Boussinesq system (e.g.Vallis, 2006).The horizontal flow then satisfies thermalwind balance (1) Here, u denotes the horizontal part of the velocity field; b = −gρ /ρ 0 is the buoyancy, where g is the gravitational acceleration and ρ is the small deviation from the reference density ρ 0 ; f = 2 sinϑ is the Coriolis parameter as a function of latitude ϑ, where is Earth's angular frequency; k is the upward unit vector; z is the local vertical coordinate; and ∇ h is the horizontal part of the gradient operator.In this approximation, the flow is incompressible and hence satisfies the continuity equation where w denotes the vertical velocity.The steady-state buoyancy equation reads where q represents sources and sinks due to diapycnal processes.Following Sandström's theorem (e.g.Defant, 1961;Huang, 1999;Wunsch and Ferrari, 2004), diapycnal mixing is believed essential for maintaining an overturning circulation (see also Wunsch, 2005) and must be parameterized in these equations.For simplicity, we resort to vertical eddy diffusion.As buoyancy gradients are expected to be much smaller in the horizontal than in the vertical direction, the eddy diffusion's horizontal component is disregarded.However crude, this representation of diapycnal mixing proves useful at this stage of conceptualization as the role of mixing in our model is merely that of transporting buoyancy across isopycnals.Following the notion of boundary mixing (Marotzke, 1997), we only prescribe nonzero diapycnal mixing at lateral boundaries.
Another process that is essential for overturning dynamics is convective mixing (Zhang et al., 1992;Marotzke and Scott, 1999): as the water column becomes statically unstable, rapid mixing occurs, removing the instability and rendering the stratification neutral (Marshall and Schott, 1999).This process is parameterized here by a convective adjustment scheme, which is discussed in more detail in Appendix A1.Buoyancy sources and sinks due to convective mixing are represented by c, and the diapycnal term in Eq. (3) reads where k v is the vertical eddy diffusivity.The model geometry is for simplicity taken to be a hemispheric sector extending from latitude ϑ t to ϑ p and longitude λ w to λ e .The basin's depth d is taken to be constant.

Approximations
To satisfy a no-normal-flow condition at the eastern boundary and still permit buoyancy variations along the wall, the zonal velocity must turn ageostrophic near the boundary.The ensuing ageostrophic boundary layer circulation opposes the geostrophic flow such that the no-normal-flow condition is satisfied, but is essentially restricted to the zonal-vertical plane (Cessi and Wolfe, 2009).This is consistent with the numerical-model result that zonal isopycnal slopes are small in the eastern boundary layer (Marotzke, 1997;Cessi and Wolfe, 2009) and with the observation and numerical-model result that the meridional flow still adheres to the thermalwind balance (Johns et al., 2005;Cessi and Wolfe, 2009) and consequently essentially vanishes.
We parameterize this eastern boundary layer structure by assuming the buoyancy to be zonally constant and the zonal velocity to linearly tend to zero over the width of the eastern boundary layer λ (Fig. 1).The associated divergence is entirely compensated by vertical motion; the continuity Eq. (2) in the eastern boundary layer reads where u i is the interior zonal velocity, and a is Earth's radius.Note that this renders the eastern vertical velocity w e zonally constant across the boundary layer.The buoyancy Eq. (3) at the eastern wall, where not only meridional, but also zonal advection vanishes, simplifies to where b e is the buoyancy at the eastern wall.As discussed earlier, this boundary buoyancy determines the buoyancy throughout the boundary layer.The buoyancy balance inside the boundary layer is hence implicitly assumed between vertical advection and eddy fluxes (cf.Cessi and Wolfe, 2009;Cessi et al., 2010).In the interior of the basin, where both components of the flow are geostrophically balanced, the linear vorticity equation is satisfied: where β = 2 cosϑ/a.This equation, combined with the buoyancy Eq. ( 3) with zero vertical mixing, has a solution that consists of a zonally constant buoyancy field -implying zero vertical, zero meridional, and zonally constant zonal flow.Observational evidence indeed indicates that interior flow is essentially zonal (Davis, 1998;Hogg and Owens, 1999).To match the eastern boundary buoyancy, we therefore assume the buoyancy to equal its eastern boundary value throughout the interior.The thermal-wind balance for the interior zonal flow hence reads At the western boundary, we again apply semigeostrophy: the zonal velocity tends to zero linearly as over the eastern boundary layer, and the meridional flow satisfies thermalwind balance.Here, however, Rossby waves cannot flatten isopycnals zonally and a slope persists.We assume this slope to be zonally constant across the boundary layer (Fig. 1), such that the thermal-wind balance reads where b w is the buoyancy at the western wall and the meridional velocity v w is zonally constant across the boundary layer.
The continuity equation in the western boundary layer and the buoyancy equation at the western wall must include the meridional flow.They hence read, respectively, 1 a cosϑ Note that also the western vertical velocity w w is zonally constant across the boundary layer.Inside the western boundary layer, the buoyancy balance is hence implicitly assumed to be between meridional and vertical advection and eddy fluxes (potentially plus convective mixing).1 The six Eqs. ( 5), ( 6), and ( 8)-( 11) constitute a closed set for determining the quantities b w , b e , u i , v w , w w , and w e -all functions of latitude and depth only.We thus need to solve this set on a meridional-vertical plane extending from ϑ t to ϑ p and having constant depth d.
Note that this model formulation does in principle allow for meridional, but not for zonal variations of topography.That the latter can be neglected is supported by the finding that observations around the Mid-Atlantic Ridge are not essential for estimating the meridional overturning, which indicates that buoyancy differences across zonally varying topography are not of first-order importance for THC dynamics (Kanzow et al., 2007(Kanzow et al., , 2010)).

Boundary conditions
We now turn to boundary conditions, noting that the geostrophic approximation already eliminates the possibility of specifying a wind stress.As noted in the introduction and in line with Marotzke (1997), we focus on the buoyancy-driven circulation.
At the surface, buoyancy is prescribed at both eastern and western boundaries by a meridional profile, b w (ϑ,0) = b e (ϑ,0) = b 0 (ϑ) for all ϑ.The profile represents the observed meridional surface density distribution dominated by temperature (Fig. 2).Prescribing surface buoyancy does not accurately represent ocean surface processes.While fixed surface temperatures are a passable representation because surface heat fluxes strongly depend on these temperatures and anomalies are thus rapidly removed (Davis, 1976), surface salinity has essentially no effect on evaporation and precipitation.Surface salinity anomalies can hence persist on much longer time scales, and a flux boundary condition is generally preferable (Bryan, 1986).These mixed-type boundary conditions -prescribed salinity fluxes and fixed temperatures -have a crucial influence on THC stability and allow for multiple steady states (e.g.Stommel, 1961;Bryan, 1986;Marotzke et al., 1988;Marotzke and Willebrand, 1991).This type of analysis is left for further studies; we content ourselves with representing the salient THC dynamics, for which fixed surface buoyancy suffices.
At the basin's flat bottom, a no-flux condition is applied for buoyancy.We thus disregard geothermal heating because associated fluxes are relatively small and the resulting flows are weak (Huang, 1999;Adcroft et al., 2001;Scott et al., 2001).We write At the surface, we apply the rigid lid approximation; at the bottom, we impose a no-normal-flow condition.We hence require for all ϑ.The exclusion of wind forcing, together with the linear momentum balance, requires the vertically integrated horizontal flow to vanish.At both polar and tropical boundaries, we require the meridional flow to be zero: for all z.These conditions imply, by thermal-wind balance (Eq.9), that for all z.Also, by the respective buoyancy equation, easternand western-boundary vertical velocities must match at meridional boundaries if stratification is not neutral.This in turn implies that for all z.This can be seen by combining the thermal-wind relations ( 8) and ( 9) with the continuity Eqs. ( 5) and ( 10) and using the condition (16).Note that condition (17) necessitates ∂ ϑ b 0 (ϑ t ) = ∂ ϑ b 0 (ϑ p ) = 0 for consistency, which is fulfilled by the surface buoyancy profile (Eq.12).

The overturning stream function
The sum of continuity Eqs. ( 5) and (10) for, respectively, the eastern and western boundaries yields demonstrating that our approximations retain the zonally integrated flow's nondivergence.Convergence of vertical velocity thus always occurs in conjunction with divergence of meridional velocity at the western boundary, and vice versa.Furthermore, this nondivergence enables us to introduce a stream function ψ that characterizes the zonally integrated overturning in the meridional-vertical plane.We define ψup to an additive constant -by and The arbitrary constant is chosen such that ψ vanishes at the boundaries of the meridional-vertical domain.

Nondimensionalization
To gain insight into parameter dependence and to facilitate computational procedures, the Eqs.( 5), ( 6), and ( 8)-( 11) are nondimensionalized.Guided by the model equations, we define and nondimensionalize using these constants (see Table 1).Note that these do not necessarily represent appropriate scales (i.e. the resulting nondimensional quantities are not necessarily of order one); for instance, the vertical scale for buoyancy variations, the pycnocline depth, is expected to be much smaller than the basin depth d.
Using carets to denote nondimensional quantities, the thermal-wind relations (8) and ( 9) reduce to sinϑ the continuity Eqs. ( 5) and (10) to and the buoyancy Eqs. ( 6) and ( 11) to where the nondimensional vertical diffusivity  is the only parameter left in the system.A sensitivity study thus only needs to consider a one-dimensional parameter space.
The surface boundary condition on buoyancy (Eq.12) transforms to and the defining equations for ψ, Eqs. ( 19) and ( 20), become

Computational procedure
Since an analytical solution to the set of Eqs. ( 22)-( 27) does not appear to be known, we resort to numerical approximations.We discretize the equations on a 128 × 128 Arakawa C-grid, reinsert time tendencies into the buoyancy Eqs. ( 26) and ( 27), and numerically integrate (time step t = 10 −2 ), starting from an arbitrary initial guess (see Appendix A2), until a steady state is reached.To ensure numerical stability, we introduce artificial diffusion in the meridional direction.The artificial diffusivity κh is chosen small enough not to substantially affect the solution.The analysis described in Appendix A3 suggests that this is the case for κh = 2×10 −4 .The time stepping is performed using a time-splitting approach incorporating a second-order Runge-Kutta method for the advective terms, a fully implicit scheme for the artificial meridional diffusion, a Crank-Nicolson scheme for the vertical diffusion, and a convective adjustment scheme (see Appendix A1).

The standard solution
The external parameters used in this standard setup are given in Table 2; they approximately correspond to κv = 4.6 × 10 −5 .Here and in all subsequent solutions, the basin is taken to extend from ϑ t = 10 • to ϑ p = 70 • .The Equator is excluded because thermal-wind balances ( 22) and ( 23) can there not be applied in the present form.The steady-state velocity fields show very distinct patterns at the eastern and western boundaries (Fig. 3).At the western boundary, there is upwelling everywhere; it is strongest in the upper high-latitude ocean (Fig. 3c).At the eastern boundary, on the other hand, there is downwelling above upwelling equatorward of 61 • and strong downwelling reaching the bottom poleward of 61 • (Fig. 3d).Note that deep downwelling is much greater in magnitude than any upwelling.
Where downwelling exists above upwelling at the eastern boundary, the zonal flow exhibits a three-layer structure; at high latitudes, a two-layer structure prevails -consistent with the vertical divergence of eastern-boundary vertical velocities (Fig. 3a).
At the western boundary, the meridional current comes into play (Fig. 3b).Zonally and vertically convergent waters near the surface induce a strong poleward current intensifying with latitude.Near the polar boundary, zonal outflow balances meridional and vertical convergence in the upper ocean.At depth, the opposite is true near the polar boundary: zonal inflow feeds an equatorward current, supplying water to upwell or flow east at low latitudes.
The overturning stream function constitutes a single, hemisphere-wide overturning cell.It has a maximum of 15.1 × 10 6 m 3 s −1 that is located at 61 • and 1600 m depth (Fig. 4a).
The latitudinal dependence of the meridional overturning is captured by the stream function's vertical maximum (Fig. 4b).Fed by upwelling at both boundaries, the overturning steadily increases with latitude until it reaches its maximum at 61 • and sharply decreases further poleward.The maximum's proximity to the polar boundary reflects the narrowness of deep sinking.
The zonally and vertically integrated meridional buoyancy transport has its maximum at 38 • (Fig. 4c).
The velocity fields well capture primitive-equation model results.The solution presented in Marotzke (1997) equally exhibits upwelling at the western boundary as well as downwelling above upwelling and deep high-latitude sinking at the eastern boundary.The overturning stream functions show much resemblance both in structure and magnitude.That our model's overturning is slightly weaker than that found in Marotzke (1997) can at least partly be attributed to the basin being slightly smaller.
The steady-state buoyancy distributions exhibit a distinct pycnocline structure (Fig. 5): Throughout the abyss, buoyancy takes a nearly constant value, which is determined by the least buoyant surface water available.In the upper ocean, buoyancy drops from its surface value to the abyssal one within a few hundred meters.
In the western-boundary pycnocline, waters are well stratified -stratification is strongest near the surface -except for a region in the upper high-latitude ocean (Fig. 5a).There, the large advective buoyancy transport accomplished by the poleward surface current destabilizes the water column.As surface buoyancy is fixed, convective mixing occurs and adjusts the entire upper ocean to surface properties.The depth of this convective layer increases with latitude.
At constant depth, buoyancy decreases with latitude everywhere at the western boundary.This is not true for the eastern boundary, where near-surface downwelling advects relatively buoyant surface waters, resulting in nearly neutral stratification reaching deeper than at the western boundary (Fig. 5b).As a consequence, isopycnals slope downward with latitude before they steeply outcrop, forming a bowltype structure.Note that, although stratification is nearly neutral, convective mixing is not present in this outcropping region.The distinct structures at eastern and western boundaries result in an east-west buoyancy difference in thermal-wind balance with the meridional velocity field (Fig. 5c).Note that where the western boundary current transports large amounts of buoyancy poleward, convective mixing is essential for maintaining the eastern boundary's excess in buoyancy, in turn essential for the western boundary current.
The buoyancy patterns found are in good agreement with the ones presented in Marotzke (1997), especially in capturing the distinct structure at the different boundaries and the associated east-west buoyancy difference.The bowltype structure prevailing at the eastern boundary, and hence throughout the interior, was also found by other studies with primitive-equation and planetary geostrophic models (e.g.Bryan, 1987;Colin de Verdière, 1988;Zhang et al., 1992;Winton, 1996).Samelson (1998) argued that eastern diapycnal mixing induces zonal buoyancy gradients near the eastern boundary and thus supports an eastern boundary current.To analyze Samelson's (1998) eastern-boundary mixing case, we reduce the western boundary diffusivity by an order of magnitude and leave everything else as in the standard case.The resulting buoyancy fields are similar to the standard solution ones, except for an approximately 30 % shallower pycnocline (not shown).The western boundary upwelling is decimated, and the overturning strength is reduced by approximately 30 %, but the patterns of all other velocity fields are very similar to the standard solution ones (not shown).That essentially the same overturning structure ensues indicates that even in the eastern-boundary mixing case, an eastern boundary current is not of first-order importance in overturning dynamics.

Parameter dependence
To gain insight into the model's parameter dependence, we now vary the parameter κv over several orders of magnitude and analyze the (nondimensional) solution.With increased κv , the pycnocline broadens and the overturning strengthens.The maximum of the overturning stream function shifts downward and equatorward, which reflects a deepening of the poleward surface western boundary current and a broadening of the downwelling region (Fig. 6).
To quantify this dependence, we consider four characteristics of the solutions: the maximum value of the overturning stream function ψmax , the maximum value of the advective meridional buoyancy transport Ĥmax , and the tropical pycnocline depth defined in two different ways: and, following Park and Bryan (2000), Both definitions correspond to the e-folding depth if buoyancy decays exponentially.Note, however, that the first definition is a measure of surface decay, whereas the second formally requires exponential decay over an infinitely deep ocean for this correspondence to hold.All characteristics mentioned above increase with increasing κv (Fig. 7).For κv ≤ 10 −3 , they follow power laws very well.Linear fits, excluding the solutions for κv > 10 −3 , yield The second definition of the tropical pycnocline depth yields a slightly higher exponent in the power law than the first definition.We argue that the overturing dynamics are determined in the pycnocline rather than in the abyss.In the following, we hence resort to the first definition because it better captures the surface buoyancy structure and does not depend on the abyssal stratification -and hence is a better measure for the pycnocline depth.
Replacing the exponents in the power laws (Eq.35) by simple common fractions, i.e. ψ ∼ κv 2 3 , Ĥ ∼ κv 2 3 , and δ ∼ κv after restoring dimensions.Analogous to Welander (1971), the above scaling relations are now derived, assuming a unique internal depth scale δ.But having noted that zonally sloping isopycnals are confined to a western boundary layer, its width a λ is taken as zonal length scale, rather than the basin width.Two further assumptions must be made: -the east-west buoyancy difference is assumed to scale with the imposed meridional surface buoyancy difference (Marotzke, 1997;Park and Bryan, 2000), and -vertical advection and diffusion are assumed to be of the same order of magnitude in the buoyancy equation (e.g.Welander, 1971;Marotzke, 1997).We hence conclude from the thermal-wind balance for the meridional flow (Eq.9), the continuity equation for the zonally integrated flow (Eq.10), and the buoyancy Eqs. ( 6) and ( 11) that Combining these yields and thus, by the definition of the overturning stream function (Eq.19), The meridional buoyancy transport then scales like All assumptions made to derive this simple scaling are a priori not clear from the model equations.That the model results reported in Eq. ( 36) agree with this simple scaling, however, indicates that the assumptions leading to the scaling are justified for κv ≤ 10 −3 .To explicitly verify that the east-west buoyancy difference scales with the imposed meridional surface buoyancy difference (cf. Park and Bryan, 2000), the maximum eastwest buoyancy difference is diagnosed from the model.As a function of κv , it shows a near-constant relationship over the parameter range the scaling holds in (Fig. 8).The increase at high values of κv is due to the emergence of a strong east-west buoyancy difference in the downwelling region not present in the lowκv regime.
For large values of κv , when the bottom of the basin is felt, the above assumptions must be relaxed.In particular, the assumption of a unique depth scale is not justified anymore.Indeed, the parameter dependence starts deviating from the scaling law for κv > 10 −3 (Fig. 7).Note that, by the definition of κv (Eq.28), we require a sufficiently small vertical diffusivity k v and a sufficiently large imposed meridional surface buoyancy difference b for a sufficiently small κv .
For primitive-equation and planetary geostrophic models (e.g.Bryan, 1987;Colin de Verdière, 1988;Marotzke, 1997), some disagreement concerning the correct power-law dependence of both overturning strength and meridional heat transport on vertical diffusivity prevailed until Park and Bryan (2000) consolidated a large portion of these discrepancies.They pointed out that a restoring boundary condition, as employed in all these studies, allows for some degree of freedom for surface density.Adjusting the restoring time scale as a function of vertical diffusivity, they found good agreement with the Welander (1971) scaling.Our model's dependence on vertical diffusivity agrees with the results presented in Park and Bryan (2000) for all three quantities considered.
The Marotzke (1997) theory, developed for the boundarymixing case, also exhibits a scaling in agreement with both our model results and the simple arguments presented above.Our model results thus also substantiate the arguments put forth to derive the Marotzke (1997) scaling.
Both the Marotzke (1997) theory and the two-plane model results predict a power-law dependence of all considered quantities on the imposed surface meridional density difference.This has not yet been found to correspond to primitiveequation model solutions.Scott (2000) varied the imposed surface temperature difference and adjusted the imposed surface salinity difference proportionally.He found a power-law ψ ∼ T 0.61 , where T is the imposed meridional temperature difference.Lucas et al. (2006) used purely thermal dynamics and varied the shape of the imposed meridional temperature profile.They only found a power-law dependence when they decreased the equatorial temperature, not when they increased the polar temperature.Due to the equation of state's nonlinearity, the temperature profiles they considered do not correspond to self-similar density profiles.In the case of decreased equatorial temperatures, the deviation from self-similarity is comparatively small; in the case of increased polar temperatures, the deviation from self-similarity is comparatively large.This might explain why they found a power-law dependence in the former and not in the latter case.

The role of convective mixing
We now examine the role of convective mixing in the overturning dynamics by allowing static instability to persist, that is by switching off convective adjustment.This is certainly an unrealistic scenario, but it is hoped to reveal -in a dialectic sense -the importance of convective mixing in the overturning dynamics. 3All parameters are chosen as in the standard run.
The most apparent feature of the ensuing steady-state circulation is a two-cell meridional-overturning structure (Fig. 9).At the western boundary, near-surface meridional currents converge around 38 • , and an eastward flow inducing deep sinking at the eastern boundary forms (Fig. 10).At depth, waters flow west, spread out both pole-and equatorward at the western boundary, and upwell at both boundaries to close the respective loop (Fig. 10).
For the standard setup, we found that the abyss is filled with water of the same buoyancy as that of surface waters at the site of deep sinking (Fig. 5).This still holds when convective mixing is switched off -with the consequence that abyssal waters are much more buoyant and that the highlatitude pycnocline is rendered statically unstable (Fig. 11).
That this instability is a necessary property of the nonconvective solution can readily be explained by considering the surface buoyancy budget.In a steady state, the integrated surface buoyancy flux must vanish.In the standard run, this is accomplished by a convective loss of buoyancy at high latitudes, compensating a low-latitude diffusive gain.In the nonconvective case, however, there must be regions in which the ocean diffusively transports buoyancy across the surface, implying statically unstable stratification.
We now argue that, without convective mixing, a one-cell structure overturning can no longer exist.If it did prevail, high-latitude downwelling would keep filling the abyss with the least buoyant surface waters available.At the western boundary, the surface boundary current, transporting buoyancy poleward, would render the buoyancy stratification statically unstable.At the eastern boundary, where only vertical advection is possible, downwelling can at most induce neutral stratification.In effect, buoyancy would be greater 3 Cf.Lindzen's (1990, p. 101) discussion of eddy-free circulation models of the atmosphere as a means to get insight into the role of eddies.at the western than at the eastern boundary, implying an equatorward shear.Western upwelling and eastern downwelling would therefore have to be enhanced in order to ensure a poleward surface current.The strength of the associated zonal circulation, however, has an upper limit.Zonal shear must comply with thermal-wind balance (Eq.8), and the meridional gradient of eastern-boundary buoyancy is limited by the surface boundary condition.
The state the nonconvective system assumes, namely two overturning cells rotating in opposite directions, advects comparably little buoyancy meridionally.Ensuring appropriate shear of the western boundary current, the strong zonal overturning is formed in the middle of the basin, where the meridional gradient of surface buoyancy is greatest.
The overturning's two-cell structure agrees with the nonconvective case of Marotzke and Scott (1999).Their overturning cells are stronger than in our model, however, and the site of deep sinking is located further poleward.
The results of the nonconvective case, which are in stark contrast to the standard solution, underscore the finding that convective mixing is essential for THC dynamics (Zhang et al., 1992).It ensures that deep sinking occurs at high latitudes and that the abyss fills with the densest surface waters available.In contrast, Marotzke et al. (1988) found convective mixing not to play a crucial role in their zonally averaged model.Zhang et al. (1992) explained this by pointing out that the closure employed for the zonal pressure gradient does not allow for its reversal; but this reversal is crucial in nonconvective dynamics.Marotzke and Scott (1999) argued that all closures employing a local and linear relation between the meridional pressure gradient and flow strength necessarily yield a collocation of convective mixing and deep sinking sites.In line with three-dimensional models, our standard solution involves convective mixing at the western but deep sinking at the eastern boundary -overcoming a fundamental deficiency of zonally averaged models.

Main model results
The two-plane model formulated in Sect. 2 exhibits a steadystate THC that shows all salient properties of velocity and density fields diagnosed from a primitive-equation model considering the boundary-mixing case (Marotzke, 1997).In its standard setup the two-plane model yields narrow deep sinking in the eastern high latitudes and broad upwelling throughout both eastern and western boundary layers, as well as a poleward surface and an equatorward deep western boundary current.A single overturning cell of strength 15.1 × 10 6 m 3 s −1 prevails.The buoyancy fields show a pycnocline structure with slightly higher values in the eastern-boundary pycnocline, resulting in a bowl-type structure there.
Varying the nondimensional parameter κv yields a powerlaw dependence of approximately ψ ∼ k v 2 3 b 1 3 for sufficiently low diffusivities k v and sufficiently large imposed meridional surface buoyancy differences b.For larger diffusivity, when the pycocline is not small compared to the basin depth anymore, the model solutions show a notable deviation from the power-law dependence.
Allowing static instability to persist alters the overturning substantially.A two-cell structure arises, and downwelling of relatively buoyant waters occurs at midlatitudes.These relatively buoyant waters fill the abyss, and an inverse pycnocline forms at high latitudes.

Conceptual picture of the overturning dynamics
We now give a conceptual picture of the THC dynamics consistent with boundary mixing.Note that causality is not an appropriate notion here -we rather try to explain the inevitable patterns of a self-consistent steady-state circulation.
To start our discussion, we assume a pycnocline structure: the abyss is filled with the least buoyant surface water available, and higher buoyancy values are confined to a thin surface layer.This is already where convective mixing comes into play: it rules out any other abyssal buoyancy by prohibiting static instability.
Diapycnal mixing diffuses buoyancy down the pycnocline, necessitating advective transports to accomplish steady state.As zonal buoyancy gradients cannot prevail at the eastern boundary for dynamical reasons (Rossby wave activity), meridional transports do not enter the balance here.In effect, upwelling of less buoyant abyssal waters alone counterbalances the downward diffusive transport (Fig. 3d).
In the absence of meridional flow, vertical convergence of upwelling waters at pycnocline-depths must be compensated for by westward flow; vertical divergence in the abyss must be compensated for by eastward flow (Fig. 3a).At the western boundary, near-surface convergence of upwelling waters is therefore enhanced by zonal inflow at pycnocline-depths, and a poleward boundary current must result (Fig. 3b).
This meridional current requires a positive east-west buoyancy difference by thermal-wind balance.Convective mixing supplies an effective buoyancy sink at the western boundary.More buoyant waters at the eastern boundary, however, can only be attained by downwelling of buoyant surface waters, resulting in a two-layer structure in the eastern-boundary vertical velocity field (Fig. 3d) and steep isopycnals (Fig. 5b).
The two-layer structure in eastern-boundary vertical velocity corresponds to a three-layer structure in the zonal velocity field.In the west, near-surface zonal outflow in turn relaxes the need for a boundary current.
The meridional surface current intensifies poleward over large parts of the basin.Only near the polar boundary, convergence must occur.Here, a strong eastward flow emerges, feeding the strong downwelling in the eastern high-latitude ocean.Mass balance at the eastern boundary requires a deep westward return flow feeding a deep equatorward flow in the west, closing mass balance there.
To understand why deep downwelling occurs at the eastern boundary, again thermal-wind balance can be invoked: downwelling at the western boundary would advect relatively buoyant surface waters downward, resulting in a negative east-west buoyancy gradient -inconsistent with a poleward shear of the western boundary current.

Outlook
The applicability of the two-plane model to the real ocean is constrained by two major omissions.First, the basin is confined to a single hemisphere; and second, temperature and salinity are not considered separately.
Extending the model equations to an inter-hemispheric basin should in principle be possible (Marotzke and Klinger, 2000).The breakdown of the geostrophic approximation at the Equator raises the question of how to incorporate equatorial dynamics.Observational studies (Lukas and Firing, 1984), however, have indicated that, for long time scales, the meridionally differentiated geostrophic balance applies at the equator, which might be exploited in the two-plane model.The equatorial thermal-wind relations where all symbols are the same as in Sect.2.2, might hence successfully be used (Lukas and Firing, 1984).Another possibility is the inclusion of viscous terms.Eventually, a representation of the Southern Ocean is highly desirable (Gnanadesikan, 1999).Considering two active tracers, temperature and salinity, and applying mixed boundary conditions is expected to modify the dynamics of the two-plane model substantially.We hope to then gain insight into questions of multiple equilibria (cf.Stommel, 1961;Bryan, 1986;Marotzke et al., 1988;Marotzke and Willebrand, 1991;Johnson et al., 2007).A thermodynamic coupling to a zonally averaged atmosphere model also appears feasible.

Conclusions
The dynamics of the THC can be cast into a two-plane model based on boundary mixing and geostrophy.In a hemispheric sector basin geometry, this two-plane model predicts self-consistent steady-state density and overturning structures that resemble primitive-equation model solutions.
The two-plane model confirms the scaling proposed by the Marotzke (1997) theory.The overturning strength has a 2/3-power-law dependence on vertical diffusivity and a 1/3power-law dependence on the imposed meridional surface density difference.The dependence on vertical diffusivity agrees with primitive-equation model results.Adapting the simple arguments by Welander (1971) to the boundarymixing case, this scaling can successfully be explained.
In contrast to previous two-dimensional models but in line with three-dimensional models, convective mixing is an essential element of the two-plane model.This indicates that the two-plane model better represents the THC dynamics than previous two-dimensional models.
errors as finite differences are taken across the boundaries of convective regions, resulting in oscillatory vertical velocity fields, especially at the western boundary.Increasing resolution reduces these oscillations, and we are confident that convergence is warranted.

A2 Initial guess
Assuming the steady-state solution to be unique, we can prescribe an arbitrary initial guess fulfilling the boundary conditions for the time-stepping scheme.For all results presented in this paper, the state with δ = 10 −3 is used.Starting from initial guesses with different choices of δ indeed yields the same steady-state solution, supporting the uniqueness assumption.
Only numerical stability puts certain constraints on the choice of the initial guess.Certainly, a sensible initial guess does not differ too much from the expected solution.Nevertheless, even highly irrational initial guesses such as one with virtually vertical isopycnals prove to converge to the same solution.

A3 Sensitivity to artificial diffusion and resolution
To assess the solution's sensitivity to numerical diffusion, we vary κh between 2 × 10 −4 and 2 × 10 −1.5 for different values of κv and analyze the solutions in terms of the four characteristics introduced in Sect.3. We find that the solutions do not substantially depend on κh for κh < 10 −3 .For our choice of κh = 2×10 −4 , the solution has hence sufficiently converged.
In order to provide some indication that the discretization scheme converges to a unique solution for ever finer grids, we analyze approximate solutions to the standard setup of varied resolution.The patterns of buoyancy and velocity do not change substantially when the resolution is varied between 32 × 32 and 256 × 256 grid points, but seem to converge to a unique solution.To quantify this convergence, again the four characteristics introduced in Sect. 3 are analyzed.Indeed, they all indicate convergence to a definite value.

Fig. 1 .
Fig. 1.Envisioned basin geometry and zonal structures of buoyancy and zonal velocity as assumed in the model.Note that the spherical geometry of the model is not represented in this sketch.

Fig. 7 .
Fig. 7. Dependence on κv of the four characteristics (a) maximum value of the overturning stream function ψmax , (b) maximum value of the advective meridional buoyancy transport Ĥmax , (c) tropical pycnocline depth δ1 , and (d) tropical pycnocline depth δ2 (all nondimensional).The red lines are least square fits to the values for κv ≤ 10 −3 , the blue lines' slopes correspond to the theoretical scaling.The numbers denote the slope of the respective lines.

Table 2 .
External parameters as employed in the standard setup.