A diapycnal diffusion algorithm for isopycnal ocean circulation models with special application to mixed layers

Ocean Modelling 5 (2003) 297–323 www.elsevier.com/locate/omodol

A diapycnal diffusion algorithm for isopycnal ocean circulation models with special application to mixed layers Roland A. de Szoeke *, Scott R. Springer College of Oceanic and Atmospheric Sciences, Oregon State University, 104 Oceanography Administrative Building, Corvallis, OR 97331-5503, USA Received 17 January 2002; received in revised form 14 August 2002; accepted 14 August 2002

Abstract An algorithm is presented for solving the one-dimensional diffusion equation for density, written in terms of density (or a like surrogate) as the independent variable. The algorithm maintains nonnegative layer thicknesses, the premise of the transformation to density as the independent coordinate, under certain restrictions. Near-zero thickness layers can be maintained at the boundaries to accommodate future inflation in response to heating from the boundary. Layers can shrink to near-zero thickness in response to cooling from the boundary. A slight modification of the algorithm permits layers to have diffusion coefficients which differ by orders of magnitude. This provides a natural framework for a surface mixed layer in an isopycnal model, in which the mixed layer is distinguished as a zone of very high turbulent diffusivity overlying an ocean interior of much smaller turbulent diffusivity. The ‘‘mixed layer’’ may be an aggregation of several isopycnal layers rather than just one. A substantial jump in density at the mixed layer base can be represented by several near-zero thickness isopycnal layers. The specification of the thickness of the mixing zone, i.e., the mixed layer depth, is external to the algorithm. An illustration is given using a Kraus–Turnertype specification. Ó 2002 Elsevier Science Ltd. All rights reserved.

1. Introduction The use of isopycnal, isentropic, and hybrid coordinate models of the oceans and atmosphere is becoming common (Bleck, 1974; Bleck and Smith, 1990; Arakawa and Hsu, 1990; Oberhuber, 1993; Hallberg and Rhines, 1996; Bleck, 2001). In such models a thermodynamic variable, such as

*

Corresponding author. E-mail address: [email protected] (R.A. de Szoeke).

1463-5003/03/$ - see front matter Ó 2002 Elsevier Science Ltd. All rights reserved. doi:10.1016/S1463-5003(02)00041-0

298

R.A. de Szoeke, S.R. Springer / Ocean Modelling 5 (2003) 297–323

potential density or potential temperature, is used as the vertical coordinate in place of the geometric height. Other variables than these two have also been suggested (Sun et al., 1999; de Szoeke, 2000; de Szoeke et al., 2000). Whatever its precise definition, we shall generically term the thermodynamic coordinate ‘‘density’’ and models discretized in it ‘‘isopycnal’’ models. The use of isopycnal models is attractive if the variable is nearly conservative, one requirement for which is that irreversible processes such as turbulent mixing are small in some sense. In that case the surfaces on which density is constant describe layers which exchange little mass from one to another, so that a model based on density as coordinate shows very little diapycnal flow compared to the vertical motion in a conventional level-surface model. However, little is not the same as none, and it is still important to take account quantitatively of cross-isopycnal mixing processes. The representation of diapycnal mixing in isopycnal coordinates poses a number of practical problems for numerical models. When density is used as the independent coordinate, the conventional linear Fickian diffusion equation for density transforms into a nonlinear equation. The layer thickness, or vertical spacing of isopycnals, becomes the dependent variable and occurs in the denominator of the terms representing turbulent diffusion. In the continuity equation the effect of density diffusion is to make the evolution of the thickness inversely proportional to the thickness itself. If layer thicknesses become vanishingly small, then very large estimates of the change in layer thickness result. As a consequence, a simple forward scheme can become unstable in a region of a sharp vertical density gradient. Hu (1996) limited diapycnal fluxes in an explicit scheme by replacing the grid point estimates for tracer gradient with vertically averaged values, and Dewar (2001) set the flux divergence to zero for very thin layers. Such artificial measures can be avoided with an implicit scheme, which may, however, require solution of a system of nonlinear equations. Oberhuber (1993) avoided this difficulty by assuming that the entrainment rate in each layer was independent of the others, effectively decoupling the layers. Hallberg (2000) developed an iterative technique for solving a fully coupled implicit formulation of diapycnal diffusion and showed how it could be used to represent Richardson-number-dependent mixing. The numerical scheme must also be able to accommodate a very wide range of diffusivities varying in time and space. Oceanic diffusivities range from perhaps 10
2 m2 s
1 in the surface mixed layer to 10
5 m2 s
1 in the deep ocean. By conceptualizing the mixed layer as a region of enhanced mixing near the surface, a scheme which accommodates such a wide range of diffusivities can be used to represent the full depth of the ocean in a unified framework. Alternatively, mixed layers in isopycnal models have been represented as homogeneous layers with a continuously variable density (Bleck et al., 1992). But these layers do not always match naturally with the underlying isopycnal model, making the detailed implementation awkward. The difficulty of including mixed layers in pure isopycnal models has motivated the use of hybrid isopycnal-pressure coordinates instead (Bleck, 2001). Recently, Dewar (2001) described an isopycnal coordinate implementation of a surface mixed layer model which can be forced by surface fluxes of heat and freshwater as well as penetrative solar radiation. The goal of the present work is to present a new algorithm for solving the density diffusion equation in isopycnal coordinates which is suitable for use in oceanic general circulation modeling. An explicit method for solving the equation is described in Section 2. In the interior layers this technique is sign-definite and stable for any size time step when the diffusivity coefficient is constant. Special methods applied to the boundary layer at the surface and bottom ensure that layer thicknesses do not vanish or become negative. Allowing for a vertically variable diffusivity

R.A. de Szoeke, S.R. Springer / Ocean Modelling 5 (2003) 297–323

299

coefficient requires a simple modification, described in Section 3, involving iteration to prevent negative layer thicknesses. The modified algorithm is demonstrated in an implementation of a mixed layer in a one-dimensional isopycnal coordinate model. The depth of the mixed layer and the size of the eddy diffusivity are specified by a supplementary model, such as the mixed layer model of Kraus and Turner (1967). The flexibility of this approach suggests that it will be useful in three-dimensional oceanic general circulation modeling, which is discussed briefly in the summary, Section 4. This work concentrates on the effects of diapycnal diffusion, which is the sole source of diapycnal mass flux in an ocean model with a linear equation of state. If a nonlinear equation of state is used, the relative variations of temperature, salinity, and pressure on an isopycnal surface give rise to nonlinear effects such as cabbeling and thermobaricity (McDougall, 1987; Davis, 1994; McDougall and Dewar, 1998). These nonlinear effects result in diapycnal velocities which cannot be represented as the Fickian diffusion of density. Mass fluxes arising from these nonlinear effects will not be addressed here. McDougall and Dewar (1998) and Hallberg (2000) discussed strategies for dealing with these difficulties in mixing parameterizations in layered ocean models.

2. Development of the algorithm The one-dimensional diffusion equation, written in terms of potential specific volume (inverse potential density) as the independent variable, leads to
 o o2 K
q ; ð1aÞ pa ¼
2 ot oa pa subject to flux boundary conditions K
q ¼
qtop ;
qbottom pa

ð1bÞ

at a ¼ atop ; abottom . (See Appendix A for a concise derivation and the definitions of parameters.) The solution of Eqs. (1a) and (1b) is invariant to the addition of an arbitrary function q0 ðtÞ to all of qða; tÞ, qtop ðtÞ, qbottom ðtÞ. This invariance represents the (non)effect of a radiative flux passing through the medium without absorption. 2.1. Density discretization Suppose that the independent variable, a, is discretized into L layers between a1 ¼ atop and aL ¼ abottom ; then Eq. (1a) may be written 

 o 1 Kl
1 Dal
1 Kl Dal Dpl ¼
þ ql
1
þ ql ot Dal
12 Dpl
1 Dpl 

 1 Kl Dal Klþ1 Dalþ1 þ ql
þ qlþ1 ; ð2Þ þ Dalþ12 Dpl Dplþ1

300

R.A. de Szoeke, S.R. Springer / Ocean Modelling 5 (2003) 297–323

