## Abstract

A recently proposed reduced-gravity model of the warm-water branch of the middepth meridional overturning circulation in a rectangular basin with a circumpolar connection is extended to include time dependence. The model describes the balance between gain of warm water through northward Ekman advection across the circumpolar current, loss of warm water through eddy fluxes southward across the current, net gain or loss of warm water through diabatic processes north of the current, and changes in the thickness of the warm-water layer. The steady solutions are the same as those found previously, when the previous parameterization of diabatic fluxes is used. Time-dependent solutions are considered for the approach of the solution to a new equilibrium when the forcing or parameters are abruptly changed and then held fixed. An initial adjustment occurs through a combination of boundary and equatorial adjustment, followed by planetary wave propagation. The longer-term adjustment to equilibrium consists primarily of the slow change in eastern boundary thickness of the warm layer, which controls the mean depth of the entire layer. An approximate analytical solution of the time-dependent equations yields an explicit expression for the intrinsic time scale of the long-term adjustment, which depends on the eddy and diabatic flux parameters and on the equilibrium solution toward which the time-dependent solution adjusts. Numerical solutions are also considered with a second, advective–diffusive diabatic flux parameterization. These solutions differ quantitatively but not qualitatively from those with the original parameterization. For the range of parameter values considered, the adjustment time scale has dimensional values of several decades to more than a century, but the meridional flux of warm water may respond to changes in external parameters or forcing much more rapidly than this time scale for equilibration of the eastern boundary thickness and thermocline structure.

## 1. Introduction

The ocean’s large-scale meridional overturning circulation is an important element of the earth’s climate system. Our understanding of basic aspects of this circulation, such as the processes and quantities that determine its intrinsic adjustment time scales, remains extremely limited. Here, these intrinsic dynamical time scales are studied in an idealized model of the middepth (Reid 1981) meridional overturning circulation. This middepth cell, which lies between the wind-driven midlatitude cells of the upper thermocline and the abyssal flow at the deepest levels, includes the northward flow of warm thermocline waters from the Southern Hemisphere to the northern subpolar gyre and the compensating southward flow of North Atlantic Deep Water from northern subpolar latitudes to the Southern Hemisphere.

The model discussed here focuses on the warm, northward-flowing thermocline waters and the associated wind-driven, upper-thermocline gyre structure; the compensating southward middepth flow is not represented explicitly. The formulation extends the simple, reduced-gravity model of the warm-water branch recently proposed by Samelson (2009) to allow time dependence of the thickness of the warm-water layer. The temporal adjustment of the resulting model is studied under conditions in which external forcing parameters are abruptly changed and then held fixed. The steady dynamics of the Samelson (2009) model reflect balances similar to the scaling arguments proposed by Gnanadesikan (1999) and Johnson et al. (2007), in which the warm-water volume and flux are established through a balance between gain of warm water through Ekman advection across a circumpolar current, loss of water through eddy fluxes southward across the circumpolar current, and net gain or loss of warm water through diabatic processes north of the circumpolar current (e.g., Deacon 1937; Stommel and Webster 1962; Veronis 1978; Tziperman 1986; Warren 1990; de Szoeke 1995; Toggweiler and Samuels 1995). The extended, time-dependent model presented here allows temporal changes in the thickness of the warm-water layer so that the volume and flux of warm water can evolve dynamically.

## 2. Model

The model formulation and notation follows Samelson (2009), except that, following the standard planetary geostrophic approximation for a reduced-gravity vertical structure, the local rate of change of layer thickness is retained in the mass balance. Thus, the warm-water branch of the meridional overturning circulation is represented by a single homogeneous layer of depth *h*, subject to the direct action of wind stress and resistance from a linear friction, on a *β* plane in a rectangular domain {*x _{W}* <

*x*<

*x*,

_{E}*y*<

_{S}*y*<

*y*} (Fig. 1). The corresponding zonal and meridional momentum equations are given by

_{N}and the mass conservation equation is

where *τ* = (*τ ^{x}*,

*τ*) is the kinematic wind stress,

^{y}*r*is a constant friction coefficient, and the form of

*W*remains to be specified. A second, more minor difference from the formulation of Samelson (2009) is that a second, alternative form of

*W*is considered here, in addition to the original form; both are described below. The geometry is basic: a rectangular basin with vertical sidewalls that extends across the equator, at

*y*= 0, into both hemispheres. Periodic zonal boundary conditions allow a zonally reentrant representation of a circumpolar current in a narrow band of high Southern Hemisphere latitudes,

*y*

_{1}<

*y*<

*y*

_{2}< 0.

where

and **∇** is the horizontal gradient operator. Conditions on the normal velocity component **u** · **n** at the boundary thus lead to a differential equation around the boundary, which relates the normal and tangential derivatives of *h* (Salmon 1986). Along the contiguous rigid boundary, the standard no-normal-flow condition **u** · **n** = 0 is applied. This is supplemented by representations of Ekman and eddy fluxes of warm water across the open boundary of the warm layer at *y* = *y*_{2}, the southern edge of the warm layer and the northern edge of the circumpolar current. As in Samelson (2009), an analytical, zonally symmetric, geostrophic circumpolar current can be introduced in the region *y*_{1} < *y* < *y*_{2}, at the southern edge of the warm-water layer, to complete the model’s large-scale thermal structure.

