G. KPP__Boundary__Layer__Mixing__Scheme_____________________________
The nonlocal K-profile parameterization of Troen and Mahrt (1986) has been adapted
for use as an oceanic boundary layer model. Important characteristics of both applications
are that they are consistent with similarity theory in the surface layer, the boundary layer
is capable of penetrating the interior stratification, and turbulent transport vanishes at the
surface. The OBL model also has additional desirable features. For example, turbulent
shear contributes to the diagnosed boundary layer depth, so as to make the entrainment
of buoyancy at the base of the OBL independent of the interior stratification. Interior
mixing at the base of the boundary layer (d = h) influences the turbulence throughout
the boundary layer. Also, in the convective limit the turbulent velocity scales for both
momentum and scalars become directly proportional to w* . The problem of determining
the vertical turbulent fluxes of momentum and both active and passive scalars throughout
the OBL is closed by adding a nonlocal transport term flx :
_____wx(d) = -K
x (@z X - flx ) : (G1)
In practice the external forcing is first prescribed, then the boundary layer depth, h, is
determined, and finally profiles of the diffusivity and nonlocal transport are computed. A
complete description of this model can be found in Large et al. (1994), and what follows
is based on Appendix D of that paper. The KPP scheme is activated in the model by
specifying the ifdef option kmix. The ifdef option implicitvmix must also be specified,
because the vertical mixing is done implicitly. Specifying this ifdef option also eliminates
the original MOM convective adjustment scheme. The ifdef option constvmix must not be
specified with kmix, so that the original MOM subroutine cnvmix is also bypassed when
the KPP scheme is used. There are an additional three sub-options that can be specified
with kmix. If ifdef option kmixcheckekmo is activated, then the code does an additional
check against the Ekman and Monin-Obukhov length scales. If ifdef option kmixnori is
specified, then the vertical viscosity and diffusivity below the boundary layer are just set
to the background values, fkpm and fkph and are not dependent on the local Richardson
number. If ifdef option kmixdd is specified, then a double diffusion contribution is added
to the vertical diffusivity, see Large et al. (1994).
Vertical discretization
Layer grid points are denoted by whole number indices, whereas layer interfaces are
denoted by half number indices. The layer thicknesses are n = dn+0:5 - dn-0:5 , and the
distances between grid levels are n+0:5 = dn+1 - dn . To define the latter for n = 0 and
n = M , where M is the total number of layers, d0 is taken to be zero and a fictitious layer,
M + 1, of zero thickness is added at the bottom. An index k can always be found such
that dk-1 h < dk . A useful variable that varies from 0 to 1 over this grid interval is
ffi = (h_-_dk-1__)____: (G2)
k-:5
G-1
Property values are defined at the grid levels, dn , for 1 n M + 1, where the M + 1
values are prescribed bottom boundary conditions. In computing near surface reference
values, Xr , as the average value between the surface and fflh, the following continuous
profile is assumed:
X (d) = X1 ; d < d1
(G3)
= Xn-1 + (Xn__-_Xn-1___)_____(d(d - dn-1 ) d1 dn-1 d < dn D:
n - dn-1 )
Property gradients are required both at interfaces and grid levels, and in the model where
partial derivatives with respect to d equal -@z these are computed, respectively, as
[@z X ]n+:5 = Xn__-_Xn+1_______ ; 1 n M
n+:5
[@z X ]n = _Xn-1___-_Xn+1________ ; 1 < n M: (G4)
n-:5 + n+:5
= X1__-_X2_____2 ; n = 1:
1:5
Interior_boundary layer matching
The convective and wind deepening cases are handled well, but not ideally by a low
vertical resolution KPP model. Several techniques were explored in an attempt to dampen
oscillations and reduce the bias in h. One approach is to use nonlinear interpolation.
Several schemes were tested and they could be made to work very well with idealized
profiles. However, in general the property profiles, and especially the bulk Richardson
number profile, are highly variable and no one scheme could be found to work well over a
wide variety of naturally occurring conditions.
The following practical scheme is used to remove the bias and dampen oscillations. It
involves enhancing the diffusivity at the k - :5 interface. Figure G1 shows the diffusivities
produced by the model's parameterizations in the vicinity of layer k for two situations:
(a) dk-1 < h < dk-:5 (left panel) and (b) dk-:5 < h < dk (right panel). The first
step in the matching process is to determine the interior diffusivities (triangles). These
are interpolated (dotted line) to give x (h) and the gradient, @z x (h), at h. The only
important requirement here is that there not be discontinuities in either quantity as h
progresses through the grid. This is true of the simple scheme used in the model:
x (h) = x (dn+:5 ) + @z x (h)(dn+:5 - h) ; dn-:5 < h dn+:5
h (d ) - (d ) i h (d ) - (d ) i (G5)
@z x (h) = (1 - R) x___n-:5______x___n+:5__________ + R x___n+:5______x___n+1:5___________ ;
n n+1
with n = k - 1 in situation (a) and n = k in situation (b). The interpolated gradient is
just a weighted average of two discrete gradients, with the weight R = (h - dn-:5 )-1n .
G-2
The continuous boundary layer diffusivity profile is the solid line in Fig. G1. The dotted
line in Fig. G1 follows the x (h) that would be computed from (G5) for different values of
h between dk-1:5 and dk+:5 .
The next step is to compute a modified diffusivity, x (cross in Fig. G1), which is to
be applied at the k - :5 interface. It is a weighted average of x (dk-:5 ) (triangle in Fig. G1)
and an enhanced boundary layer diffusivity, K*x, such that for case (b) situations
x = (1 - ffi) x (dk-:5 ) + ffi K*x ;
(G6)
K*x = (1 - ffi)2 K(dk-1 ) + ffi2 K(dk-:5 ) ;
where ffi, as defined in (G2), varies from 0 to 1. For case (a), x (dk-:5 ) replaces K(dk-:5 )
in (G6), because the latter is not defined for d > h. The important feature is that
the dependency on K(dk-1 ) (square in Fig. G1) leads to an enhanced diffusivity at the
k - :5 interface as soon as h becomes greater than dk-1 . The increased deepening that
results greatly reduces the boundary layer depth bias of low, relative to high, resolution
simulations.
Figure G2 compares the diffusivity, x , used by the model to the boundary layer
diffusivity at dk-:5 = 15m, as ffi varies from 0 to 1. The latter are kept small by using a
small u* = 0:006ms-1 and by using the grid of Fig. G1 with k = 4, so that h does not get
very large. Two extreme cases are shown: no interior mixing (dashed versus dot-dashed
lines) and substantial interior mixing (solid versus dotted lines). In the latter case, the
interior diffusivities are 1; 2; 4, and 8 x 10-4 m2 =s at 25; 20; 15, and 10m depth, respectively.
The boundary layer diffusivities at the k - :5 interface in the range 0 ffi 0:5 are constant
at the interior value. The importance of using x is to enhance these diffusivities in this
range of ffi, when the interior diffusivity is small. When there is substantial interior mixing
as in the cases shown in Figs. G1 and G2, the enhancement is not necessary and is much
reduced. As ffi approaches 1, the form of (G6) makes x less than Kx (dk-:5 ) in an attempt
to reduce biases.
Semi-implicit time integration
The prognostic equations of the model are solved by a semi-implicit integration scheme
whose general matrix form is:
Aix Xi+1t+1= Xt + Hix ; (G7)
where A is an M by M tridiagonal matrix, and X and H are vectors of length M .
Integration over a time step, t, is accomplished by inversion of A, which allows the
properties at the new time t + 1 to be computed from past values at t. The integration
can only be semi-implicit, because A and H depend on quantities like the diffusivity that
depend on h, which in turn depend on the profiles Xt+1 that are themselves computed
from A and H. The superscripts in (G7) denote various choices of when and how A and
H are calculated. The simplest method, denoted by i = 0, would be to compute all the
required quantities, including the forcing, at time t using Xt values. In some numerical
G-3
schemes the prognostic variables are updated by advection prior to vertical diffusion, so
that these updated X0t+1 values could also be used.
Equation (G7) has been iterated until the new boundary layer depth from the ith
iteration, hi+1 , differs from the previous hi by less than a specified tolerance, jh :
|hi_-_hi+1__|___ < j ; (G8)
k h
where again dk-1 hi+1 dk , defines the vertical index, k, and k is the local vertical
resolution. Iteration allows the Xit+1 values used to determine Aix and Hix to be close to
the final Xi+1t+1values. Either Xt or X0t+1 from above could be used for the first iteration,
i = 1. Alternatively, linear extrapolations of Xt-1 and Xt can be used to give the first
iteration X1t+1 values. Since it is undesirable to use the extrapolated values in the final
iteration, at least two iterations are always performed. Over an annual cycle at OWS
Papa less than 1% of all time integrations required more than two iterations. At most 12
iterations were needed, but 0:3 percent of all iterations (when the stratification was weak)
failed to converge and the results of the twentieth iteration were used.
The elements of A have the same form for all properties, so the subscript x can be
dropped for convenience. The elements of A then become An;m , where n is the row index
and m the column. The non-zero elements are:
A1;1 = (1 + +1)
An;n-1 = --n 2 n M
(G9)
An;n = (1 + -n + +n) 2 n M
An;n+1 = -+n 1 n M - 1 ;
where
-n = t____ Kx_(dn-:5__)____ and +n = t____ Kx_(dn+:5__)___ : (G10)
n n-:5 n n+:5
The vector H is very different for scalars than for velocity. In general, the scalar H includes
the boundary conditions, countergradient terms, and non-turbulent forcing. Let J be the
non-turbulent forcing at the interfaces, Qn (d) for temperature and -Fn (d) for salinity,
with Jo the surface value. For a surface turbulent flux of _____wso, the elements of Hs are :
H1 = t____ (Ks fls )1:5 - _____wso+ J1:5 - Jo (G11)
1
Hn = t____ (Ks fls )n+:5 - (Ks fls )n-:5 + Jn+:5 - Jn-:5 2 n M - 1
n
HM = _t____ (Ks fls )M +:5 - (Ks fls )M -:5 + JM +:5 - JM -:5 + SM +1 Ks_(dM_+:5__)___ ;
M M +:5
where SM +1 is the prescribed bottom boundary value of the scalar and subscripts n + :5
and n - :5 denote evaluation at dn-:5 and dn+:5 , respectively.
G-4
For velocity components, H includes only boundary conditions and Coriolis terms. For
the zonal velocity component the elements of Hu are:
H1 = t f V1t+:5 - _t___ _____wuo
1
Hn = t f Vnt+:5 2 n M - 1 (G12)
HM = t f VMt+:5 + _t____ UM_+1____Km__(dM_+:5__)____ :
M M +:5
For the meridional velocity component the elements of Hv are:
H1 = -t f U1t+:5 - _t___ _____wvo
1
Hn = -t f Unt+:5 2 n M - 1 (G13)
HM = -t f UMt+:5 + _t____ VM_+1____Km__(dM_+:5__)____ :
M M +:5
In (G12) and (G13) the fixed bottom boundary conditions are UM +1 and VM +1 . The
Coriolis terms are estimates of their average value over the time step:
Unt+:5 = 0:5 (Unt + Uni)
(G14)
Vnt+:5 = 0:5 (Vnt + Vni) ;
where the layer n velocity components at time t are Unt and Vnt.
Large, W.G., J.C. McWilliams, and S.C. Doney, 1994: Oceanic vertical mixing: A review
and a model with a nonlocal boundary layer parameterization. Rev. of Geophys., 32,
363-403.
Troen, I.B., and L. Mahrt, 1986: A simple model of the atmospheric boundary layer;
sensitivity to surface evaporation. Boundary-Layer Meteor., 37, 129-148.
G-5
Figure G1. Schematic of the diffusivities required to match the interior and boundary
layer mixing for two cases: (a) h between dk-1 and dk-1=2 , and (b) h between dk-1=2 and
dk . Shown are the discrete interior diffusivities (triangles) and their interpolation (dashed
curve) from (G6); the continuous boundary layer diffusivity profile Kx (d) (solid trace),
including its value at dk-1 (square); and finally the modified diffusivity x at dk-1=2 ,
which is used by the model (cross).
Figure G2. Comparison of x used by the model and Kx (dk-0:5 ) as ffi varies from 0 to 1.
Shown are a case of no interior mixing (x , dashed trace; Kx , dot-dashed trace) and a
case with substantial interior mixing (x , solid trace; Kx , dotted trace).
G-6