where variables with integer subscripts are defined in the layers and variables with noninteger subscripts are defined at the interfaces between layers. Layer thicknesses and specific volume intervals are defined so that all Dpl > 0 and Dal > 0. With ql  0, Eq. (2) is equivalent to the onedimensional form of the continuity equations considered by McDougall and Dewar (1998) and Hallberg (2000) (with minor differences arising from the fact that these authors made the Boussinesq approximation and discretized in density rather than specific volume). If, for simplicity, a uniform specific volume spacing, Da ¼ ða1
aL Þ=ðL
1Þ > 0, is assumed then Eq. (1a) is o Kl
1 2Kl Klþ1 Dpl ¼
þ

Ql
1 þ 2Ql
Qlþ1 ot Dpl
1 Dpl Dplþ1

ð3aÞ

for 2 6 l 6 L
1. (Take K1 ¼ 0 and KL ¼ 0.) For l ¼ 1, L special forms of (3a) take account of boundary conditions (1b): o K2 þ Q1
Q2 ; Dp1 ¼
ot Dp2 o KL
1 DpL ¼

QL
1 þ QL : ot DpL
1

ð3bÞ ð3cÞ

The Ql for 2 6 l 6 L
1 are discrete analogues of qðpÞ=Da, while Q1 , QL are qtop =Da, qbottom =Da. Eqs. (3) may be written in a neater form by defining the entrainment mass flux, elþ12 =g (in kg m
2 s
1 ) across the interface between layers l, l þ 1: elþ12 ¼

Kl Klþ1 þ Ql

Qlþ1 Dpl Dplþ1

ð4Þ

for 1 6 l 6 L
1, and e12 ¼ eLþ12 ¼ 0. Sometimes the so-called ‘‘dual entrainment’’ hypothesis (Oberhuber, 1993; McDougall and Dewar, 1998; Hallberg, 2000) is used to interpret or motivate the contributions to (4). According to this view the l-subscripted terms of (4) are an upward mass flux, and the ðl þ 1Þ-subscripted terms are a downward mass flux across the interface. For the purposes here it is important to keep the dependence of the diapycnal fluxes on Dpl at the forefront, so the substitution of an entrainment mass flux is not used extensively. The system of equations represented by (3) has two integral conservation properties. When all these equations are added for 1 6 l 6 L, the right-hand sides cancel exactly, ensuring that mass is conserved, L o X Dpl ¼ 0: ot l

ð5Þ

By weighting Eqs. (3) with their respective specific volumes, and summing over all layers, it also may be shown that the total heat content of the system changes only as a result of the difference between qtop and qbottom : L o X al Dpl ¼ qtop
qbottom : ot l¼1

ð6Þ

Desirable features of a discrete time scheme for (3) are that it preserve the properties (5) and (6).

R.A. de Szoeke, S.R. Springer / Ocean Modelling 5 (2003) 297–323

301

2.2. Temporal discretization The premise of writing the system (3) is that layer thicknesses Dpl remain nonnegative. A scheme like   Kl
1 2Kl Klþ1 ð7aÞ Dplnþ1 ¼ Dpln þ Dt
 þ 

Ql
1 þ 2Ql
Qlþ1 Dpl
1 Dpl Dplþ1 for 2 6 l 6 L
1, and Dp1nþ1 DpLnþ1

 K2 ¼ þ Dt
 þ Q1
Q2 ; Dp2   KL
1 n ¼ DpL þ Dt

QL
1 þ QL DpL
1 Dp1n



ð7bÞ ð7cÞ

for l ¼ 1; L, satisfies the time discrete versions of (5) and (6). The explicit scheme obtained by setting Dpl ¼ Dpln does not guarantee nonnegative Dplnþ1 . An implicit scheme with Dpl ¼ Dplnþ1 seems to require iteration to solve it (Hallberg, 2000). The following two-step scheme, intermediate between these extremes, is therefore proposed: Step 1: 1 1  n  n
Dpl
1 Þ þ ðDplþ1
Dplþ1 Þ Dpl ¼ Dpln þ ðDpl
1 2 2   Kl
1 2Kl Klþ1 þ Dt
 þ 

Ql
1 þ 2Ql
Qlþ1 Dpl
1 Dpl Dplþ1

ð8aÞ

for 2 6 l 6 L
1, and

  1 1 K1 K2  n  n ¼ þ ðDp2
Dp2 Þ þ ðDp1
Dp1 Þ þ Dt þ Q1
Q2 ;
2 2 Dp1 Dp2   1 1 KL
1 KL  n  n  n DpL ¼ DpL þ ðDpL
1
DpL
1 Þ þ ðDpL
DpL Þ þ Dt
 þ 
QL
1 þ QL 2 2 DpL
1 DpL

Dp1

Dp1n

ð8bÞ ð8cÞ

for l ¼ 1; L (remember K1 ¼ KL ¼ 0); Step 2: n  n  Dplnþ1 ¼ Dpl þ 12ðDpl
1
Dpl
1 Þ þ 12ðDplþ1
Dplþ1 Þ

ð9aÞ

for 2 6 l 6 L
1, and Dp1nþ1 ¼ Dp1 þ 12ðDp2n
Dp2 Þ þ 12ðDp1n
Dp1 Þ; DpLnþ1

¼

DpL

þ

n 1 ðDpL
1 2




 DpL
1 Þ

þ

1 ðDpLn 2




DpL Þ

ð9bÞ ð9cÞ

for l ¼ 1; L. By adding together the corresponding equations in these two sets, (8) and (9), it is easily seen that the set (7) is recovered. By rearranging Eqs. (8), we see that they are solved by Dpl


2Kl Dt ¼ Dpln þ 2Ql Dt  xl Dpl

ð10Þ

302

R.A. de Szoeke, S.R. Springer / Ocean Modelling 5 (2003) 297–323

for 1 6 l 6 L. This is a quadratic equation for each Dpl , whose solution is Dpl ¼ 12½xl þ ðx2l þ 8Kl DtÞ1=2 :

ð11Þ

The positive branch of the square root in (11) is chosen because only this ensures that Dpl P 0. (This solution is not unique: an arbitrary constant, perhaps different at each time n, may be added to the right of (10). This nonuniqueness is linked to the invariance of Eqs. (7) to the addition of a constant to all Ql .) Eq. (11) is used in the second step, (9), to produce the Dplnþ1 . 2.3. Boundary layers Na€ıvely implemented, the equation for Dp1nþ1 permits negative values. To forestall this, two devices are proposed. First, the total boundary (either top or bottom) buoyancy flux, Q0 , is assumed to be absorbed linearly with depth: ð12aÞ Ql ¼ Q0 q0 ðpl
12 Þ; where q0 ðpÞ

¼1


p ; pQ

¼ 0;

p < pQ ;

) ð12bÞ

p P pQ

with plþ12 ¼

l X

Dpjnþ1 ;

p12 ¼ 0:

ð13Þ

j¼1

It is important that thicknesses at time-level n þ 1 be used in (12) and (13), and thus in (10) and (11), even though they are not yet known. How this implicit dependence of the amount of heat absorption on the anticipated layer thicknesses is taken into account will become clear in a moment. The effect of (12) is that the amount of downward buoyancy flux Q0 absorbed by a layer within a distance pQ of the surface is proportional to the thickness of that layer. A vanishingly thin layer absorbs no radiation but transmits it to the next layer. Second, the turbulent exchange coefficient grows with distance from the surface according to Kl ¼ K0 k 0 ðpl
12 Þ; where k 0 ðpÞ

) p ; p 6 pK ; pK ¼ 1; p > pK :

¼

ð14aÞ

ð14bÞ

Thus, the Kl of layers within pK of the surface depend implicitly on the thickness, Dplnþ1 . They grow linearly to a value K0 beyond pK of the boundary (chosen to be similar to pQ ). Classically, K  kv u pq0 g in a turbulent shear flow along a wall, where u is the friction velocity and kv is K arm anÕs constant. This provides guidance for the choice of K0 and pK (Large et al., 1994). To see the effect of the implicit dependencies of Ql and Kl , suppose that the first layer is thinner than pQ or pK , i.e., that Dp1nþ1 < minðpQ ; pK Þ. In this case, the first choices for q0 and k 0 in (12b) and (14b) apply. Explicitly writing down x2 and Dp2 from the right-hand side of (10) and (11) gives

