KPP Boundary Layer Mixing

Table of Contents

1. Overview

2. Requirements

3. Algorithmic Formulation

4. Design

5. Verification and Testing

1 Overview

This document describes the Omega implementation of the K Profile Parameterization (KPP) ocean boundary layer mixing. KPP computes a boundary-layer depth, vertical viscosity, vertical diffusivity, and a non-local tracer flux shape implemented outside the implicit vertical mixing routine. The implementation follows that in MPAS-Ocean and uses direct ports of the functions defined in the CVMix version of KPP.

The implementation is in KPPMix and is integrated with the Omega tendency and RK2, RK4, and Forward-backward stepping routines. Relative to broad vertical mixing documentation, this page focuses specifically on KPP theory, algorithmic choices, and testing.

2 Requirements

2.1 Requirement: Boundary-layer depth from bulk Richardson criterion

Following Large et al (1994), the OBL depth must be diagnosed from a bulk Richardson criterion so that mixing depth responds to evolving stratification, shear, and surface forcing. It also must include a unresolved turbulent shear contribution.

2.2 Requirement: Coefficients must be computable in parallel over columns

The KPP implementation must operate over many columns in parallel using Omega array/kernels, rather than serial single-column calls.

2.3 Requirement: Compatible with additive vertical-mixing framework

KPP viscosity/diffusivity fields must be compatible with existing Omega vertical mixing infrastructure so that other chosen vertical mixing sources can be merged with KPP.

2.4 Desired: Non-local flux and profile matching controls

KPP will support a non-local tracer flux from LMD94 and include configurable viscosity/diffusivity matching at the base of the boundary layer.

3 Algorithmic Formulation

The implementation follows a two-stage KPP structure.

3.1 Stage 1: OBL depth search

For each water column, OBL depth \(h\) is diagnosed by searching downward until bulk Richardson number reaches a critical value:

\[ Ri_b(z) = \frac{\Delta b(z)\, z}{|\Delta \mathbf{U}(z)|^2 + V_t^2(z)} \]

with threshold

\[ Ri_b(h) = Ri_{crit}. \]

Here, \(\Delta b\) is buoyancy jump relative to a surface layer average, \(|\Delta \mathbf{U}|^2\) is shear contribution again computed relative to the surface layer average, and \(V_t^2\) is unresolved shear.

When the bulk Richardson number falls between model layers, quadratic interpolation is utilized to find the depth. In addition the boundary layer depth is constrained to fall between a configurable minimum OBL under sea ice and a maximum set by the water column depth.

The boundary layer depth search is a growing inner loop. The outer loop iterates over all model layers. The current model layer is set as a boundary layer depth candidate and \(Ri_b\) is calculated for all model layers shallower than the current depth. If any layer in the inner loop has an \(Ri_b*StopOBL\) that exceeds \(Ri_{crit}\) the loop terminates.

3.2 Stage 2: KPP coefficients and non-local flux

Given a diagnosed \(h\), KPP computes viscosity and diffusivity coefficients at the top of every Omega cell except the surface and the bottom using shape functions in normalized depth \(\sigma = -d/h\), where \(d\) is the depth relative to the sea surface height, not the physical depth:

\[ K_m(\sigma) = h\, w_m(\sigma)\, M_1(\sigma), \]

\[ K_s(\sigma) = h\, w_s(\sigma)\, S_1(\sigma), \]

where \(w_m\) and \(w_s\) are turbulent velocity scales from Monin-Obukhov style stability functions, see Appendix B of Large et al, 1994. \(M_1\) and \(S_1\) are shape functions. The generic form of the shape function is given by

\[ X(\sigma) = c_1 \sigma^3 + c_2 \sigma^2 + c_3 \sigma + c_4 \]

