G. \underline{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 γ_{x}:

\overline{wx}(d) = −K_{x}(∂_{z}X − γ_{x}) .  (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} = d_{n+0.5} − d_{n−0.5}, and the distances between grid levels are Δ_{n+0.5} = d_{n+1} − d_{n}. To define the latter for n = 0 and n = M, where M is the total number of layers, d_{0} 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 d_{k−1} ≤ h < d_{k}. A useful variable that varies from 0 to 1 over this grid interval is

δ = \frac{(h − d_{k−1})}{Δ_{k−.5}}.  (G2)

G-1

Property values are defined at the grid levels, d_{n}, for 1 ≤ n ≤ M + 1, where the M + 1 values are prescribed bottom boundary conditions. In computing near surface reference values, X_{r}, as the average value between the surface and ϵh, the following continuous profile is assumed:

\begin{aligned}X(d) & = X_{1} ,  d < d_{1} \\ & = X_{n−1} + \frac{(X_{n} − X_{n−1})}{(d_{n} − d_{n−1})}(d − d_{n−1}) d_{1} ≤ d_{n−1} ≤ d < d_{n} ≤ D. \end{aligned}  (G3)

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

\begin{aligned}[∂_{z}X]_{n+.5} & = \frac{X_{n} − X_{n+1}}{Δ_{n+.5}}  , 1 ≤ n ≤ M \\ [∂_{z}X]_{n} & = \frac{X_{n−1} − X_{n+1}}{Δ_{n−.5} + Δ_{n+.5}} , 1 < n ≤ M. \\ & = \frac{X_{1} − X_{2}}{2Δ_{1.5}}  , n = 1. \end{aligned}  (G4)

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) d_{k−1} < h < d_{k−.5} (left panel) and (b) d_{k−.5} < h < d_{k} (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:

\begin{aligned}ν_{x}(h) & = ν_{x}(d_{n+.5}) + ∂_{z}ν_{x}(h)(d_{n+.5} − h) , d_{n−.5} < h ≤ d_{n+.5} \\ ∂_{z}ν_{x}(h) & = (1 − R) [\frac{ν_{x}(d_{n−.5}) − ν_{x}(d_{n+.5})}{Δ_{n}}] + R [\frac{ν_{x}(d_{n+.5}) − ν_{x}(d_{n+1.5})}{Δ_{n+1}}] , \end{aligned} (G5)

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 − d_{n−.5})Δ_{n}^{−1}.

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 d_{k−1.5} and d_{k+.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}(d_{k−.5}) (triangle in Fig. G1) and an enhanced boundary layer diffusivity, K_{x}^{∗}, such that for case (b) situations

\begin{aligned}Λ_{x} & = (1 − δ) ν_{x}(d_{k−.5}) + δ K_{x}^{∗} , \\ K_{x}^{∗} & = (1 − δ)^{2} K(d_{k−1}) + δ^{2} K(d_{k−.5}) , \end{aligned}  (G6)