R.A. de Szoeke, S.R. Springer / Ocean Modelling 5 (2003) 297–323

x2 ¼ Dp2n
2Q0 and

Dp1nþ1 Dt þ 2Q0 Dt pQ

( 1=2 )
nþ1 1 Dp x2 þ x22 þ 8K0 1 Dt Dp2 ¼ : 2 pK

Combining these with the left-hand side of (10) gives ( 1=2 )
K2 Dt 1 Dp1nþ1 2 x2
x2 þ 8K0 : ¼
Dt Dp2 4 pK Now substitute (12), (14) and (17) into (7b): ( 1=2 )
nþ1 nþ1 Dp 1 Dp ; Dp1nþ1 ¼ Dp1n þ Q0 Dt 1 þ x2
x22 þ 8K0 1 Dt 4 pQ pK or, substituting x2 , which is linear in Dp1nþ1 , from (15), 

1=2 1 Q0 Dt 1 1 1 2 Dpnþ1 Dp1nþ1 þ Dp1n þ Dp2n þ Q0 Dt ¼
1
: x2 þ 8K0 1 Dt 2 pQ 4 2 4 pK

303

ð15Þ

ð16Þ

ð17Þ

ð18Þ

ð19Þ

Squaring both sides of this yields a quadratic equation for Dp1nþ1 , whose solution can be written .  1=2 a; Dp1nþ1 ¼ b
ðb2
acÞ where a ¼ 1
Q0 Dt=pQ ; b ¼ 12ð1 þ aÞDp1n þ 14ðDp2n þ 2Q0 DtÞa þ K0 Dt=ð4pK Þ; c ¼ Dp1n ðDp1n þ 12Dp2n þ Q0 DtÞ:

ð20Þ

The properties of this algorithm, and in particular the fairly mild conditions which guarantee positive Dp1nþ1 , are discussed in Appendix B. The negative sign of the radical in (20) must be taken to ensure that Dp2 (Eq. (16)) be positive. The algorithm (20) replaces (10) (or (11)), for l ¼ 1, and (9b) for the determination of Dp1nþ1 . Should the next layer also be contained within pQ or pK of the surface, i.e., Dp1nþ1 þ Dp2nþ1 < minðpQ ; pK Þ, then Dp2nþ1 is obtained in a similar way, for x3 and Dp3 depend on Dp2nþ1 through Q3 and K3 in a manner similar to (15) and (16). Thus an expression for Dp2nþ1 similar to (20) can be developed. The procedure can be extended to obtain any number of layer thicknesses contained within pQ or pK . 2.4. Positiveness of interior layer thicknesses A criterion exists which, under certain conditions, guarantees that Dplnþ1 P 0, for interior layers. Together with the integral constraints (5) and (6) this criterion ensures the stability of the algorithm posed by (9)–(11).

304

R.A. de Szoeke, S.R. Springer / Ocean Modelling 5 (2003) 297–323

Using (10) and (11) one may show that 1 1 Kl Dt
Ql Dt: ðDpln
Dpl Þ ¼ fxl
4Ql Dt
ðx2l þ 8Kl DtÞ1=2 g ¼
2 4 Dpl

ð21Þ

Substituting this into (9a), one obtains     Dplnþ1 Kl
1 Dt Ql
1 Dt Klþ1 Dt Qlþ1 Dt ¼1
þ
þ :   Dpl
1 Dpl Dpl Dplþ1 Dpl Dpl Dpl

ð22Þ

For layers completely beyond pQ of the boundary, Ql
1 ¼ Qlþ1 ¼ 0. Eq. (11) may be written ! Dpl xl ; ¼f ð2Kl DtÞ1=2 ð2Kl DtÞ1=2

ð23aÞ

where f ðfÞ ¼ 12½f þ ðf2 þ 4Þ1=2 :

ð23bÞ

For f P 0, f ðfÞ P 1 þ 12f > 1. Hence Dpl P 12xl þ ð2Kl DtÞ1=2 > ð2Kl DtÞ1=2 ;

ð24Þ

and it follows that Dplnþ1 Kl
1 Dt Klþ1 Dt 1 ¼1

  P1
  Dpl
1 Dpl Dpl Dplþ1 2 Dpl




Kl
1 Kl

1=2

1
2




Klþ1 Kl

1=2 :

ð25Þ

If all Kl are the same, the right side of this inequality is zero. Thus inequality (24) assures the positiveness of Dplnþ1 for interior layers with uniform diffusivities, Kl ¼ constant. Exceptional layers contained wholly within pQ or pK of the boundary, or straddling this transition region and the interior, are dealt with in Appendix B. The case of variable diffusivities will be considered in the discussion of the mixed layer implementation. 2.5. Examples Two classic diffusion problems demonstrate the features of the algorithm just developed. These solutions are shown, not merely to demonstrate that the algorithm well reproduces classical diffusion problems, but to indicate how inflation and shrinkage of specific volume layers, to and from nearly zero thickness, is employed to represent the solutions. 2.5.1. Diffusion of a discontinuity The solution of the one-dimensional diffusion equation in conventional pressure coordinates, oa o2 a ¼K 2 ot op

ð26Þ

R.A. de Szoeke, S.R. Springer / Ocean Modelling 5 (2003) 297–323

with insulating boundary conditions oa ¼ 0 at p ¼ 0; pb ; K op and a step function initial profile 1 a ¼ atop ; 0 < p 6 pb ; 2 1 ¼ abottom ; pb < p < pb 2 is " 1 a X
a bottom top aðp; tÞ ¼ atop þ ð
1Þn erf 1þ 2 n¼
1

305

ð27Þ

ð28Þ

ðp
12pb Þ þ npb ð4KtÞ1=2

!# :

ð29Þ

The asymptotic steady state is aðp; 1Þ ¼ ðabottom þ atop Þ=2. To represent this initial value problem in the numerical model, the interval between atop and abottom was divided into L ¼ 32 layers, such that initially Dp1 ¼ DpL ¼ 12 pb ¼ 0:5 MPa (50 m), while Dpl ¼ 1  10
2 Pa (practically zero thickness) for 1 < l < L. All Ql were set to zero for the duration of the experiment, and the turbulent diffusivity was K ¼ 104 Pa2 s
1 (10
4 m2 s
1 ). The analytic and discrete solutions obtained using a time step of 100 s at several times are shown in Fig. 1a. As the experiment progresses, layers of initially near-zero thickness at mid-level and intermediate a inflate spontaneously while layers representing ranges of a toward the extremities, atop ; abottom , gradually shrink to zero at the boundaries p ¼ 0; pb . Differences between the solutions are most pronounced in the early stages and near the center of the domain. The closeness of the discrete solutions to the analytic solution can be seen better by comparing the slopes (or inverse buoyancy fluxes) of the two solutions at the midpoint in pressure (Fig. 1b). Also shown is the effect of using a different Dt. If the time step is chosen to be too large, then the initial heat flux takes too long to get established, and, in the long term, diffusion proceeds too slowly, as if the diffusion coefficient were effectively reduced. Both empirical studies and theoretical considerations indicate that convergence is first order in time. Convergence also depends on the number of layers, which is represented in the scaling of the ordinate of Fig. 1b, though all experiments presented here have the same number of layers. 2.5.2. Cavity heated from above, cooled from below The heat conduction problem considered above was modified to represent the boundary conditions oa ¼ q at p ¼ 0; pb ð30Þ K op with the initial condition ðatop þ abottom Þ at t ¼ 0: ð31Þ a¼ 2 The asymptotic steady state is the linear profile  a þ a q
1 top bottom ð32Þ p
pb : þ a¼ K 2 2

306

R.A. de Szoeke, S.R. Springer / Ocean Modelling 5 (2003) 297–323

