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