The Ekman transport per unit longitude across *y* = *y*_{2} is given by

where *f*_{2} = *f*(*y*_{2}). A simple prescription of the eddy flux *V _{e}* of warm water across

*y*=

*y*

_{2}is adopted,

where *α _{e}* is a constant of proportionality and is the zonal mean of the squared depth of the warm-water layer on the north side of the current, at

*y*=

*y*

_{2};

*V*is thus independent of longitude

_{e}*x*. The primary motivation for the prescription (7) is analytical convenience, but it may also be interpreted as the zonally averaged depth integral of fluxes proportional to isopycnal slopes

*z*/

_{ρ}*L*

_{acc}, where

*z*ranges from 0 to

_{ρ}*h*

_{m}in a local, stratified representation of the warm layer and

*L*

_{acc}is a fixed width for the circumpolar current; the latter interpretation shows the connection to the linear gradient-flux schemes that are widely used to parameterize eddy fluxes in large-scale numerical models. In the solutions considered here,

*h*(

*x*,

*y*

_{2}) is itself nearly independent of

*x*, so a local prescription

*V*= −

_{e}*α*(

_{e}h*x*,

*y*

_{2})

^{2}would yield a nearly identical flux.

With the conditions (6) and (7) at *y* = *y*_{2} and the no-normal-flow condition along the rigid boundary, (4) may be solved to obtain **u**(*x*, *y*) from a given *h*(*x*, *y*). The remaining equation (3) then determines *h _{t}*, once the diabatic velocity

*W*(

*x*,

*y*) is specified in terms of

*h*. Here,

*W*is first taken to be proportional to the difference of the square of the layer thickness

*h*and an imposed function , as in Samelson (2009):

where *α _{w}* is a given constant of proportionality and is a given function of

*x*and

*y*. A second form is also briefly explored,

where *κ* is a given constant with the dimensions of a diffusivity. The form (9) is motivated by the scaling for the vertical advective–diffusive balance in the thermodynamic equation,

according to which *W _{κ}* =

*κ*/

*δ*(Welander 1971) and is the parameterization employed by Gnanadesikan (1999), Johnson et al. (2007), Allison (2009), and Allison et al. (2011, manuscript submitted to

*J. Mar. Res*). The relation between (8) and (9) is discussed below in section 5.

The equations may be made dimensionless using the same dimensional scales as in Samelson (2009): horizontal length scale *L* = 5000 km, depth scale *H* = 1000 m, gravitational acceleration *g* = 9.8 m s^{−2}, density *ρ*_{0} = 1027 kg m^{−3}, Coriolis parameter *f*_{0} = 10^{−4} s^{−1} and its meridional gradient *β*_{0} = 2 × 10^{−11} m^{−1} s^{−1}, wind stress scale *τ*^{*} = 0.1 N m^{−2}, horizontal velocity scale *U* = *τ*^{*}/(*ρ*_{0}*f*_{0}*H*^{*}) = 10^{−3} m s^{−1}, and time scale *T*_{adv} = *L*/*U*. The scale for volume transport is then *UHL* = 5 × 10^{6} m^{3} s^{−1}. The resulting dimensionless equations have the same form as the dimensional equations, with the parameters *γ*, , *α _{e}*,

*α*, and

_{w}*r*scaled by

*f*

_{0}

*UL*/

*H*,

*H*

^{2},

*U*/

*H*,

*U*/(

*HL*), and

*f*

_{0}, respectively. A density difference Δ

*ρ*/

*ρ*

_{0}= 10

^{−3}is assumed, giving a dimensionless value

*γ*= 20.

The spatial structure of the forcing functions is the same as in Samelson (2009): the imposed wind stress *τ* is taken to be purely zonal,

and *τ ^{y}* = 0, and the squared equivalent thickness function for the diabatic forcing north of the circumpolar channel is given by

The dimensionless domain has *x _{W}* = 0,

*x*= 1,

_{E}*y*

_{2}= −0.6,

*y*= 1, and

_{N}*f*= 0 at

*y*= 0. For numerical solution, the equations were discretized on a 51 × 81 grid and solved by a simple second-order Runge–Kutta scheme, with a tridiagonal matrix solver for the differential equation around the boundary.

Steady-state solutions of these equations may be obtained by setting *h _{t}* = 0 and solving for

*h*(

*x*,

*y*). With the choice (8) for

*W*, the steady-state equations have the same form as those solved by Samelson (2009). The analysis here focuses on the time-dependent adjustment to these steady states, in response to abruptly changed forcing or parameter values, which are then held fixed until the new steady solution is reached.

## 3. Numerical solutions

### a. Initial state