Fig. 1. (a) The solid line represents the discrete solution obtained using the algorithm described in the text and the dashed line is the exact solution obtained from (29). The initial distribution of the specific volume is 0:84  10
6 m3 kg
1 for p 6 0:5 MPa and 1  10
6 m3 kg
1 for 0:5 < p 6 1 MPa. After t ¼ 0 the specific volume is allowed to diffuse. A sequence of intermediate steps (t ¼ 5002 ; 10002 ; 15002 ; 20002 s) is represented. In the final state (not shown) the specific volume is uniformly distributed throughout the domain with a value equal to the average value at the start, 9:2  10
6 m3 kg
1 . (b) The layer thickness (which is inversely proportional to the specific volume flux) at the center of the domain as a function of time. The results for three different sized time steps (Dt ¼ 10; 100; 1000 s) as well as the analytic solution are shown. The numerical solution converges to the analytical solution as the time step shrinks. Convergence depends also on the number of layers, which is represented by the scaling of the ordinate, although all examples shown used L ¼ 32.

The initial condition was represented in the numerical model by setting DpL=2 ¼ pb ¼ 1 MPa for the layer representing a ¼ ðatop þ abottom Þ=2, and all other Dpl ¼ 1  10
2 Pa. The boundary heat fluxes were Q1 Da ¼ QL Da ¼ q ¼ 1  10
8 m2 s
3 , while all other Ql were zero. The increment Da

R.A. de Szoeke, S.R. Springer / Ocean Modelling 5 (2003) 297–323

307

Fig. 2. The initial condition is a uniform value of a ¼ 9:2  10
6 m3 kg
1 throughout the domain. Beginning at t ¼ 0 buoyancy fluxes q ¼ 1  10
8 m2 s
3 are applied at p ¼ 0 and at pb ¼ 1 MPa. With the passage of time, layers which were initally collapsed at the boundaries inflate and grow toward the middle of the domain. A series of intermediate states at times t ¼ ð10002 ; 20002 ; 30002 ; 60002 Þ s shows the eventual development of a uniform density gradient. In the final state the gradient allows for a constant buoyancy flux from one boundary to the other. Note that some layers near the boundaries remain collapsed.

was chosen to be at least qpb =ðKðL
1ÞÞ (for L ¼ 32), so that ðL
1ÞDa accommodates the expected range of the solution. Fig. 2 shows the numerical solution at successive times for Dt ¼ 1000 s. (The solution was not as sensitive to the size of the time step as in the previous experiment.) Initially massless layers near the boundaries inflate sequentially to represent the solution. If ðL
1ÞDa spans a wider range than necessary to represent the solution, as in Fig. 2, the excess layers are never inflated. If q is doubled (Fig. 3) and the other parameters are held constant, then an insufficient range of specific volumes is available to represent the solution, i.e., ðL
1ÞDa < qpb =K. In that case the correct linear profile is approached in the center of the cavity, but it is bracketed by excessively thick end layers for l ¼ 1; L. The heat flux short-circuits across these layers (for which K1 ¼ KL ¼ 0) to establish the linear gradient over the portion of the cavity where it can be represented. These experiments demonstrate that it is important to provide a sufficient range of specific volume layers to represent the final solution, but that if that is not done (perhaps because the solution is unknown), the algorithm does not fail.

3. One-dimensional mixed layer modeling The framework that has been developed for representing diffusion between discrete density layers may be adapted for modeling mixed layers. Mixed layers at the ocean surface are regions of

308

R.A. de Szoeke, S.R. Springer / Ocean Modelling 5 (2003) 297–323

Fig. 3. The experiment depicted in Fig. 2 was repeated with twice the buoyancy flux applied at the boundaries. In this case the span of densities specified was inadequate to represent the linear gradient required at steady state. The buoyancy flux ‘‘short circuits’’ across the excessively thick end layers to establish the requisite linear gradient in the middle of the domain.

intensive turbulence maintained by the agitation of the wind or convection due to cooling. These layers are usually very homogeneous in their properties and often marked by an abrupt transition to a much more quiescent ocean interior. The feature of mixed layers germane to the layer representation of diffusion is that  K ¼ KML ; for p 6 pML ; ð33Þ ¼ KI ; for p > pML ; where KML  KI . The thickness of the mixed layer, pML , may be a function of time, and may be thought of as controlled by external parameters such as wind stress and surface cooling and internal parameters such as ocean stratification and current shear. A range of semi-empirical theories is available for specifying pML . It is not necessary for the present purposes to canvass any particular choice, but merely to suppose that pML is given. Similarly, the diffusion coefficients KML and KI will be supposed given. Detailed mixed layer theories aspire to predict these parameters, which for simplicity we take to be constants, though vertical variation is possible. The noteworthy feature of these parameters is that in the mixed layer KML is very large––perhaps 10
3 to 10
2 m2 s
1 ––while in the interior KI is orders of magnitude smaller––perhaps 10
6 to 10
4 m2 s
1 . 3.1. Variable diffusivity profile Let us examine the consequences of such choices of KML , KI in the algorithm outlined above. Let l1 be the index of the layer below that containing pML (which is assumed to be bigger than pQ

R.A. de Szoeke, S.R. Springer / Ocean Modelling 5 (2003) 297–323

309

Table 1 Layer thickness criterion around the mixed layer base Layer index

Kl

RHS inequality (25)

l < l1
1 l ¼ l1
1

KML KML

0 1 ½1
ðKI =KML Þ1=2 > 0 2

l ¼ l1 l P l1 þ 1

KI KI

1 ½1 2


ðKML =KI Þ1=2 < 0

0

or pK ). That is, Kl ¼ KI for l P l1 , while Kl ¼ KML for l < l1 . (For layers contained within pQ or pK , set K0 ¼ KML in Eqs. (14) and what follows.) The condition that a layer thickness remain nonnegative during a time step is that the right-hand side of inequality (25) is P 0. The situation for the layers around l1 is summarized in Table 1. The condition is satisfied for all layers except by l1 itself. A moments reflection shows that it is reasonable to expect layer l1 to be losing mass, even to the vanishing point, to l1
1, the layer just above it that constitutes the mixed layer. To avert negative Dpl1 the following simple stratagem is proposed. If, in a given time step, Dpl1 becomes negative, that step is repeated with Kl1 reset to KML . This has the effect of incrementing l1 in the discussion just given. Then the concern becomes whether layer l1 þ 1 disappears in the given time step. But if so l1 is incremented yet again. Eventually a layer will be found whose thickness does not become negative (excluding only the possibility where the mixed layer occupies the entire water column), and the time step will be successfully advanced. The final index l1 may differ somewhat from the initial l1 . On the next time step the initial index l1 is again chosen to be the layer below that containing pML . This procedure guarantees positive layer thicknesses. More generally, a high-diffusivity layer grows at the expense of its low-diffusivity neighbors. This is why condition (25) guarantees positivity of layer thicknesses only when diffusivity is uniform among neighbors. One hardly expects a contrast of diffusivity anywhere else as large as across the mixed layer base. But a straightforward elaboration of the strategy outlined will ensure positive thicknesses everywhere when diffusivity varies between layers. 3.2. Bulk mixed layer depth A method for determining the mixed layer thickness pML used in (33) to specify the diffusivity profile will now be described. It is based on the model of Kraus and Turner (1967) for the effects of wind stirring and buoyant forcing on the ocean surface layers. The Kraus–Turner mixed layer parameterization has been used in isopycnal ocean models (Bleck et al., 1989), though not in the way proposed here. Modifications for the effects of energy sources for mixing from surface current shears can be made (Pollard et al., 1973; Niiler, 1975; Niiler and Kraus, 1977; de Szoeke and Rhines, 1976; Turner, 1981). The Kraus–Turner equation for the deepening of the mixed layer is 1
1 o 1 3
1 ð34Þ g dapML pML ¼ a
1 0 m0 u
pML g qtop nQ : 2 ot 2 The terms in this equation have units of W m
2 . The term on the left is the rate of change of the vertically integrated potential energy of the surface layers of the ocean as the mass field is redistributed by mixing; pML is the thickness of the mixed layer in pressure units (Pascals: multiply by a0 g
1 to get meters); da ¼ aML
aðpML Þ is the specific volume difference between the mixed

310