The coefficients are determined by various conditions, e.g., zero viscosity and diffusivity at the surface, assuming a linear reduction of the turbulent flux with distance from the surface in the surface layer. As in MPAS-Ocean, we include two options to determine the final coefficients. The original version of KPP matches predicted viscosities and diffusivities to those predicted by other schemes (e.g., shear instability driven mixing) and a second option where viscosities and diffusivities are instead additive. In the latter case, the shape function greatly simplifies to $X(\sigma) = \sigma(1-\sigma)^2. When diffusivities and viscosities are matched at the OSBL base, only KPP is active in the OSBL. When diffusivities and viscosities are not matched, we assume that convective diffusivities are not applied in the OSBL. KPP and convective mixing target the same physical processes in the OSBL, so including both is double counting.

For either shape function, enhanced diffusivity can be included near the boundary layer base. This can smooth boundary layer deepening in time.

The non-local tracer flux uses the scalar shape function, scaled by the constant \(C_s\) from Eq. (20) of Large et al. (1994) rather than by \(h\, w_s\):

\[ \gamma_s(\sigma) = C_s\, S_1(\sigma). \]

Here \(S_1\) is always the unmatched scalar shape \(\sigma(1-\sigma)^2\), regardless of the MatchTechnique setting. Matching is a property of the diffusivity profile only: the matched shape is non-zero at \(\sigma = -1\) by construction, so reusing it for \(\gamma_s\) would leave a finite non-local flux at the boundary layer base that drops discontinuously to zero immediately below it. CVMix draws the same distinction, exposing the non-local shape as a separate choice from the matching option.

4 Design

4.1 Data types and parameters

4.1.1 Parameters

KPP is configured from the VertMix: KPP YAML group. Key parameters include:

• Enable

• CriticalBulkRichardsonNumber

• MatchTechnique (SimpleShapes or MatchBoth)

• InterpType2 (LMD94, Linear, Quadratic, Cubic)

• UseEnhancedDiffusion

• IceFractionThresholdForLangmuir

• IceFractionThresholdForMinimumOSBL

• MinimumOSBLUnderSeaIce

• DebugDiagnostics

Defaults and usage examples are documented in the user guide page: KPP in the User Guide.

4.1.2 Class/data structure

KPPMix is a singleton that owns persistent output fields and KPP-only workspaces, including:

• OSBLDepth, OSBLDepthIndex

• VertDiff, VertVisc

• VertNonLocalFlux

• diagnostics such as BulkRichardsonNumber, BulkRichardsonShear, UnresolvedShear, BuoyancyJump, and TurbulentVelocityScale

• surface-referenced pressure, reconstructed edge tangential velocity, and temporary ice fraction

4.2 Methods

The production interface is update(...), which receives canonical tracer, velocity, EOS, and forcing state, prepares KPP-derived state, validates required non-local tracer fluxes, and computes the KPP coefficients. It uses zero-copy temperature and salinity tracer subviews; the forcing state retains ownership of surface tracer fluxes.

The lower-level numerical stage is:

void computeKPPMix(const Array2DReal &PotentialDensity,
                   const Array2DReal &NormalVelocity,
                   const Array2DReal &TangentialVelocity,
                   const Array1DReal &SurfaceFrictionVelocity,
                   const Array1DReal &SurfaceBuoyancyFlux,
                   const Array2DReal &BruntVaisalaFreqSq,
                   const Array1DReal &IceFraction,
                   const Array1DReal &WindSpeed10m = Array1DReal());

Internal stages:

• computeOSBLDepth(...)

• computeMixingCoefficients(...), using surface forcing and optional interior mixing coefficients

4.3 Time stepper coupling behavior

KPP is coupled to all three Omega time steppers – Forward-Backward, RungeKutta2, and RungeKutta4. For every stepper, KPP is evaluated exactly once per time step, at the start of the step on the state at time \(n\), before any tendency is evaluated. The resulting boundary-layer depth, viscosity, diffusivity, and non-local flux profile are then held fixed for the remainder of the step.

This design differs from MPAS-Ocean, where the boundary layer depth, diffusivity, and viscosity are computed at the end of the time step and the non local flux is applied on the following timestep. The primary advantage of this new approach is:

1. Consistency. The non-local flux \(\gamma_s\) applied in the tracer tendency at each stage and the diffusivity \(K_s\) used by the end-of-step implicit vertical mixing solve are derived from the same OBL depth and the same shape function \(S_1(\sigma)\). Recomputing KPP at each stage would pair a stage-dependent \(\gamma_s\) with a different \(K_s\), breaking the correspondence described in section 3.2.

Because KPP is evaluated before the step advances, the KPP diagnostics written for a step describe the state at the beginning of that step. The coefficients are lagged relative to the state during the implicit solve at the end of the step. Full call-flow detail per stepper is described for developers and users in:

• Developer KPP workflow

• User runtime notes

5 Verification and Testing

5.1 Unit-level checks

The KPP unit-test driver is KPPMixTest.cpp. It supports focused test groups as well as an all group that runs all 28 unit tests:

• profiles (10 tests) checks the stability-function (used for the turbulent velocity scales) branches and continuity, simple shape and matched-shape functions, Langmuir utilities, OSBL depth clamping and indexing, turbulent velocity scales, and non-local flux scaling.

• osbl (7 tests) checks bulk-Richardson threshold crossing and quadratic interpolation, nonuniform and horizontally varying layer pseudo thickness, SSH-offset invariance to confirm positive definite geometric depths, partial columns, and area-weighted OSBL smoothing.

• vmix (8 tests) checks wind-only and convection-only profiles, stable and zero-forcing behavior, unmatched non-local profiles, MatchBoth interior coefficient matching, enhanced diffusion at the OSBL base, valid vertical domain edges, and invalid wet bounds.

• integration (3 tests) checks configuration parameters, the enabled end-to-end KPP call and diagnostics, and preservation of existing output fields when KPP is disabled.

The driver also provides three separate expected-failure cases: config-gradient, config-unsupported, and config-parabolic. These inject unsupported MatchTechnique values and verify that KPPMix::init() rejects them; they are not included in the all group. The test registration and group dispatch are shown at the end of KPPMixTest.cpp.

Together, these tests verify:

• OBL depth search threshold crossing, interpolation, diagnostics, and edge handling

• KPP profile, matching, enhanced-diffusion, and non-local-flux invariants

• Langmuir, sea-ice, smoothing, stable-forcing, and zero-forcing behavior

• KPP’s separation from shared VertMix background mixing and its handling of supplied interior coefficients

• Correct enabled, disabled, and invalid-configuration behavior

Tests cover requirements: 2.1, 2.2, 2.3, 2.4.

5.2 Polaris testing

The single column test case in Polaris has been extended to include a new set of regimes to exercise KPP across a set of stratifications and forcings. These tests are primarily drawn from Wagner et al., 2025 and Van Roekel et al., 2018, hereafter W25 and VR18 respectively.

The new tasks include

Test	Description
Wind-only, no Coriolis.	Tests mechanical shear-driven OSBL growth. Can be compared to solution from Kato and Phillips
Free convection (cooling)	VR18 FC case including a comparison to the analytic solution in Eq. F11
Mixed forcing convection	Wind, Evaporation, and cooling (all destabilizing)
Nonlocal supression	Stable heating with evaporation (VR18 CEW case). Tests suppression of nonlocal forcing under stable buoyancy forcing
Langmuir	Uses a fixed Stokes drift to test langmuir turbulence enabled and disabled
Sea ice	Wind forcing with a sea-ice fraction. Tests OSBL clamping under sea ice.
Convection with Evaporation	VR18 FCE case, cooling with evaporation. Tests the salinity contribution to surface buoyancy forcing.
Cooling with ML	VR18 FCML case with surface cooling on an existing ML.
Strong convection and stratification	W25 case showing KPP creating heating in the presence of cooling This is expected behavior

These tests include conservation and physics checks for each regime. These tests also include plots of critical quantities like the bulk Richardson number, ocean surface boundary layer depth, vertical diffusivity and viscosity, and the nonlocal flux.

Global stand-alone Omega cases, forced by annual averaged ERA-5 net surface heat, freshwater, and momentum fluxes provides a qualitative assessment of KPP behavior.