A sequence of two time-dependent numerical solutions with different, fixed parameters illustrates the basic characteristics of the evolution of the system toward steady state. The sequence is initialized with the steady solution for a third set of parameters, which, except for the friction parameter *r*, are the same as the first steady solution described by Samelson (2009). This initial state has weak midlatitude upwelling and equatorially symmetric zonal wind forcing. For the first time-dependent solution, the Northern Hemisphere cooling parameter is abruptly increased and then held fixed. This solution models the large-scale, time-dependent response of the warm-water branch of middepth meridional overturning to increased deep-water formation at the high latitudes of the Northern Hemisphere. For the second time-dependent solution, after the system reaches the modified steady state with increased Northern Hemisphere cooling, the Southern Hemisphere winds at the latitude of the model circumpolar current are abruptly increased and then held fixed. This solution models the large-scale, time-dependent response of the warm-water branch of middepth meridional overturning to increased Southern Hemisphere wind forcing. For these solutions, the form (8) is used for the diabatic velocity *W*(*x*, *y*).

The initial state is obtained for dimensionless parameter values *τ*_{0} = −1, *τ*_{1} = 0, , *α _{w}* = 2,

*α*= 1, ,

_{e}*γ*= 20, and

*β*= 1. To enhance stability of the numerical integration scheme, the frictional parameter is increased in the vicinity of the equator so that the constant

*r*in (1) and (2) is replaced by . For the solutions discussed here, the corresponding dimensionless parameter values are

*r*

_{0}= 0.05,

*r*

_{eq}= 0.25, and

*L*= 0.1. The initial state thus differs from the first solution discussed by Samelson (2009) only by this stronger, spatially variable friction , which replaces the previous constant value

_{r}*r*= 0.02.

The basic characteristics of the first steady solution discussed by Samelson (2009; Figs. 2, 3) are maintained in this modified numerical solution, with some minor differences because of the stronger friction. There are wind-driven subtropical gyres in both hemispheres and a wind-driven subpolar gyre in the Northern Hemisphere, with the accompanying deformations of the warm-layer thickness *h* (Fig. 2a). The zonally integrated meridional flow *V*(*y*), where

is northward at all latitudes in the warm layer (Fig. 2c). The northward Ekman transport across the circumpolar current is partially compensated by opposing eddy transport of warm water so that the overturning flow at the southern edge of the warm layer is roughly one-half of the northward Ekman transport. This overturning flow travels northward to the Northern Hemisphere subpolar gyre, where it cools and sinks. There is weak subtropical and tropical upwelling (Fig. 2b), so the northward warm-water transport increases northward and the high-latitude Northern Hemisphere downwelling is slightly more than the partially compensated Ekman transport into the southern edge of the warm layer (Fig. 2c). The layer thickness increases westward from the eastern boundary in the interior of the subtropical gyres, resulting in slightly decreased diabatic upwelling toward the west, and decreases westward in the subpolar gyre. Because of the larger friction, the midlatitude thickness *h* is slightly smaller and the weak midlatitude upwelling is consequently slightly stronger than in the Samelson (2009) *r* =0.2 solution (Fig. 2b). Similarly, the large value of *r*_{eq} results in relatively viscous equatorial dynamics, but this does not cause dramatic changes in the solution structure.

### b. Adjustment to strong NH cooling

For the first time-dependent numerical solution, the Northern Hemisphere cooling is abruptly intensified by changing the diabatic forcing parameter from to at time *t* = 0 and then holding it fixed at the new value. The initial state for this calculation is the steady solution described above (Fig. 2). Except for this change to and the increased, spatially variable friction, the parameters for the time-dependent solution are again the same as for the first steady solution described by Samelson (2009; Figs. 2b, 3a,c).

The final state from this time-dependent adjustment again has the standard subtropical and subpolar-gyre thickness structure, downwelling in the northern subpolar gyre, and northward meridional transport at all latitudes (Fig. 3). Relative to the initial state (Fig. 2), it has enhanced meridional warm-water transport and Northern Hemisphere sinking (Fig. 4). There is stronger subtropical and tropical upwelling than in the initial state, so the northward warm-water transport increases northward by roughly one-third from the circumpolar current latitudes to the Northern Hemisphere subpolar-gyre boundary (Fig. 4c), with strong cooling and sinking in the northern part of the subpolar gyre. There is a mean decrease in the warm-layer thickness (Fig. 4a), relative to the initial state, which supports the subtropical and tropical upwelling (Fig. 4b). This decrease also reduces the southward eddy transport across the circumpolar current so that the net inflow at *y* = *y*_{2} is increased relative to the initial state (Fig. 4c). Thus, the greater overturning and Northern Hemisphere sinking induced by the intensified Northern Hemisphere cooling is supported by a combination of increased net inflow across the circumpolar current and weak, broad upwelling.

The time-dependent adjustment from the initial to the final state proceeds in several stages, which in general outline are familiar from previous analyses of large-scale adjustment, such as that of Kawase (1987). There are essentially two distinctive time scales, the planetary wave basin-crossing time *T*_{pwave} ≈ −*L*/*c _{R}* (where

*c*= −

_{R}*βγh*/

*f*

^{2}; see section 4), and the long-term advective adjustment time scale

*T*

_{adv}≈

*L*/

*U*. For characteristic scales as above,

*T*

_{pwave}corresponds to dimensional times of several years to a decade and

*T*

_{adv}is on the order of a century or longer.