R.A. de Szoeke, S.R. Springer / Ocean Modelling 5 (2003) 297–323

Fig. 4. Schematic drawing of Kraus–Turner bulk mixed layer model and discrete density mixed layer model. In the latter, the ‘‘mixed layer’’ may be weakly stratified, represented by more than one discrete specific volume layer. In this region the high diffusion coefficient of the mixed layer is applied. Beneath the mixed layer there may lie a thin transition layer, which has the weaker, intermediate diffusion coefficients but relatively high buoyancy flux because of the thinness of the layers. The specific volume jump analogous to that used in the Kraus–Turner model is calculated between the mean specific volume of the mixed layer and the specific volume of the layer just beneath the mixed layer base, where the buoyancy fluxes diminish to values typical of the deep ocean.

layer aML and the fluid below aðpML Þ (see Fig. 4). The first term on the right is a parameterization of the rate of mechanical stirring by the surface wind stress s, where u ¼ ða0 jsjÞ1=2 is the friction velocity, and m0 is a dimensionless parameter taken to be 1.0 (Davis et al., 1981). The second term is the effect of surface buoyancy flux qtop on mixed layer deepening: downward flux of positive buoyancy, qtop > 0, tends to lower mixed layer potential energy, and the opposite sense, qtop < 0, tends to increase it. The parameter nQ is taken to be 1 when qtop > 0 and a smaller value when qtop < 0, to account for dissipation of buoyant energy production when this is destabilizing. Though nQ  0:2 is the customary choice based on laboratory experiments, we shall use nQ ¼ 0:5 to produce a more emphatic density step at the mixed layer base (see below). When qtop > 0, the right side of (34) is physically meaningful only as long as it is positive. Otherwise pML is determined to be the value which makes the right side zero, 2m0 u3 g pML ¼ ; ð35Þ a0 qtop this is (proportional to) the Monin–Obukhov depth. The conservation of buoyancy requires that the integral from p ¼ 0 to pML of the difference between aML at time t and the initial aðpÞ must be supplied by the accumulation of surface buoyancy flux, i.e., Z t Z pML faML
aðpÞgdp ¼ qtop dt: ð36Þ 0

0

R.A. de Szoeke, S.R. Springer / Ocean Modelling 5 (2003) 297–323

311

For the choice aðpÞ ¼ a0 þ a0 p, this gives da ¼ aML
aðpML Þ ¼

qtop t 1 0
a pML pML 2

ð37Þ

for the specific volume jump at the base of the mixed layer. Eq. (34) may be discretized in time by means of nþ1 n ¼ pML þ Dt pML

n 2m0 u3
pML a0 g
1 qtop nQ : n n g
1 a0 pML faML
aðpML Þg

On the right we may take X lðpX 6 pML Þ aML ¼ Dpl ; al Dpl

ð38Þ

ð39Þ

l¼1

the sum being taken over the layers down to pML , while aðpML Þ is the al of the layer below the mixed layer base. The mixed layer base is defined as the depth at which the magnitude of the buoyancy flux transitions from the high values of the mixed layer to the lower values of the deep ocean. 3.3. Examples A number of calculations will be shown that illustrate the use of (38) in determining pML in (33). In all the examples the initial state consisted of a linear specific volume gradient giving a proportional specific volume contrast of ða1
aL Þ=a0 ¼ 10
3 over a pressure interval of 1 MPa (100 m), divided equally over L ¼ 32 specific volume layers (so that Da=a0 ¼ ða1
aL Þ=ððL
1ÞaÞ), except as noted below. The diffusivities used in (33) are KML ¼ 1  105 Pa2 s
1 (1  10
3 m2 s
1 ) and KI ¼ 2:5  101 Pa2 s
1 (2:5  10
7 m2 s
1 . The latter is perhaps an unrealistically small value even for the deep ocean, but that does not affect the short-term solutions presented here.) The time step was 1000 s. 3.3.1. Forced convection In this example a wind stress of magnitude u2 ¼ a0 s ¼ 10
4 m2 s
2 commences to act on the sea surface at t ¼ 0, stirring the surface layers, with no surface buoyancy flux, qtop ¼ 0. Fig. 5 shows a time series of profiles constructed from the solutions for the Dpl produced by the algorithm. Large diffusivity defines the ‘‘mixed layer’’ which is given by advancing (38), shown in the middle panel of Fig. 5. The figure illustrates how the mixed layer is represented within the isopycnal-layer model. Note in particular the thin layers at the surface representing light water whose presence is no longer required as time proceeds, and the vanishingly thin layers at the mixed layer base representing the sharp specific volume jump that occurs there. While there is no restriction on the thinness of these layers, the procedures described above guarantee that they remain nonnegative. In the continuous analogue given by (34), the specific volume jump at the base of the mixed layer is da ¼ aML
aðpML Þ ¼
12a0 pML ;

ð40Þ

312

R.A. de Szoeke, S.R. Springer / Ocean Modelling 5 (2003) 297–323

Fig. 5. The response to wind stress started impulsively at t ¼ 0 involves deepening and decreasing the specific volume of the mixed layer. The top panel shows the mixed layer specific volume as a function of time. The dashed line is the analytical solution to the Kraus–Turner equation (41) and the solid line is the solution from the discrete model. The middle panel similarly shows the mixed layer depth as a function of time for which the analytic and discrete solutions are indistinguishable. The bottom panel shows the specific volume profiles at four times (in 106 s), each offset by Da ¼ 1  10
6 m3 kg
1 . Note that the mixed layer is weakly stratified in the discrete model. This accounts for the slight difference between the specific volumes in the discrete and analytical model in the top panel.

where a0 =a0 ¼ 10
3 MPa
1 . (The last equality in (40) follows from the principle that the vertically integrated buoyancy of the mixed layer must not change in the absence of surface buoyancy flux. Hence the mixed layer specific volume aML must be the average of aðpML Þ and að0Þ.) Thus the solution of (34), with initial condition pML ð0Þ ¼ 0, is gu pML ¼ ð41Þ ð12m0 NtÞ1=3 ; a0 N 1=2

where N ¼ gð
a0 Þ =a0 ¼ 2pð5:6 cphÞ is the buoyancy frequency. This analytical solution is shown in Fig. 5, together with the numerical solution of (38), calculated along with numerical solutions for Dpl obtained by the methods described above. In the instantaneous profiles of specific volume one sees a number of features of the densitylayer representation of the mixed layer. The mixed layer itself may be made up of a number of density layers, typically three in the example of Fig. 5. Choosing a larger KML will tend to reduce this number, allowing the central layers to dominate and expelling the marginal layers to the surface and to the mixed layer base. At the surface, as the specific volume of the mixed layer decreases (cooling), more and more lighter layers shrink to vanishing thickness. As mixing proceeds and the mixed layer entrains and advances, layers in the interior at first shrink as they become part of the mixed layer transition, and then grow to nearly fill the mixed layer as their value of specific volume for a time dominates it, and perhaps eventually vanish in the massless

R.A. de Szoeke, S.R. Springer / Ocean Modelling 5 (2003) 297–323

313

stack of layers at the surface. The mixed layer base and transition to the interior, often modeled as a discontinuity, may be seen, given adequate specific volume resolution, as a high but finite gradient region composed of several thin layers. 3.3.2. Natural convection Consider the opposite extreme with mixed layer convection driven solely by surface cooling. For this case we set the buoyancy flux qtop ¼
1  10
7 m2 s
3 ()200 W m
2 ), and u ¼ 0 (no wind). Using (37) in (34), one obtains pML

½ð2 þ 4nQ Þjqtop jt 1=2 : ¼ a0 N =g

The specific volume jump at the base of the mixed layer is nQ da ¼ ð
a0 ÞpML : 1 þ 2nQ

ð42Þ

ð43Þ

These analytical solutions for pML are shown in Fig. 6, for nQ ¼ 0:5, together with the corresponding numerical solution of (38). Also shown are profiles of specific volume obtained from the numerical solutions at several times. Many of the remarks made about Fig. 5 apply here. Nonzero nQ gives penetrative convection, for which a nonzero specific volume jump da occurs at the mixed layer base (Fig. 6); for nQ ¼ 0, or nonpenetrative convection, no such jump occurs (Fig. 7).

