Energetically-constrained Planetary Boundary Layer

We here describe a scheme for modeling the ocean surface boundary layer (OSBL) suitable for use in global climate models. It builds on the ideas in Bulk Surface Mixed Layer, bringing in some of the ideas from Shear-driven mixing in Jackson, to make an energetically consistent boundary layer suitable for use with a generalized vertical coordinate. Unlike in Bulk Surface Mixed Layer, variables are allowed to have vertical structure within the boundary layer. The downward turbulent flux of buoyant water by OSBL turbulence converts mechanical energy into potential energy as it mixes with less buoyant water at the base of the OSBL. As described in [53], we focus on OSBL parameterizations that constrain this integrated potential energy conversion due to turbulent mixing.

The leading-order mean OSBL equation for arbitrary scalar \(\phi\) is:

\[\frac{\partial \overline{\phi}}{\partial t} = - \frac{\partial}{\partial z} \overline{w^\prime \phi^\prime} + \nu_\phi \frac{\partial^2 \overline{\phi}}{\partial z^2}\]

where the symbols are as follows:

Symbol	Meaning
\(u_i\)	horizontal components of the velocity
\(\phi\)	arbitrary scalar (tracer) quantity
\(w\)	vertical component of the velocity
\(\overline{w}\)	ensemble average \(w\)
\(w^\prime\)	fluctuations from \(\overline{w}\)
\(k\)	turbulent kinetic energy (TKE)
\(K_M\)	turbulent mixing coefficient for momentum
\(K_\phi\)	turbulent mixing coefficient for \(\phi\)
\(\sigma_k\)	turbulent Schmidt number
\(b\)	buoyancy
\(\epsilon\)	buoyancy turbulent dissipation rate

This equation describes the evolution of mean quantity \(\overline{\phi}\) due to vertical processes, including the often negligible molecular mixing. We would like to parameterize the vertical mixing since we won’t be resolving all the relevant time and space scales.

We use the Boussinesq hypothesis for turbulence closure. This approximates the Reynolds stress terms using an eddy viscosity (eddy diffusivity for turbulent scalar fluxes):

\[\overline{u_i^\prime w^\prime} = - K_M \frac{\partial \overline{u_i}}{\partial z} ,\]

Similarly, the eddy diffusivity is used to parameterize turbulent scalar fluxes as:

\[\overline{\phi^\prime w^\prime} = - K_\phi \frac{\partial \overline{\phi}}{\partial z} ,\]

The parameters needed to close the system of equations are then reduced to the turbulent mixing coefficients, \(K_\phi\) and \(K_M\).

We start with an equation for the turbulent kinetic energy (TKE):

\[\frac{\partial k}{\partial t} = \frac{\partial}{\partial z} \left( \frac{K_M}{\sigma_k} \frac{\partial k}{\partial z} \right) - \overline{u_i^\prime w^\prime} \frac{\partial \overline{u_i}} {\partial z} + \overline{w^\prime b^\prime} - \epsilon\]

Terms in this equation represent TKE storage (LHS), TKE flux convergence, shear production, buoyancy production, and dissipation.

Well-mixed Boundary Layers (WMBL)

Assuming steady state and other parameterizations, integrating vertically over the surface boundary layer, [53] obtains the form:

\[\frac{1}{2} H_{bl} w_e \Delta b = m_\ast u_\ast^3 - n_\ast \frac{H_{bl}}{2} B(H_{bl}) ,\]

with the following variables:

Symbol	Meaning
\(H_{bl}\)	boundary layer thickness
\(w_e\)	entrainment velocity
\(\Delta b\)	change in buoyancy at base of mixed layer
\(m_\ast\)	sum of mechanical coefficients
\(u_\ast\)	friction velocity ( \(u_\ast = (|\tau| / \rho_0)^{1/2}\))
\(\tau\)	wind stress
\(n_\ast\)	convective proportionality coefficient
	1 for stabilizing surface buoyancy flux, less otherwise
\(B(H_{bl})\)	surface buoyancy flux

Energetics-based Planetary Boundary Layer

Once again, the goal is to formulate a surface mixing scheme to find the turbulent eddy diffusivity (and viscosity) in a way that is suitable for use in a global climate model, using long timesteps and large grid spacing. After evaluating a well-mixed boundary layer (WMBL), the shear mixing of [36] (JHL, Shear-driven mixing in Jackson), as well as a more complete boundary layer scheme, it was decided to combine a number of these ideas into a new scheme:

\[K(z) = F_x(K_{ePBL}(z), K_{JHL}(z), K_n(z))\]

where \(F_x\) is some unknown function of a new \(K_{ePBL}\), \(K_{JHL}\), the diffusivity due to shear as determined by [36], and \(K_n\), the diffusivity from other ideas. We start by specifying the form of \(K_{ePBL}\) as being:

\[K_{ePBL}(z) = C_K w_t l ,\]

where \(w_t\) is a turbulent velocity scale, \(C_K\) is a coefficient, and \(l\) is a length scale.

Turbulent length scale

We propose a form for the length scale as follows:

\[l = (z_0 + |z|) \times \max \left[ \frac{l_b}{H_{bl}} , \left( \frac{H_{bl} - |z|}{H_{bl}} \right)^\gamma \, \right] ,\]

where we have the following variables:

Symbol	Meaning
\(H_{bl}\)	boundary layer thickness
\(z_0\)	roughness length
\(\gamma\)	coefficient, 2 is as in KPP, [40]
\(l_b\)	bottom length scale

Turbulent velocity scale

We do not predict the TKE prognostically and therefore approximate the vertical TKE profile to estimate \(w_t\). An estimate for the mechanical contribution to the velocity scale follows the standard two-equation approach. In one and two-equation second-order \(K\) parameterizations the boundary condition for the TKE is typically employed as a flux boundary condition.

\[K \left. \frac{\partial k}{\partial z} \right|_{z=0} = c_\mu^0 u_\ast^3 .\]

The profile of \(k\) decays in the vertical from \(k \propto (c_\mu^0)^{2/3} u_\ast^2\) toward the base of the OSBL. Here we assume a similar relationship to estimate the mechanical contribution to the TKE profile. The value of \(w_t\) due to mechanical sources, \(v_\ast\), is estimate as \(v_\ast (z=0) \propto (c_\mu^0)^{1/3} u_\ast\) at the surface. Since we only parameterize OSBL turbulent mixing due to surface forcing, the value of the velocity scale is assumed to decay moving away from the surface. For simplicity we employ a linear decay in depth:

\[v_\ast (z) = (c_\mu^0)^{1/3} u_\ast \left( 1 - a \cdot \min \left[ 1, \frac{|z|}{H_{bl}} \right] \right) ,\]

where \(1 > a > 0\) has the effect of making \(v_\ast(z=H_{bl}) > 0\). Making the constant coefficient \(a\) close to one has the effect of reducing the mixing rate near the base of the boundary layer, thus producing a more diffuse entrainment region. Making \(a\) close to zero has the effect of increasing the mixing at the base of the boundary layer, producing a more ‘step-like’ entrainment region.

An estimate for the buoyancy contribution is found utilizing the convective velocity scale:

\[w_\ast (z) = C_{w_\ast} \left( \int_z^0 \overline{w^\prime b^\prime} dz \right)^{1/3} ,\]

where \(C_{w_\ast}\) is a non-dimensional empirical coefficient. Convection in one and two-equation closure causes a TKE profile that peaks below the surface. The quantity \(\overline{w^\prime b^\prime}\) is solved for in ePBL as \(KN^2\).

These choices for the convective and mechanical components of the velocity scale in the OSBL are then added together to get an estimate for the total turbulent velocity scale:

\[w_t (z) = w_\ast (z) + v_\ast (z) .\]

The value of \(a\) is arbitrarily chosen to be 0.95 here.

Summarizing the ePBL implementation

The ePBL mixing coefficient is found by multiplying a velocity scale (Turbulent velocity scale) by a length scale (Turbulent length scale). The precise value of the coefficient \(C_K\) used does not significantly alter the prescribed potential energy change constraint. A reasonable value is \(C_K \approx 0.55\) to be consistent with other approaches (e.g. [64]).

The boundary layer thickness ( \(H_{bl}\)) within ePBL is based on the depth where the energy requirement for turbulent mixing of density exceeds the available energy (Well-mixed Boundary Layers (WMBL)). \(H_{bl}\) is determined by the energetic constraint imposed using the value of \(m_\ast\) and \(n_\ast\). An iterative solver is required because \(m_\ast\) and the mixing length are dependent on \(H_{bl}\).

We use a constant value for convectively driven TKE of \(n_\ast = 0.066\). The parameterizations for \(m_\ast\) are formulated specifically for the regimes where \(K_{JHL}\) is sensitive to model numerics \((|f| \Delta t \approx 1)\) ([53]).

Combining ePBL and JHL mixing coefficients

We now address the combination of the ePBL mixing coefficient and the JHL mixing coefficient. The function \(F_x\) above cannot be the linear sum of \(K_{ePBL}\) and \(K_{JHL}\). One reason this sum is not valid is because the JHL mixing coefficient is determined by resolved current shear, including that driven by the surface wind. The wind-driven current is also included in the ePBL mixing coefficient formulation. An alternative approach is therefore needed to avoid double counting.

\(K_{ePBL}\) is not used at the equator as scalings are only investigated when \(|f| > 0\). The solution we employ is to use the maximum mixing coefficient of the two contributions,

\[K (z) = \max (K_{ePBL} (z), K_{JHL} (z)),\]

where \(m_\ast\) (and hence \(K_{ePBL}\)) is constrained to be small as \(|f| \rightarrow 0\). This form uses the JHL mixing coefficient when the ePBL coefficient is small.

This approach is reasonable when the wind-driven mixing dominates, since both JHL and ePBL give a similar solution when deployed optimally. One weakness of this approach is the tropical region, where the shear-driven ePBL \(m_\ast\) coefficient is not formulated. The JHL parameterization is skillful to simulate this mixing, but does not include the contribution of convection. The convective portion of \(K_{ePBL}\) should be combined with \(K_{JHL}\) in the equatorial region when shear and convection occur together. Future research is warranted.

Finally, one should note that the mixing coefficient here ( \(K\)) is used for both diffusivity and viscosity, implying a turbulent Prandtl number of 1.0.

Langmuir circulation

While only briefly alluded to in [53], the MOM6 code implementing ePBL does support the option to add a Langmuir parameterization. There are in fact two options, both adjusting \(m_\ast\).