On short time scales *t* ≪ *T*_{pwave}, there is a zonally symmetric increase in Northern Hemisphere cooling and sinking (*W* < 0) and a corresponding zonally symmetric decrease in layer thickness, which extends southward along the western boundary to the equator (Fig. 5). The enhanced divergence of the northward transport, which supplies the enhanced northern sinking, is supported by this southward extension of the region of decreasing layer thickness. An eastward penetration along the equator of the leading edge of the disturbance is visible in the height field during this period (Fig. 5a). Note that, for a given subpolar-gyre western boundary disturbance, the propagation southward along the western boundary is effected by the solution of the differential equation around the boundary that results from the no-normal-flow condition, which supports the steady, damped model representation of Kelvin-type coastal boundary wave dynamics (see appendix).

On the longer but still short (*t* ≈ *T*_{pwave} ≪ *T*_{adv}) time scales of planetary wave basin-crossing times, the familiar zonal slope of the layer thickness associated with the Sverdrup balance is established (Fig. 6). The pathway of the signal eastward along the equator to the eastern boundary, then poleward along the eastern boundary, and westward via planetary waves is apparent in the height disturbance field during this period (Fig. 6a).

Finally, on the full (*t* ≈ *T*_{adv}) adjustment time scale, the mean layer thickness evolves toward the equilibrium state, whereas the horizontal structure of the layer thickness variations remains roughly constant. In view of the analytical solution discussed below, it is natural to measure this change in mean layer thickness in terms of the eastern boundary layer depth *h*_{E} = *h*(*x _{E}*,

*y*). This adjustment is essentially an exponential decay, with a time scale of order 0.2, corresponding to dimensional times of several decades (Fig. 7a). As demonstrated below by an approximate analytic solution, it is the eddy and diabatic flux parameters, along with the steady solution itself, that determine this adjustment time scale.

### c. Adjustment to strong SH winds

For the second time-dependent numerical solution, the high-latitude Southern Hemisphere winds are abruptly intensified, by changing the wind forcing parameter *τ*_{1} from *τ*_{1} = 0 to *τ*_{1} = 0.5 at time *t* = 0 and then holding it fixed at the new value. The initial state for this calculation is the steady solution that is the final state of the first time-dependent numerical solution (Fig. 3). Except for the increased value of *τ*_{1}, the parameters for this second time-dependent solution are the same as for the first time-dependent solution.

The final state from this time-dependent adjustment (Fig. 8) is similar to the original initial state (Fig. 2), but with a larger meridional overturning circulation. As in the original initial state, there is weak midlatitude and equatorial upwelling (Figs. 8b,c), with the meridional overturning circulation strengthening progressively northward toward the Northern Hemisphere subpolar gyre so that the high-latitude Northern Hemisphere downwelling is slightly greater than the partially compensated Ekman transport into the southern edge of the warm layer. The latter is nearly twice that for the original initial state because of the increased Southern Hemisphere winds. Relative to the initial state for the time-dependent adjustment, the additional meridional transport divergence for the final state is nearly constant with latitude (Fig. 9c). Thus, the additional Ekman transport into the warm layer, from the intensified Southern Hemisphere winds, is balanced by a nearly uniform decrease in the weak midlatitude upwelling, supported by a nearly uniform increase of 0.1 units (100 m) in the layer thickness (Figs. 9a,b). There is a small region in the western part of the Southern Hemisphere subtropical gyre, in which this decrease is sufficient to cause local downwelling; this is probably not physically realistic, but it has only a small effect on the solution.

As in the previous case of cooling, the time-dependent adjustment from the initial to the final state for the intensified Southern Hemisphere winds proceeds in several stages (Kawase 1987). On short time scales, *t* ≪ *T*_{pwave}, there is a zonally symmetric increase in the thickness of the Southern Hemisphere subtropical gyre because of Ekman convergence from the intensified winds. This signal extends northward along the western boundary to the equator, although the extension is less dramatic than in the cooling case, because the forced response itself directly affects the subtropical gyre (Fig. 10). An enhanced convergence of the northward transport, induced by the enhanced northward Ekman transport at *y* = *y*_{2}, extends northward to the equator and into the Northern Hemisphere. An eastward penetration along the equator of the leading edge of the disturbance is again visible in the height field during this period (Fig. 10a).

On the longer but still short (*t* ≈ *T*_{pwave} ≪ *T*_{adv}) time scales of planetary wave basin-crossing times, with *T*_{pwave} again corresponding to dimensional times of several years to a decade, the zonal slope of the layer thickness is again established (Fig. 11). The pathway of the signal is again eastward along the equator to the eastern boundary, then poleward along the eastern boundary, and westward via planetary waves.

Finally, on the full (*t* ≈ *T*_{adv}) adjustment time scale, the mean layer thickness again evolves toward the equilibrium state, whereas the horizontal structure of the layer thickness variations remains roughly constant. This adjustment is again exponential, with a time scale on the order of 0.2, corresponding to dimensional times of several decades (Fig. 7b).

## 4. Asymptotic analysis

### a. Approximate solution

The main object of the asymptotic analysis is to derive an approximate solution for the approach to steady state of the time-dependent solution with fixed forcing and parameter values. This analysis extends the results of Samelson (2009), which gave analytical approximations for the corresponding steady-state solutions, which were accurate when the friction coefficient *r* and a parameter proportional to *α _{w}* were both small. One additional condition is necessary to obtain the reduced equation in the present case: namely, the cross-basin propagation time scale for the reduced-gravity planetary waves must be short compared to the time scale for approach of the large-scale solution toward steady state. With an additional approximation that the volume of the warm-water layer can be estimated from the product of the eastern boundary thickness