Fig. 6. The response to a negative buoyancy flux applied at the surface in the absence of wind is a strong decrease in the specific volume and deepening of the mixed layer. The panels are as described in Fig. 5. As the mixed layer deepens, sequentially denser layers constitute the mixed layer in the discrete model. These layers collapse to near-zero thickness at the oceanÕs surface. This case of penetrative convection, nQ ¼ 0:5, has a clear jump in specific volume at the base of the mixed layer.

314

R.A. de Szoeke, S.R. Springer / Ocean Modelling 5 (2003) 297–323

Fig. 7. As in the lower panel of Fig. 6, but for nonpenetrative convection nQ ¼ 0. No jump in specific volume forms at the base of the mixed layer.

In response to cooling a less dense (warmer) layer near the surface shrinks until it becomes smaller than the boundary layer thickness, pQ , when the boundary layer strategy limits its thickness to slightly greater than zero and the cooling is transferred to the next layer below. Over time, the less dense layers collapse to the surface, exposing more dense layers to the surface buoyancy flux. When the mixed layer reaches the bottom and the densest layer is exposed to the surface, the onedimensional model fails. The cooling process cannot continue indefinitely. 3.3.3. Forced convection with heating A numerical experiment was run with both wind stirring, a0 s ¼ u2 ¼ 10
4 m2 s
2 , and positive buoyancy flux, qtop ¼ 0:5  10
7 m2 s
3 . In order to allow for warming of the mixed layer, 32 collapsed layers were established initially at the surface in addition to the 32 layers used in other experiments. The mixed layer in both the analytic and numerical solutions deepens (Fig. 8) until

Fig. 8. When wind stress and positive buoyancy flux were applied at the surface, the mixed layer tends to the Monin– Obukhov depth and continues to increase in specific volume. The panels are as described for Fig. 5.

R.A. de Szoeke, S.R. Springer / Ocean Modelling 5 (2003) 297–323

315

they asymptotically approach the Monin–Obukhov depth of 0.388 MPa given by (35). With no entrainment cooling the bottom of the mixed layer, the continued positive buoyancy flux at the surface causes the mixed layer to continue warming. This heating process is represented in the model by ever lighter layers, which were vanishingly thin at the surface at t ¼ 0 s, successively inflating and occupying the column down to the Monin–Obukhov depth. If the model were to run out of these layers, or if they were not provided in the first place, then the mixed layer would begin to deepen again, much as in Fig. 3, even though this deepening is not part of the analytic solution. As in the previous case, if the least dense layer reaches the bottom of the domain, then continued heating will cause the algorithm to fail. 3.3.4. Restratification The forced convection case of Fig. 5 was run out to t ¼ 2:5  105 s as described above with qtop ¼ 0, whereupon positive buoyancy flux, qtop ¼ 0:5  10
7 m2 s
3 , was applied. This immediately limited pML to the same Monin–Obukhov depth as in Fig. 8. The water above this level became lighter as this buoyancy flux continued. As in the previous case, ever lighter layers, which were collapsed at the surface before t ¼ 2:5  105 s, successively grew to occupy the water column down to the Monin–Obukhov depth (Fig. 9). Below this depth the pycnostad of the previous mixed layer and the strong pycnocline at its base are preserved . This density jump diffuses away in a process like that shown in Fig. 1a, but only slowly because of the very weak diffusivity below the Monin–Obukhov depth. This process is called restratification, or (somewhat misleadingly) mixed layer retreat, because the pycnocline at the mixed layer base appears to retreat to shallower layers. 3.3.5. Shear production Niiler (1975) modified the derivation of (34) to include effects of turbulent energy production by wind-driven sheared currents at the mixed layer base. Pollard et al. (1973) devised a model of shear-driven deepening based on a bulk Richardson number criterion. De Szoeke and Rhines (1976) showed quantitatively the circumstances under which the Niiler model could reproduce either the Pollard et al. result or the Kraus–Turner deepening. Because the Niiler modification of (34) has features that make it undesirable for numerical implementation––a denominator which approaches zero when shear mixing dominates––an alternative approach was adopted that harks back to the Pollard et al. model. To calculate shear, the low-frequency current must be calculated. This is done by solving onedimensional momentum equations for the system of layers. In these equations, the momentum

Fig. 9. The final state of Fig. 5 was used as an initial condition, and a surface buoyancy flux was applied. The mixed layer retreats to the Monin–Obukhov depth (0.388 MPa) as in Fig. 8. The specific volume gradient of the previous mixed layer base decays slowly owing to the weak diffusivity in the interior ocean.

316

R.A. de Szoeke, S.R. Springer / Ocean Modelling 5 (2003) 297–323

Fig. 10. The configuration of Fig. 5 was repeated with shear mixing enabled. Inertial oscillations initiated by the onset of the wind stress cause shear between the mixed layer and deeper ocean. The effect of this shear mixing is to enhance the deepening of the mixed layer during the first half-inertial period. Two inertial periods are shown.

exchange between layers associated with entrainment mass fluxes, Eq. (4), are taken into account, along with turbulent interfacial stresses. A bulk Richardson number is calculated from Rib ¼

pML ½aML
aðpML Þ

½uML
uðpML Þ

2

;

ð44Þ

where uML is the average mixed layer velocity, calculated from the ul by a formula just like (39), and uðpML Þ is the velocity below the mixed layer. Whenever Rib < 1, pML is incremented to include the next discrete layer. Over several time steps this procedure has the effect of reducing Rib below 1, while pML is driven deeper than countenanced by the Kraus–Turner deepening formula (38). By these means a mixed layer thickness pML is calculated for use in specifying the turbulent diffusivity profile according to (33). For an impulsively started wind stress, inertial motions are generated with magnitude given by jpML uML j ¼ gf
1 js0 j21=2 ð1
cos ftÞ1=2 :

ð45Þ

In the absence of surface buoyancy flux, the specific volume jump is given by (40). Substituting this in (44) where Rib is set to 1, one may solve for pML (assuming uðpML Þ ¼ 0) (Pollard et al., 1973), gu 1=4 21=2 ð1
cos ftÞ : ð46Þ pML ¼ 1=2 a0 ðfN Þ This analytic solution is valid for the first inertial half-period, pf
1 ¼ 3:1  104 s, for which it is increasing, and is shown in Fig. 10. Also shown is the numerical solution for shear-modified pML calculated with the Dpl and ul obtained by the methods described. This solution, though jerky, follows the analytic solution quite well. Beyond t ¼ 3:1  104 s, pML is well ahead of the Kraus– Turner mixed layer thickness (Fig. 5) that it would have assumed but for the shear mixing effects, and does not increase much.

4. Summary and discussion An algorithm has been propounded for solving diapycnal diffusion problems in isopycnal coordinates. The algorithm proceeds in two steps: a preliminary step obtains provisional time advances of isopycnal layer thicknesses, which are used for calculating the diffusion terms; a second

R.A. de Szoeke, S.R. Springer / Ocean Modelling 5 (2003) 297–323

317