where δ, as defined in (G2), varies from 0 to 1. For case (a), ν_{x}(d_{k−.5}) replaces K(d_{k−.5}) in (G6), because the latter is not defined for d > h. The important feature is that the dependency on K(d_{k−1}) (square in Fig. G1) leads to an enhanced diffusivity at the k − .5 interface as soon as h becomes greater than d_{k−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 d_{k−.5} = 15m, as δ 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×10^{−4}m^{2}/s at 25, 20, 15, and 10m depth, respectively. The boundary layer diffusivities at the k−.5 interface in the range 0 ≤ δ ≤ 0.5 are constant at the interior value. The importance of using Λ_{x} is to enhance these diffusivities in this range of δ, 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 δ approaches 1, the form of (G6) makes Λ_{x} less than K_{x}(d_{k−.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:

A_{x}^{i} X_{t+1}^{i+1} = X_{t} + H_{x}^{i} ,  (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 X_{t+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 X_{t} values. In some numerical

G-3 schemes the prognostic variables are updated by advection prior to vertical diffusion, so that these updated X_{t+1}^{0} values could also be used.

Equation (G7) has been iterated until the new boundary layer depth from the i^{th} iteration, h^{i+1}, differs from the previous h^{i} by less than a specified tolerance, η_{h} :

\frac{|h^{i} − h^{i+1}|}{Δ_{k}} < η_{h} ,  (G8)

where again d_{k−1} ≤ h^{i+1} ≤ d_{k}, defines the vertical index, k, and Δ_{k} is the local vertical resolution. Iteration allows the X_{t+1}^{i} values used to determine A_{x}^{i} and H_{x}^{i} to be close to the final X_{t+1}^{i+1} values. Either X_{t} or X_{t+1}^{0} from above could be used for the first iteration, i = 1. Alternatively, linear extrapolations of X_{t−1} and X_{t} can be used to give the first iteration X_{t+1}^{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 A_{n,m}, where n is the row index and m the column. The non–zero elements are:

\begin{aligned}A_{1,1} & = (1 + Ω_{1}^{+}) \\ A_{n,n−1} & = −Ω_{n}^{−}  2 ≤ n ≤ M \\ A_{n,n} & = (1 + Ω_{n}^{−} + Ω_{n}^{+})  2 ≤ n ≤ M \\ A_{n,n+1} & = −Ω_{n}^{+}  1 ≤ n ≤ M − 1 , \end{aligned}  (G9)

where

Ω_{n}^{−} = \frac{Δt}{Δ_{n}} \frac{K_{x}(d_{n−.5})}{Δ_{n−.5}} and Ω_{n}^{+} = \frac{Δt}{Δ_{n}} \frac{K_{x}(d_{n+.5})}{Δ_{n+.5}}.  (G10)

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, Q_{n}(d) for temperature and −F_{n}(d) for salinity, with J_{o} the surface value. For a surface turbulent flux of \overline{ws}_{o}, the elements of H_{s} are :

H_{1} = \frac{Δt}{Δ_{1}} ((K_{s} γ_{s})_{1.5} − \overline{ws}_{o} + J_{1.5} − J_{o})  (G11)
H_{n} = \frac{Δt}{Δ_{n}} ((K_{s} γ_{s})_{n+.5} − (K_{s} γ_{s})_{n−.5} + J_{n+.5} − J_{n−.5})  2 ≤ n ≤ M − 1
H_{M} = \frac{Δt}{Δ_{M}} ((K_{s} γ_{s})_{M+.5} − (K_{s} γ_{s})_{M−.5} + J_{M+.5} − J_{M−.5} + S_{M+1} \frac{K_{s}(d_{M+.5})}{Δ_{M+.5}}) ,

where S_{M+1} is the prescribed bottom boundary value of the scalar and subscripts n + .5 and n − .5 denote evaluation at d_{n−.5} and d_{n+.5}, respectively.

G-4

For velocity components, H includes only boundary conditions and Coriolis terms. For the zonal velocity component the elements of H_{u} are:

H_{1} = Δt f V_{1}^{t+.5} − \frac{Δt}{Δ_{1}} \overline{wu}_{o}
H_{n} = Δt f V_{n}^{t+.5}  2 ≤ n ≤ M − 1
H_{M} = Δt f V_{M}^{t+.5} + \frac{Δt}{Δ_{M}} \frac{U_{M+1} K_{m}(d_{M+.5})}{Δ_{M+.5}} .

(G12)

For the meridional velocity component the elements of H_{v} are:

H_{1} = −Δt f U_{1}^{t+.5} − \frac{Δt}{Δ_{1}} \overline{wv}_{o}
H_{n} = −Δt f U_{n}^{t+.5}  2 ≤ n ≤ M − 1
H_{M} = −Δt f U_{M}^{t+.5} + \frac{Δt}{Δ_{M}} \frac{V_{M+1} K_{m}(d_{M+.5})}{Δ_{M+.5}} .

(G13)

In (G12) and (G13) the fixed bottom boundary conditions are U_{M+1} and V_{M+1}. The Coriolis terms are estimates of their average value over the time step:

\begin{aligned}U_{n}^{t+.5} & = 0.5 (U_{n}^{t} + U_{n}^{i}) \\ V_{n}^{t+.5} & = 0.5 (V_{n}^{t} + V_{n}^{i}) , \end{aligned}  (G14)

where the layer n velocity components at time t are U_{n}^{t} and V_{n}^{t}.

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 d_{k−1} and d_{k−1/2}, and (b) h between d_{k−1/2} and d_{k}. Shown are the discrete interior diffusivities (triangles) and their interpolation (dashed curve) from (G6); the continuous boundary layer diffusivity profile K_{x}(d) (solid trace), including its value at d_{k−1} (square); and finally the modified diffusivity Λ_{x} at d_{k−1/2}, which is used by the model (cross).

Figure G2. Comparison of Λ_{x} used by the model and K_{x}(d_{k−0.5}) as δ varies from 0 to 1. Shown are a case of no interior mixing (Λ_{x}, dashed trace; K_{x}, dot-dashed trace) and a case with substantial interior mixing (Λ_{x}, solid trace; K_{x}, dotted trace).

G-6