*h*

_{E}and the basin area, an explicit analytical solution can be readily obtained that describes the approach to steady state and its dependence on the forcing and parameter values.

As *r* → 0, with *τ ^{y}* = 0, the no-normal-flow condition and (2) together imply that the eastern boundary depth

*h*

_{E}of the warm layer is constant along the boundary and thus can depend only on time

*t*,

It is the eastern boundary depth that appears naturally in the steady-state solutions, because it controls the mean value of the layer thickness *h* and so is determined by the warm-water volume balance that closes the steady problem (Samelson 2009; see also de Szoeke 1995; Johnson et al. 2007). It will similarly play a central role in the time-dependent solution derived here.

As *r* → 0, again, the momentum equations (1) and (2) may be divided by *f* and substituted into the mass conservation equation (3), yielding the modified, nonlinear long planetary wave equation,

where

Now, suppose that *h*^{2} may be written as

where is a function to be determined, and substitute this expression into (15) to obtain

Because *dh _{E}*/

*dt*should be less than

*H*/

*T*=

*HU*/

*L*= 1000 m/100 yr = 7 m yr

^{−1}, whereas typical midlatitude values of

*W*

_{Ek}are on the order of 30 m yr

^{−1}, it is consistent to neglect

*dh*/

_{E}*dt*in (18). This is equivalent to the assumption that the time scale for adjustment of

*h*

_{E}is long compared to the basin-crossing time at the long planetary wave speed

*c*= −

_{R}*βγh*/

*f*

^{2}. In this case, the function is determined to first order by the steady-state balance,

where *λ* = 2*f*^{2}α_{w}/(*βγ*). In addition, if *α _{w}* is small enough that

*λ*(

*x*−

_{E}*x*) ≪ 1, the terms proportional to

_{W}*α*in (19) may also be neglected (Samelson 2009), yielding

_{w}Thus, the instantaneous structure of *h*(*x*, *y*, *t*) may be rationally approximated in this limit by the familiar steady-state, adiabatic, wind-driven Sverdrup theory, provided that the instantaneous value *h*_{E}(*t*) of the eastern boundary thickness is known.

An equation for the rate of change of *h*_{E} may be obtained by integrating the mass conservation equation (3) over the extent of the warm-water layer. With the boundary conditions (6) and (7) on transport across *y* = *y*_{2}, this yields

where

Thus, provided that ∂(*τ ^{x}*/

*f*)/∂

*y*= 0 at

*y*=

*y*

_{2}so that , as in Samelson (2009), it follows that

where

For *h _{t}* = 0, (23) yields the steady-state solutions of Samelson (2009), whereas, in general,

The latter integral depends on *h*_{E} and is not easily evaluated analytically. However, a reasonable first approximation for the relatively large-scale thicknesses *H* = 1000 m of the warm-water layer is

That is, the effective average of the wind-driven deformations over the warm-layer extent may be taken to be small relative to the eastern boundary thickness *h*_{E}.

The resulting ordinary differential equation for the evolution of *h*_{E} is

where the parameter *μ* is given by

and *h _{Es}* is the steady solution of Samelson (2009),

Equation (28) is nonlinear but may be readily solved for *h*_{E}(*t*),

where the constant Δ*h*_{E} is the difference of the initial value of *h*_{E}(*t*) from *h*_{Es},

Thus, the time-dependent solution approaches the steady solution exponentially, with decay time scale

### b. Discussion

The approximate analytical solution in (29)–(31) has several immediate implications. First, under conditions of steady forcing and parameter values, the approach to the unique steady-state solution is monotonic; despite the nonlinearity, there are no oscillations and no multiple equilibria. Second, the nonlinearity does cause a dependence of the approach time scale *T*_{MOC} [(33)] on the steady-state solution value of the eastern boundary thickness, such that the approach to solutions with relatively larger eastern boundary thicknesses will be relatively faster. Third and perhaps most important, the approach time scale *T _{μ}* is directly proportional to the eddy flux and diabatic parameters

*α*and

_{e}*α*.

_{w}Despite the approximations involved, the analytic solution accurately describes the time-dependent adjustment of both of the numerical solutions described in the preceding section (Fig. 7). A small error is incurred if the analytic solution is used to determine the steady-state values of *h*_{E}. If the functional form (31), with *μ* given by (29), is used with the initial and final values of *h*_{E} taken from the numerical solution, the resulting time evolution of *h*_{E} reproduces closely that from the numerical solution (Fig. 7).

The parameterizations of eddy and diabatic fluxes used here are crude, so the estimates of adjustment time scale for the thermocline structure, as represented by the layer thickness, may not be quantitatively accurate. The primary significance of the results is the qualitative demonstration of the direct dependence of the adjustment time scale on these parameterizations. This result should be robust and should not depend on the particular parameterizations chosen here. Similar dependencies of adjustment time scale on parameterized processes can therefore be anticipated in much more sophisticated numerical models of global ocean circulation, such as those used in coupled climate system models.