step calculates final estimates of layer thicknesses at the next time step. Under fairly mild restrictions on the distribution of interior buoyancy sources, nonnegative layer thicknesses may be guaranteed in the interior layers. The restrictions constitute restraints on the range of physical problems that can be considered: there is no question that buoyancy sources could be arranged in such a way as to form specific volume inversions, i.e., negative layer thicknesses. This would vitiate the assumption underlying the basic transformation to isopycnal coordinates. It is not that such situations are impossible in the oceans, but rather that they are rare. It is because of this rarity, and because of the value of expressing oceanic properties relative to the distribution of specific volume, that it is useful to solve the diffusion equation with specific volume as the independent variable. Layers at the boundary surfaces, especially the top, may become vanishingly thin––for example, in response to negative buoyancy flux. To avert negative layer thicknesses at the boundaries, measures are adopted which short-circuit the surface buoyancy input to deeper layers as the surface layer becomes thin and reduce the turbulent diffusivity in proportion to the distance from the surface. These measures are easily implemented and allow layers at the surface to approach near zero without any other intervention. Conversely, under different forcing conditions (e.g., when surface buoyancy flux is positive), collapsed, extremely thin layers at the surface may inflate as necessary in a natural, self-regulating way. In either case, care must be taken to ensure that the range of specific volumes represented by the layers at the beginning of the experiment is sufficient to represent the solution. It may be possible to enhance the algorithm to create new layers as they are required and to eliminate collapsed layers when they are no longer needed. The algorithm permits a very natural representation of the surface mixed layer, a part of the water column in which surface-driven processes (wind stirring, buoyancy flux, surface current shears) generate vigorous turbulence. In the diapycnal mixing model the mixed layer is distinguished by high turbulent diffusivity, in contrast to the ocean interior, which possesses a diffusivity orders of magnitude smaller. The specification of the thickness of this zone of high diffusivity (the mixed layer depth), and the value of the large diffusivity, may be regarded as extraneous to the diapycnal diffusion algorithm. The only measure that needs to be taken within the diapycnal diffusion algorithm is to allow for the creation and maintenance of vanishingly thin layers at the base of the mixed layer. This is done very simply and transparently. The ‘‘mixed layer’’ in this model representation is no single one of the isopycnal layers: several finite thickness isopycnal layers may participate in the zone of high diffusivity that marks the mixed layer, although at any given time usually only a small number dominate, and the identities of the participants may evolve. For example, the identity of the dominant layer or layers may change as a thin layer from the mixed layer base (or at the surface) inflates to represent decrease (increase) of mixed layer specific volume. We illustrated a mixed layer implementation using Kraus–Turner–Niiler ideas. But we emphasize that the mixed layer implementation is modular. Other formulations of mixed layer models exist: for example, the closure model of Mellor and Yamada (1982), or the K-profile parameterization of Large et al. (1994). Such models explicitly determine profiles of eddy exchange coefficients, and as such could be made to fit naturally into the diapycnal diffusion framework developed herein. The one-dimensional algorithm for diapycnal diffusion could be readily adapted into a threedimensional general circulation model based on isopycnal layers. Indeed it is for this purpose that

318

R.A. de Szoeke, S.R. Springer / Ocean Modelling 5 (2003) 297–323

the rather contrary transformation of the familiar linear diffusion equation in level coordinates to its nonlinear form in isopycnal coordinates was carried out. The natural adaptation of mixed layer parameterizations into the diapycnal diffusion algorithm and so into a general circulation model commends this approach. Implementation of these ideas in an ocean model with a nonlinear equation of state also requires a method of dealing with diapycnal mass fluxes arising from thermobaricity, cabbeling, and densification. Work with such a three-dimensional circulation model will be reported in subsequent papers.

Acknowledgements This research was supported by National Science Foundation grants 9402891, 9907008, the Jet Propulsion Laboratory under the Topex/Poseidon Announcement of Opportunity, and by NASA grant NAGS-4947.

Appendix A. Representation of diapycnal mixing in isopycnal coordinates Consider the one-dimensional continuity equation, turbulent heat diffusion equation, and hydrostatic balance: oq o þ ðqwÞ ¼ 0; ot oz 
oh oh o oh 0 þw ¼ j þ q  h_; ot oz oz oz op ¼
gq: oz

ðA:1Þ ðA:2Þ ðA:3Þ

Here z is the vertical space coordinate, w is vertical velocity, q is in situ density, h is potential temperature, p is pressure, j is a turbulent diffusivity coefficient, q0 is a nondiffusive heat flux (a radiative flux, say; in units of K m s
1 ). The right side of (A.2) is denoted by h_ for convenience. Pressure, potential temperature, and density are related by an equation of state, q ¼ qðp; hÞ:

ðA:4Þ

Potential specific volume is the reciprocal of density when a water parcel is removed adiabatically to the reference pressure (zero, say):
1

a ¼ ½qð0; hÞ : Note that in situ density, q, appearing in (A.1) and (A.3) differs from a
1 . Multiplying (A.2) by da=dh, one obtains
 oa oa o oa q þw ¼ j þ
ahh d  a_ ; ot oz oz oz qg

ðA:5Þ

ðA:6Þ

R.A. de Szoeke, S.R. Springer / Ocean Modelling 5 (2003) 297–323

319

where q ¼ qgah q0 is the buoyancy flux associated with radiation and
2 oh oh d¼j þ q0 oz oz

ðA:7Þ

ðA:8Þ

is the dissipation of temperature. Dissipation, coupled with the convexity of the equation of state, ahh > 0, gives the cabbeling, or densification, effect (McDougall and Garrett, 1992; Davis, 1994). The transformation of partial derivatives when the independent variable is changed from z to a is given by !

 9 > o o zt o > ¼
;> = ot ot a za oa t ðA:9Þ
z
 > o 1 o > > ; ¼ ; oz t za oa t as long as za > 0. Thus (A.6) becomes
 o j q w ¼ zt þ za a_ ¼ zt þ þ
za ahh d: oa za qg

ðA:10Þ

The continuity equation (A.1) becomes, using the first equality of (A.10), o o ðqza Þ þ ðqza a_ Þ ¼ 0 ot oa or, using the hydrostatic relation (A.3), o o pa þ ðpa a_ Þ ¼ 0; ot oa where, from the second equality of (A.10), ! o jðqgÞ2 e  pa a_ ¼
q
pa ahh d: oa pa

ðA:11Þ

ðA:12Þ

ðA:13Þ

When divided by g, e is the diapycnal mass flux, in units of kg m
2 s
1 (positive when upward). If densification is negligible, Eqs. (A.12) and (A.13) give
 o o2 K pa ¼
2
q ; ðA:14Þ ot oa pa where K ¼ jðqgÞ2 is the coefficient of diffusion expressed in units of Pa2 s
1 . If not negligible, densification can be absorbed into the ‘‘radiative flux’’ term, q. An equation identical in form to (A.14) could have been obtained using h as independent coordinate, regardless of densification. However, the use of potential density (specific volume) is standard in physical oceanography. The extension of the present treatment to include effects of salinity (here suppressed) leads naturally to

320

R.A. de Szoeke, S.R. Springer / Ocean Modelling 5 (2003) 297–323

the use of potential density (or some other appropriate variable, such as orthobaric density (de Szoeke et al., 2000)) as vertical coordinate. Eq. (A.14) is the nonlinear diffusion equation quoted in the main text. The derivation given may be readily extended to three space dimensions (de Szoeke and Bennett, 1993). In one dimension, we might have inferred that w ¼ 0 from the Boussinesq approximation of (A.1) and (A.14) would follow immediately from (A.10). The derivation given, only slightly longer, is more correct, does not depend on the Boussinesq approximation, and shows the interplay of mass continuity and heat balance in the argument. The boundary conditions on (A.14) are K
q ¼
qtop ; p ¼ 0; a ¼ atop ; ðA15aÞ pa K
q ¼
qbottom ; p ¼ pb ; a ¼ abottom ; ðA15bÞ pa where qtop , qbottom are the prescribed total buoyancy fluxes at the top and bottom boundaries.

Appendix B The conditions under which a positive result for Dp1nþ1 can be guaranteed from algorithm (20) will be described. The discriminant in (20) may be written as 1 Dpn b2
ac ¼ C 2 þ K0 Dt 1 > C 2 ; 2 pK

ðB:1Þ

where


1 1 1 n Dp2 þ Q0 Dt : C ¼ ð1
aÞDp1n
a 2 2 2

ðB:2Þ

This shows that only a real root of (20) is possible. If the negative sign of the radical in (20) is not taken, then the left side of (19) is certainly negative, contradicting that equality. Hence, the choice of positive sign in (20) is required. To see this, imagine taking the positive sign of the radical in (20). Then the left side of (19) may be written n o9
2a1 ð1
aÞC þ ð1 þ aÞ Kp0KDt þ ð1 þ aÞ½b2
ac 1=2 > > > = 1 6
2a fð1
aÞC þ ð1 þ aÞjCjg ðB:3Þ > > 6
Ca if C > 0 > ; 6 C if C < 0: In either case, the left side of (19) is negative, which establishes the result. To ensure that (20) gives a positive value, Dp1nþ1 P 0;