A related extension of the scaling arguments of Gnanadesikan (1999) and Johnson et al. (2007) to include time dependence has recently been proposed by Allison (2009) and L. C. Allison et al. (2010, personal communication). Although the form of the scaling solution is different in detail, the qualitative result is similar. Namely, the adjustment time scale for the large-scale, middepth ocean thermocline—framed in terms of the spin-up of the circumpolar current, rather than the warm-water volume as here—is found as here to depend on the eddy and diapycnal flux parameters.

It is important to note that, although (29)–(31) describe the intrinsic time scale of thermocline adjustment associated with the middepth meridional overturning, the adjustment of the northward flux of warm water that composes the warm-water branch of the middepth overturning may occur more rapidly. For example, almost half of the change in northward flux occurs during the first half year in the Northern Hemisphere cooling case (Figs. 4c, 5c) and even more in the Southern Hemisphere winds case (Figs. 9c, 10c). This is due to the strong control exerted on the meridional flux of warm water by changes in external forcing or parameters, and it is also due to the essentially instantaneous development of meridional flow in the western boundary layer. In the Northern Hemisphere cooling case, the effect on *W* of the instantaneous change in *h*_{*} is comparable in magnitude to that due to the adjustment of *h*. In the Southern Hemisphere winds case, the instantaneous change in *T*_{Ek} is generally larger than the change in *W* or *T _{e}* that is due to the adjustment of

*h*. These rapid local changes in warm-water volume formation and flux are communicated zonally by long planetary waves and meridionally by western boundary dynamics on shorter time scales than the basin-wide adjustment of the eastern boundary depth. Thus, the volume and heat fluxes associated with the warm-water branch of meridional overturning may respond much more rapidly than the thermocline structure.

## 5. Advective–diffusive midlatitude diabatic flux

An alternative form for the diabatic flux parameterization is *W _{κ}* =

*κ*/

*h*, as given in (9). This form is motivated by scaling for heat diffusion through the subtropical main thermocline and is less appropriate for the cooling and sinking regions of the Northern Hemisphere subpolar gyre. Thus, a combination of the two forms, which retains the original (8) in the subpolar gyre, is considered here,

where

so that *A*(*y*) = 1 south of the subpolar gyre (*y* < *y*_{SP}) and *A*(*y*) = 0 in the subpolar gyre (*y* > *y*_{SP}), where Δ*y* ≪ 1 is a transition distance that is small compared to the gyre scale.

For the solutions considered above, the deviations Δ*h* of *h* in the low and middle latitudes from a reference value *h _{r}* are generally small, Δ

*h*/

*h*≤ 0.1. Thus,

_{r}*W*may be approximated by

_{κ}Similarly, the original form (8) of *W* in the tropical and subtropical regions may be approximated by

Thus, the two parameterizations will be essentially equivalent if

The dimensional scale for the corresponding diffusivity *κ* is *UH*^{2}/*L* = 2 × 10^{−4} m^{2} s^{−1}, so appropriate values of the dimensionless *κ* for the subtropical main thermocline are 0.05 ≤ *κ* ≤ 0.5. By (38), for *h _{r}* ≈ 1, the equivalent range of values for

*α*would be 0.025 ≤

_{w}*α*≤ 0.25. This is substantially smaller than the value

_{w}*α*= 2 used in the solutions above. However, the larger value is arguably appropriate for the warm-to-cold conversion in the subpolar gyre and thus was retained in the solutions above. Although the scaling leading to the form

_{w}*W*is crude, this comparison does hint that the intensity of midlatitude diabatic exchange may be exaggerated in those solutions.

_{κ}Numerical solutions with the modified diabatic flux from (34) are qualitatively similar to those with the original flux (8) but have some quantitative differences. With *κ* = 0.05, corresponding to a dimensional *κ* = 10^{−5} m^{2} s^{−1}, and all other parameters as in the second solution discussed above (, *τ*_{1} = 0), the time scale for approach to equilibrium is more than doubled (Fig. 12) because of the weaker midlatitude diabatic fluxes. To maintain numerical stability, the friction parameter *r*_{0} is increased from 0.05 to 0.06 for this solution, but this has a negligible impact on the adjustment time scale and only a small effect on the steady solution. The layer thickness for the resulting *κ* = 0.05 steady solution is shallower, with little diabatic upwelling in the tropics and subtropics so that the meridional overturning transport is nearly equal to the net flux into the warm layer across the circumpolar current (Fig. 13). Because the layer is shallower than in the previous solution (Fig. 8), the southward eddy transport is reduced, and this net flux is larger than in the previous solution. Consequently, the total overturning flux of warm water into the Northern Hemisphere subpolar gyre is close to that in the solution with *τ*_{1} = 0.5 (Figs. 13c, 8c). The pattern of sinking in the Northern Hemisphere subpolar gyre, where the diabatic flux parameterization is unchanged and the squared layer thicknesses differ only slightly relative to the large negative value of , is very similar in the two calculations.

The form in (9) leads (essentially as in Gnanadesikan 1999) to a cubic equation in the approximate solution, which is less convenient for analysis than the original form in (8). However, it is straightforward to include both the weak midlatitude diabatic fluxes and the more rapid high-latitude diabatic exchange time scale by generalizing the original form in (8) of *W* to allow *α _{w}* to take different, constant values for

*y*<

*y*

_{SP}and

*y*<

*y*

_{SP}so that, following (34),

For Δ*y* ≪ 1, this modified model is essentially no more complex than the original model and allows a similar analysis. For the approximate analytical solution, the only resulting changes are to the quantities derived from the area integral of *W* in (21). Equation (28) and the solution in (30) take the same form, but now with

If *α*_{1} ≪ *α*_{2} in (40) and (41), as suggested by the scaling from *κ*, the diabatic fluxes will be dominated by the high-latitude part, proportional to *α*_{2}. In this limit, similar to the case with flux and small *κ*, the small midlatitude diabatic fluxes will have little effect on the adjustment time scale *μ* and eastern boundary depth *h _{Es}*. Thus, with the additional simplification

*α*

_{1}= 0, the approximate solutions (40) and (41) provide an accurate estimate of the time-dependent adjustment to equilibrium for the solution with (Fig. 12). The corresponding modified values (40) and (41) of

*μ*and

*h*are

_{Es}*μ*= 1.0375 and

*h*= 0.67, relative to the values of

_{Es}*μ*= 2.625 and

*h*= 0.92 for the original analytical solution with uniform

_{Es}*α*= 2. This results in the longer, modified time scale

_{w}*T*

_{MOC}= 0.7, corresponding to a dimensional time of more than a century, relative to the original value

*T*

_{MOC}= 0.25, with a dimensional time of only several decades (Figs. 7a, 12). This longer, modified time scale appears to be somewhat shorter than the dimensional adjustment time scales of 10

^{3}yr from the related work by Allison (2009) and L. C. Allison et al. (2010, personal communication). The difference may be the lack of a subpolar-gyre diabatic flux representation in the latter work, the effect of which in the present model would be to reduce

*μ*by setting

*α*

_{2}= 0, with

*α*

_{1}≪ 1, so that the adjustment time scale is primarily dependent on

*α*

_{e}alone.

## 6. Summary

A simple model of the warm-water branch of the middepth meridional overturning cell (Samelson 2009) has been extended here by restoring the local time rate of change of the layer thickness in the mass conservation equation, consistent with the planetary geostrophic approximation for large-scale dynamics, to explore the time-dependent response to changes in external forcing or parameters. Numerical solutions for two adjustment scenarios were considered, first to intensified Northern Hemisphere cooling and then to intensified Southern Hemisphere winds. An analytical solution gives an explicit expression for the thermocline adjustment time scale, which depends on diabatic and eddy flux parameters and on the final steady state. Much of the adjustment of the meridional flux of warm water may occur much faster than this time scale as an essentially instantaneous response to changes in forcing or parameters followed by relatively rapid adjustment through boundary, equatorial, and planetary waves. For the solutions considered here, the dimensional time scale for planetary wave adjustment is several years, whereas that for thermocline adjustment is several decades to a century or more.

The specific adjustment scenarios considered here were chosen as illustrative on general grounds rather than as representative of changes in large-scale ocean forcing that may be anticipated to occur over the next century in association with predicted global climate change. In addition, the anticipated climate-related changes in large-scale forcing are expected to occur on time scales of decades to centuries, comparable to the intrinsic model adjustment time scale so that, for those conditions, it would be more appropriate to analyze the model response to time-dependent forcing parameters rather than the time-dependent adjustment to equilibrium for fixed forcing parameters that is described here. Nonetheless, some aspects of the model response under those changing conditions can be inferred from the existing analysis. General considerations and some global climate models suggest that changes in large-scale ocean forcing are likely to include increased freshwater flux into the northern North Atlantic, as well as related regional surface warming, and increasing Southern Hemisphere westerlies near the latitudes of the Antarctic Circumpolar Current. The former can be represented in the present model as a reduced tendency for Northern Hemisphere cooling, a change opposite of the first adjustment scenario (Figs. 2–4). The latter may be represented by a change in wind forcing similar to the second adjustment scenario (Figs. 8, 9).

The analysis of steady solutions indicates that both of these representative changes for anticipated future forcing conditions would tend to induce in the model an increase in the eastern boundary thickness *h*_{E} and thus a warming and deepening of the midlatitude and tropical thermocline and an increase in southward eddy fluxes of warm water across the circumpolar current. The reduced Northern Hemisphere cooling would cause an increase in *h*_{E} in response to reduced northward warm-water transport into the Northern Hemisphere subpolar gyre, leading to convergence at lower latitudes. This deepening reduces the diabatic upwelling and increases the southward eddy transport across the circumpolar current, resulting in a state with larger *h*_{E} and reduced meridional overturning. The increased Southern Hemisphere winds also induce in the model an increase in *h*_{E}, again representing a warming and deepening of the midlatitude and tropical thermocline. In this case, however, the meridional overturning transport increases at all latitudes (Fig. 9c). Thus, these two anticipated changes in large-scale forcing have opposite effects on the strength of the model midlatitude overturning transport. For the increased Southern Hemisphere winds, the associated increase in meridional overturning transport convergence is spread nearly uniformly across all latitudes so that the increased Southern Hemisphere Ekman transport is balanced primarily by decreased interior upwelling (Fig. 9b). With the weaker and presumably more realistic midlatitude diabatic fluxes described by (34), this midlatitude convergence would be decreased in favor of meridional transport into the subpolar gyre and increased cross-hemisphere overturning. The main short-term changes would be an increase in the meridional convergence of the midlatitude and tropical overturning transport, temporally correlated to the forcing with a lag of only a few years, which would in turn slowly force a deepening of the midlatitude and tropical thermocline on the longer (decades to centuries) adjustment time scale *T*_{MOC}.

The present model and analysis focus entirely on the dynamics of the warm-water branch and its interaction with surface forcing. The interface at the base of the warm layer in the present model can be taken to represent an isothermal surface around 7°–10°C, which outcrops on the equatorward flank of the Antarctic Circumpolar Current and extends well into the North Atlantic subpolar seas. With this model, it is not possible to address questions regarding the interaction of the deep return flow with the warm-water flow or the pathways of the warm-water flow through the warmer, stratified midlatitude gyres and upper thermocline. A more complete model would include a representation of upper-thermocline structure in the midlatitudes and tropics, with a two-layer or multilayer system replacing the present, single, homogenous warm layer. The formulation and analysis of such a model would require attention to the treatment of layer outcrops (Huang and Flierl 1987), which would likely preclude analytical solution and thus substantially limit the conceptual accessibility of the model. However, such an approach would have the advantage of allowing a more complete representation of the dynamics and structure of the midlatitude and tropical components of the warm-water branch of the middepth overturning, as well as the interaction of the warm-water flow with the deeper return flow.

## Acknowledgments

### APPENDIX

#### Frictional Boundary Layers

Some elementary aspects of the boundary layers supported by the model equations in (1)–(3) are briefly summarized here. Related analyses include those by Killworth (1985), Salmon (1986), and Winton (1996). It is assumed that *τ* and *W* may be neglected in the boundary layers.

For the steady problem (*h _{t}* = 0), the equation for the transport streamfunction (e.g., Samelson 2009) has the standard exponential western boundary layer, with decay scale

*r*/