ðB:4Þ

it is enough to require that both a; c > 0, i.e.,
Dp1n
12Dp2n < Q0 Dt < pQ :

ðB:5Þ

R.A. de Szoeke, S.R. Springer / Ocean Modelling 5 (2003) 297–323

If Q0 > 0, it must be required that pQ Dt < : Q0

321

ðB:6Þ

If Q0 < 0, it is required that Dt <

Dp1n þ 12Dp2n :
Q0

ðB:7Þ

Suppose layer 2 straddles the boundary transition region pQ . Then Eq. (9a) for this layer and layer 3 may be written, using (21), Dp2nþ1 ¼ Dp2


K1 Dt K3 Dt
Q1 Dt

Q3 Dt; Dp1 Dp3

ðB:8Þ

Dp3nþ1 ¼ Dp3


K2 Dt K4 Dt
Q2 Dt

Q4 Dt; Dp2 Dp4

ðB:9Þ

where, from (12) and (14), K1 ¼ 0;

Dp1nþ1 ; K3 ¼ K4 ¼ K0 ; pK
 Dp1nþ1 Q2 ¼ Q0 1
; Q3 ¼ Q4 ¼ 0: pQ

K2 ¼ K0

Q1 ¼ Q0 ;

ðB:10Þ ðB:11Þ

It is readily shown that Dp3nþ1 P 0. From inequality (24) and the equalities given above, 9 > Dp3 P ð2K0 DtÞ1=2 ; > > > >
K4 Dt = 1=2 1 P
ð K DtÞ ; 2 0  ðB:12Þ Dp4 > > > K2 Dt 1=2 > >
ð12K0 DtDp1nþ1 Þ ; > ;  Dp2 so that 0

11=2

C B1 C Dp3nþ1 P B @ 2 K0 DtA 2 P 41


2 41


Dp1nþ1 pQ

Dp1nþ1 pQ

!1=2 3 ! " nþ1 5
Q0 Dt 1
Dp1 pQ

2 !1=2 38 !1=2 < 1 5 K Dt
Q0 Dt41 þ : 2 0

Dp1nþ1 pQ

!1=2 39 = 5 : ;

ðB:13Þ

Since, by assumption, Dp1nþ1 =pQ 6 1, it is easy to see that Dp3nþ1 is positive when Q0 < 0. For Q0 > 0, " ) 1=2
nþ1 1=2 #(
Dp 1 1
2Q0 Dt Dp3nþ1 P 1
ðB:14Þ K0 Dt 2 pQ

322

R.A. de Szoeke, S.R. Springer / Ocean Modelling 5 (2003) 297–323

which is positive if K0 Dt < 2 : 8Q0 It cannot infallibly be guaranteed that Dp2nþ1 P 0. This is for the same basic reason that Dp1nþ1 from (9b) cannot be guaranteed positive. Should Dp2nþ1 become negative, it is recalculated as though it is expected that Dp1nþ1 þ Dp2nþ1 < pQ . This ensures a positive result for Dp2nþ1 .

References Arakawa, A., Hsu, Y.G, 1990. Energy conserving and potential enstrophy dissipating schemes for the shallow water equations. Mon. Wea. Rev. 118, 1960–1969. Bleck, R., 1974. Short-range prediction in isentropic coordinates with filtered and unfiltered numerical models. Mon. Wea. Rev. 102, 813–929. Bleck, R., 2001. An oceanic general circulation model framed in hybrid isopycnic-Cartesian coordinates. Ocean Model. 4, 55–88. Bleck, R., Smith, L.T., 1990. A wind-driven isopycnic model of the North and Equatorial Atlantic Ocean. 1. Model development and supporting experiments. J. Geophys. Res. 95, 3273–3285. Bleck, R., Hanson, H.P., Hu, D., Kraus, E.B., 1989. Mixed layer–thermocline interaction in a three-dimensional isopycnic coordinate model. J. Phys. Oceanogr. 29, 1417–1439. Bleck, R., Rooth, C., Hu, D., Smith, L.T., 1992. Salinity-driven transients in a wind- and thermohaline-forced isopycnic coordinate model of the North Atlantic. J. Phys. Oceanogr. 22, 1486–1505. Davis, R.E., 1994. Diapycnal mixing in the ocean: Equations for large-scale budgets. J. Phys. Oceanogr. 24, 777–800. Davis, R.E., de Szoeke, R., Niiler, P., 1981. Variability in the upper ocean during MILE. Part II: Modeling the mixed layer response. Deep-Sea Res. 28, 1453–1475. de Szoeke, R.A., 2000. Equations of motion using thermodynamic coordinates. J. Phys. Oceanogr. 30, 2814–2829. de Szoeke, R.A., Bennett, A.F., 1993. Microstructure fluxes across density surfaces. J. Phys. Oceanogr. 23, 2254–2264. de Szoeke, R.A., Rhines, P.B., 1976. Asymptotic regimes in mixed-layer deepening. J. Mar. Res. 34, 111–116. de Szoeke, R.A., Springer, S.R., Oxilia, D.M., 2000. Orthobaric density: a thermodynamic variable for ocean circulation studies. J. Phys. Oceanogr. 30, 2830–2852. Dewar, W.K., 2001. Density coordinate mixed layer models. Mon. Wea. Rev. 129, 237–253. Hallberg, R., 2000. Time integration of diapycnal diffusion and Richardson number dependent mixing in isopycnal coordinate ocean models. Mon. Wea. Rev. 128, 1402–1419. Hallberg, R., Rhines, P., 1996. Buoyancy-driven circulation in an ocean basin with isopycnals intersecting the sloping boundary. J. Phys. Oceanogr. 26, 913–940. Hu, D., 1996. The computation of diapycnal diffusive and advective scalar fluxes in multilayer isopycnic-coordinate ocean models. Mon. Wea. Rev. 124, 1834–1851. Kraus, E.B., Turner, J.S., 1967. A one-dimensional model of the seasonal thermocline: II. The general theory and its consequences. Tellus 19, 98–106. Large, W.G., McWilliams, J.C., Doney, S.C., 1994. Oceanic vertical mixing: a review and a model with nonlocal boundary layer parameterizations. Rev. Geophys. 32, 363–403. McDougall, T.J., 1987. Thermobaricity, Cabbeling and water mass conversion. J. Geophys. Res. 92, 5448–5464. McDougall, T., Dewar, W., 1998. Vertical mixing and cabbeling in layered models. J. Phys. Oceanogr. 28, 1458–1480. McDougall, T., Garrett, C.J.R., 1992. Scalar conservation equations in a turbulent ocean. Deep-Sea Res. 39, 1953– 1966. Mellor, G.L., Yamada, T., 1982. Development of a turbulence closure model for geophysical fluid problems. Rev. Geophys. Space Phys. 20, 851–875. Niiler, P.P., 1975. Deepening of the wind-mixed layer. J. Mar. Res. 33, 405–422.

R.A. de Szoeke, S.R. Springer / Ocean Modelling 5 (2003) 297–323

323

Niiler, P.P., Kraus, E.B., 1977. One-dimensional models of the upper ocean. In: Kraus, E.B. (Ed.), Modeling Prediction of the Upper Layers of the Ocean. Pergamon Press, pp. 143–172. Oberhuber, J.M., 1993. Simulation of the Atlantic circulation with a coupled sea ice-mixed layer-isopycnal general circulation model. Part I: Model description. J. Phys. Oceanogr. 23, 808–829. Pollard, R.T., Rhines, P.B., Thompson, R.O.R.Y., 1973. The deepening of the wind-mixed layer. Geophys. Fluid Dyn. 3, 381–404. Sun, S., Bleck, R., Rooth, C., Dukowicz, J., Chassignet, E., Killworth, P., 1999. Inclusion of thermobaricity in isopycnic-coordinate ocean models. J. Phys. Oceanogr. 29, 2719–2729. Turner, J.S., 1981. Small-scale mixing processes. In: Wunsch, C., Warren, B.A. (Eds.), Evolution of Physical Oceanography. MIT Press, Cambridge, MA, pp. 236–262.