*β*. The no-normal-flow condition at the western boundary leads to a momentum balance resembling that of a Kelvin wave,

with the linear damping term replacing the tangential acceleration. The tangential velocity *υ* may be eliminated to give the condition

which, along with analogous conditions on the other boundaries, is a differential equation that must be solved along the boundary (Salmon 1986); the numerical implementation used here follows Samelson and Vallis (1997), where (A2) is written in terms of the components of **∇***h*^{2} instead of **∇***h*. In the present case, as in Samelson (2009), the condition *hυ* = *V*_{Ek} is imposed in the momentum condition (A1) and (A2) at the southern domain boundary *y* = *y*_{2}, in place of the strict no-normal-flow condition (A2) that would apply at a rigid boundary. [Note that, because the inviscid definition (6) of *V*_{Ek} is used, this induces a small zonal thickness gradient *h _{x}* ~

*r*

^{2}≪ 1 at

*y*=

*y*

_{2}in the numerical solutions.] The sum of

*V*

_{Ek}and the eddy flux

*V*are imposed at

_{e}*y*=

*y*

_{2}as the southern boundary condition for the mass balance equation (3), but

*V*is not included in the equivalent to (A2), because the associated dynamics are not resolved; tests indicate that the solutions are not sensitive to these details of the southern boundary momentum equation.

_{e}The boundary differential equation, (A2) and its analogs, is coupled to the interior through the normal derivative of *h* and the constraint arising from an integration around the boundary of (A2) and its analogs, divided by *f*, in which the integral of the tangential derivative of *h* vanishes exactly. In the case of steady flow (e.g., Samelson 2009), the interior equations and the boundary equation can be solved simultaneously through an iterative process. In the case of unsteady flow, the combined system can support propagating boundary disturbances analogous to those discussed by Killworth (1985) and Winton (1996); these propagating disturbances are not Kelvin-type waves, which as noted above appear only in steady, damped form, but involve instead a combination of frictional and geostrophic-advective effects.

The local structure of the western boundary solution for *h* can be deduced from (A1), (A2), and the exponential boundary layer solution for the transport streamfunction. The latter translates into a similar exponential boundary layer for *h*^{2} so that

where *h _{B}* and

*h*are the local boundary and interior values of

_{I}*h*, respectively, which may vary with

*y*on scales much larger than

*r*/

*β*. With (A2) and

*f*=

*βy*, this means that

which can be integrated directly from an arbitrary point *y* = *y*_{*} to *y*, giving

Thus, will be approximately linear in *y* in regions where ∂*h _{I}*/∂

*y*≈ 0. Near the equator (i.e., for |

*y*| ≪

*y*=

_{eq}*r*/

_{eq}*β*; recall that

*r*≈ 0.25 for |

*y*| ≤ 0.1 in the numerical solutions), the momentum balance may be taken to be ageostrophic,

leading in the steady case to a potential flow ∇^{2}*h* = 0 and, in view of the western boundary condition *u* = 0, to penetration of the boundary value of *h* into the equatorial interior. Alternatively, an approximate balance resembling that of a damped equatorial Kelvin wave follows if the zonal flows are presumed sufficiently strong to be balanced geostrophically. Thickness advection may become significant in either the western or equatorial boundary layers if the tangential velocities become sufficiently large.