Discontinuous Galerkin discretization for two-equation turbulence closure model

TL;DR

This paper introduces a Discontinuous Galerkin discretization for two-equation turbulence models in ocean simulations, ensuring stability and positivity, and validates it through benchmarks and real-world river plume modeling.

Contribution

It presents a novel DG discretization method for turbulence closure models that maintains stability and positivity in ocean modeling applications.

Findings

Stable simulation of vertical mixing in stratified conditions

Accurate river plume characteristics captured

Validated with benchmarks and estuary simulation

Abstract

Accurate representation of vertical turbulent fluxes is crucial for numerical ocean modelling, both in global and coastal applications. The state-of-the-art approach is to use two-equation turbulence closure models which introduces two dynamic equations to the system. Solving these equations numerically, however, is challenging due to the strict requirement of positivity of the turbulent quantities (e.g., turbulence kinetic energy and its dissipation rate), and the non-linear source terms that may render the numerical system unstable. In this paper, we present a Discontinuous Galerkin (DG) finite element discretization of the Generic Length Scale (GLS) equations designed to be incorporated in a DG coastal ocean model, Thetis. To ensure numerical stability, the function space for turbulent quantities must be chosen carefully. In this work, we propose to use zeroth degree elements for the turbulent quantities and linear discontinuous elements for the tracers and velocity. The spatial discretization is completed with a positivity preserving semi-implicit time integration scheme. We validate the implementation with standard turbulence closure model benchmarks and an idealized estuary simulation. Finally, we use the full three-dimensional model to simulate the Columbia River plume. The results confirm that the coupled model generates realistic vertical mixing, and remains stable under strongly stratified conditions and strong tidal forcing. River plume characteristics are well captured.

Tables

Table 1. Table 1: List of parameter values for three different turbulence closure models. k − ε 𝑘 𝜀 k-\varepsilon and k − ω 𝑘 𝜔 k-\omega values (indicated by † † \dagger ) are from Tables 1 and 2 in Umlauf and Burchard ( 2003 ) ; The g ​ e ​ n 𝑔 𝑒 𝑛 gen model values (indicated by ∗ ∗ \ast ) corresponds to the first line in Table 7 in Umlauf and Burchard ( 2003 ) . In all cases, c l ​ i ​ m subscript 𝑐 𝑙 𝑖 𝑚 c_{lim} is computed with ( 36 ); c 3 − superscript subscript 𝑐 3 c_{3}^{-} , c μ 0 superscript subscript 𝑐 𝜇 0 c_{\mu}^{0} and σ Ψ subscript 𝜎 Ψ \sigma_{\Psi} are computed as detailed in B . The shown values correspond to the Canuto A stability functions.

Parameter	$k-ϵ$	$k-\omega$	$gen$
$p$	3†	-1.0†	2.0∗
$m$	1.5†	0.5†	1.0∗
$n$	-1.0†	-1.0†	-0.67∗
${\sigma }_k$	1.0†	2.0†	0.8∗
${\sigma }_{\Psi }$	1.20†	2.072†	1.18∗
$c_1$	1.44†	0.555†	1.0∗
$c_2$	1.92†	0.833†	1.22∗
$c_3^{+}$	1.0	1.0	1.0
$c_3^{-}$	-0.629	-0.643	0.052
$R_i^{st}$	0.25	0.25	0.25
$c_{\mu }^0$	0.5270	0.5270	0.5270
$\kappa$	0.4	0.4	0.4
$c_{lim}$	0.267	0.267	0.267
$k_{\text{min}}$	1.0$\times {10}^{-6}$	7.6$\times {10}^{-6}$	1.0$\times {10}^{-6}$
${\Psi }_{\text{min}}$	1.0$\times {10}^{-14}$	1.0$\times {10}^{-14}$	1.0$\times {10}^{-14}$

Table 2. Table 2: Prognostic and diagnostic variables and their function spaces. The shear and buoyancy frequencies, M 𝑀 M and N 𝑁 N are computed diagnostically from the state variables of the hydrodynamical model, and fed to the GLS model as inputs. Viscosity and diffusivity, ν 𝜈 \nu and ν ′ superscript 𝜈 ′ \nu^{\prime} , on the other hand, are fed to the hydrodynamical model.

GLS Model
Prognostic variables
Field	Symbol	Equation	Function space
Turbulent kinetic energy	$k$	(1)	${\text{P}}_0^{\text{DG}}\times {\text{P}}_0^{\text{DG}}$
Aux. turbulent variable	$\Psi$	(2)	${\text{P}}_0^{\text{DG}}\times {\text{P}}_0^{\text{DG}}$
Diagnostic variables
TKE dissipation rate	$\epsilon$	(11)	${\text{P}}_0^{\text{DG}}\times {\text{P}}_0^{\text{DG}}$
Turbulent length scale	$l$	(12)	${\text{P}}_0^{\text{DG}}\times {\text{P}}_0^{\text{DG}}$
Eddy viscosity	$\nu$	(13)	${\text{P}}_0^{\text{DG}}\times {\text{P}}_0^{\text{DG}}$
Eddy diffusivity	${\nu }^{\prime }$	(14)	${\text{P}}_0^{\text{DG}}\times {\text{P}}_0^{\text{DG}}$
Hydrodynamical model
Prognostic variables
Horizontal velocity	$𝒖$	(24-25) in TK2018	${\left[{\text{P}}_1^{\text{DG}}\times {\text{P}}_1^{\text{DG}}\right]}^2$
Temperature	$T$	(26) in TK2018	${\text{P}}_1^{\text{DG}}\times {\text{P}}_1^{\text{DG}}$
Salinity	$S$	(26) in TK2018	${\text{P}}_1^{\text{DG}}\times {\text{P}}_1^{\text{DG}}$
Diagnostic variables
Density	$\rho$	(14) in TK2018	${\text{P}}_1^{\text{DG}}\times {\text{P}}_1^{\text{DG}}$
Shear frequency	$M$	(5)	${\text{P}}_0^{\text{DG}}\times {\text{P}}_0^{\text{DG}}$
Buoyancy frequency	$N$	(6)	${\text{P}}_0^{\text{DG}}\times {\text{P}}_0^{\text{DG}}$

Equations

$$\frac{\partial k}{\partial t}+{\nabla }_h\cdot (\mathbf{u}k)+\frac{\partial (wk)}{\partial z}$$

$$\frac{\partial \Psi }{\partial t}+{\nabla }_h\cdot (\mathbf{u}\Psi )+\frac{\partial (w\Psi )}{\partial z}$$

$$P$$

$$B$$

$$M^2$$

$$N^2$$

$$\Psi ={\left(c_{\mu }^0\right)}^p{k}^m{l}^n,$$

$$c_3$$

$$\epsilon$$

$$l$$

$$\nu$$

$${\nu }^{\prime }$$

$$l_s$$

$$l_b$$

$$k_s$$

$$k_b$$

$${\left(\frac{\nu }{{\sigma }_k}\frac{\partial k}{\partial z}\right)|}_{{\Gamma }_b\cup {\Gamma }_s}$$

$${\Psi }_s$$

$${\Psi }_b$$

$$\frac{\partial {\Psi }_s}{\partial z}$$

$$\frac{\partial {\Psi }_b}{\partial z}$$

$${\left(\frac{\nu }{{\sigma }_{\Psi }}\frac{\partial \Psi }{\partial z}\right)|}_{{\Gamma }_s}$$

$${\left(\frac{\nu }{{\sigma }_{\Psi }}\frac{\partial \Psi }{\partial z}\right)|}_{{\Gamma }_b}$$

$${\left(\nu \frac{\partial \mathbf{u}}{\partial z}\right)|}_{{\Gamma }_b}$$

$$C_d$$

$${\left(\nu \frac{\partial \mathbf{u}}{\partial z}\right)|}_{{\Gamma }_s}$$

$$P+B-\epsilon$$

$$c_{1}P+c_{3}B-c_2\epsilon$$

$${\alpha }_M$$

$${\alpha }_N$$

$$R_i=\frac{N^2}{M^2}=\frac{{\alpha }_N}{{\alpha }_M},$$

$$R_i^{st}=\frac{c_{\mu }(R_i^{st})}{c_{\mu }^{\prime }(R_i^{st})}\frac{c_2-c_1}{c_2-c_3^{-}},$$

$$l\le l_{max}=c_{lim}\frac{\sqrt{2k}}{N}.$$

$$c_{lim}$$

$${\Psi }^{1/n}$$

$$\frac{\partial k}{\partial t}+A_k(k,\mathbf{u},w)$$

$$\frac{\partial \Psi }{\partial t}+A_{\Psi }(\Psi ,\mathbf{u},w)$$

$$S_k$$

Full text

††journal: Ocean Modelling

1 Introduction

Ocean models rely on parametrizations to account for sub-grid scale vertical mixing processes. The relevant eddy coefficients are obtained by means of turbulence closure models. While in some applications simple parametrizations, such as algebraic expressions of eddy viscosity and diffusivity, can be sufficient, in general more sophisticated schemes are needed in order to model the space-time evolution of turbulent fluxes.

Representing the dynamic evolution of turbulent fields is crucial especially in coastal, buoyancy-driven flows. Vertical mixing plays a major role in bottom boundary layer dynamics, evolution of stratification, and formation of the surface mixed layer, for example. It is also essential feature in estuarine circulation and river plume dynamics, which in general cannot be simulated without a sophisticated turbulence closure model.

Zero-equation parametrizations, such as Pacanowski and Philander (1981), the K-Profile Parametrization (KPP; Large et al., 1994; Van Roekel et al., 2018), and the recent ePBL parametrization (Reichl and Hallberg, 2018) are widely used in ocean models due to their simplicity and low computational cost. These methods do not include a prognostic turbulent variable but parametrize the eddy viscosity as a function of the mean flow state. Consequently, they are of limited applicability: the KPP and ePBL models only parametrize the the surface boundary layer, and none of the schemes can represent the temporal dynamics of turbulence.

The most sophisticated turbulence closures are two-equation models. These models consist of two partial differential equations, one for the turbulent kinetic energy (TKE), and another one for an auxiliary turbulent variable that defines the turbulent length scale. Such models include the Mellor-Yamada (level 2.5) model (Mellor and Yamada, 1982), $k-\varepsilon$ (Rodi, 1987), $k-\omega$ model (Wilcox, 1988), and the Generic Length Scale ($gen$; Umlauf and Burchard, 2003) model. The benefit of the GLS formulation is that all of the above closures can be obtained merely by changing parameters.

The turbulent eddy viscosity and diffusivity are obtained from the state variables, scaled by so-called (non-dimensional) stability functions. The most common stability functions are full-equilibrium functions by Canuto et al. (2001) (Canuto A and B), and by Cheng et al. (2002). Quasi-equilibrium functions, such as by (Kantha and Clayson, 1994), are also often used but have shown to be of limited applicability (Umlauf and Burchard, 2005).

Many existing ocean models implement two-equation turbulence closure models. The generic GLS model has become popular in recent years; it has been implemented in SELFE (Zhang and Baptista, 2008), SCHISM (Zhang et al., 2016), ROMS (Warner et al., 2005) and NEMO (Reffray et al., 2015) , for example, in addition to its original implementation in the GOTM library (Burchard et al., 1999). In most cases, the GLS equations have been implemented in a 1D vertical finite difference or finite volume context (e.g., GOTM); There have been only a few finite element implementations of two-equations turbulence closure models (e.g., Hill et al., 2012).

In this paper, we present a finite element discretization of the GLS equations intended to be incorporated in a Discontinuous Galerkin (DG) three-dimensional circulation model. The $k-\varepsilon$, $k-\omega$, and $gen$ models are tested with a series of standard benchmark test cases, and a realistic simulation of the Columbia River plume. We consider both the Canuto and Cheng stability functions. For the sake of simplicity, we do not consider convective adjustment methods.

The weak formulation of the GLS equations in a three-dimensional domain is presented. The function space for turbulent quantities must be chosen carefully to ensure numerical stability of the coupled turbulence–mean flow system (Burchard, 2002; Kärnä et al., 2012). We propose using a zeroth order DG space for the turbulent quantities as a natural choice for a linear DG hydrodynamical model; this choice also allows a straightforward implementation of the nonlinear source terms. The spatial discretization is completed by a positivity preserving, semi-implicit time integration scheme. Positivity of the turbulent quantities is ensured by using a $l$-stable implicit solver and Patankar treatment of the source terms. We present a time integration scheme for the full coupled turbulence–mean flow system.

The GLS model is implemented in the Thetis ocean model (Kärnä et al., 2018). Thetis is built on the Firedrake finite element modeling framework (Rathgeber et al., 2016) which uses a domain-specific language (Unified Form Language; Alnæs et al. 2014) to describe the weak forms; A just-in-time code generator is used produce computationally efficient C kernels at runtime. As such, Firedrake offers high flexibility in terms of describing the mathematical formulation without sacrificing computational efficiency. It also provides flexible support for various high performance computing platforms, as the code generator can apply hardware specific optimizations automatically. The present work demonstrates that turbulence closure models can be implemented in such domain specific language frameworks.

The GLS equations are presented in Section 2. The positivity preserving time integration scheme is presented in Section 3, followed by the finite element discretization in 4. Test cases are presented in Section 5, including a realistic application to the Columbia River plume in Section5.4, followed by Discussion and Conclusions in Sections 6 and 7.

2 Two equation turbulence closure models

2.1 Generic Length Scale equations

Let $\Omega$ be the three-dimensional domain in Cartesian coordinates $(x,y,z)$ (Figure 1 a). $\Omega$ spans from the sea floor $z=-h(x,y)$ to the free surface $z=\eta(x,y)$; the bottom and top surfaces are denoted by $\Gamma_{b}$ and $\Gamma_{s}$, respectively.

The Generic Length Scale (GLS) turbulence closure model (Umlauf and Burchard, 2003) solves the two equations for turbulent kinetic energy (TKE), $k$, and an auxiliary turbulent variable, $\Psi$,

[TABLE]

where $t$ denotes the time, $\bm{u}=(u,v,0)$ and $w$ are the horizontal and vertical velocity, respectively, and $\bm{\nabla}_{h}=(\partial/\partial x,\partial/\partial y,0)$ is the horizontal gradient operator. The vertical eddy viscosity is denoted by $\nu$; $\sigma_{k}$ and $\sigma_{\Psi}$ are the Schmidt numbers for the two state variables, respectively.

The second term in equations (1)-(2) represents horizontal advection. In many ocean models this term is neglected by assuming horizontal homogeneity of turbulence. In the coastal ocean, however, turbulence can change drastically in short distances, and this assumption is not valid. A river plume is a good example: Turbulent conditions change drastically across the plume front which separates the strongly stratified freshwater plume and weakly stratified coastal waters. As the front advances in the surface layer, failing to transport the turbulent quantities with the front would inevitably lead into erroneus mixing at the plume front.

The source terms in (1) are production of TKE, $P$; buoyancy production, $B$; and TKE dissipation rate $\varepsilon$. $P$ and $B$ depend on the vertical shear frequency $M$ and buoyancy (Brunt–Väisälä) frequency $N$, respectively:

[TABLE]

where $\nu^{\prime}$ is the eddy diffusivity of tracers, $g$ is the gravitational acceleration, $\rho$ is the density of the water, and $\rho_{0}$ denotes a constant reference density.

The term $P$ is always positive indicating that vertical velocity shear generates turbulence; $P$ transforms kinetic energy to turbulent kinetic energy (Burchard, 2002). The sign of $B$, on the other hand, depends on $N^{2}$: for stable stratification ($N^{2}>0$), $B$ is negative and transforms turbulent kinetic energy to potential energy. Stratification thus inhibits turbulence. In the case of unstable stratification, however, $B$ is positive and transforms potential energy to turbulence. Burchard (2002) showed that the numerical implementation of $B$ and $P$ should be energy conservative, such that increase of TKE by $P$, for example, matches the corresponding loss of kinetic energy. The dissipation term, $\varepsilon$, is always positive and represents loss of TKE into heat.

The auxiliary turbulent variable, $\Psi$, is defined by the parameters $p,m,n$:

[TABLE]

where $c_{\mu}^{0}$ is a dimensionless empirical parameter, and $l$ is a turbulent length scale. The source terms of $\Psi$ in (2) also depend on $P$, $B$, and $\varepsilon$, but are scaled by $\Psi/k$ and the empirical constants $c_{1}$, $c_{2}$, and $c_{3}$ which depend on the chosen closure.

The parameter $c_{3}$ controls buoyancy production of $\Psi$. It is split into two different values for stable and unstable stratification, respectively:

[TABLE]

The parameter $c_{3}^{+}$ is used to control mixing in unstably stratified conditions; for the closures considered herein $c_{3}^{+}=1.0$. The value of $c_{3}^{-}$ depends on the closure as discussed in B.1.

The TKE dissipation rate, $\varepsilon$, and turbulent length scale, $l$, are computed diagnostically from the state variables:

[TABLE]

Based on (11) it is clear that choosing $p=3,m=3/2,n=-1$ yields $\Psi=\varepsilon$ and results in the $k-\varepsilon$ model (Rodi, 1987).

Finally, the eddy viscosity and diffusivity are given by

[TABLE]

respectively, where $c_{\mu}$ and $c^{\prime}_{\mu}$ are non-dimensional stability functions that depend on $N^{2}$ and $M^{2}$. In this work, we consider three different stability functions, Canuto A and B, and Cheng, as defined in A.

The values of the empirical parameter $c_{\mu}^{0}$ and the Schmidt number $\sigma_{\Psi}$ are discussed in B.2 and B.3, respectively. The used parameter values are listed in Table 1.

2.1.1 Boundary conditions for $k$ and $\Psi$

Boundary conditions at the surface and bottom boundaries can be derived from law of the wall conditions. Within the boundary layer, the turbulent length scale grows linearly with the distance to the boundary, i.e.,

[TABLE]

where $\kappa=0.4$ denotes the von Karman constant and $z_{0,s}$ and $z_{0,b}$ are the corresponding roughness length scales. TKE in the (unresolved) boundary layer is (Burchard et al., 1999):

[TABLE]

where $u_{s}^{*}$ and $u_{b}^{*}$ are the friction velocities. From (17)-(18) the Neumann condition for $k$ can be derived,

[TABLE]

which is used as the boundary condition for $k$. Note that this is a natural boundary condition in DG finite element formulation as the (diffusive) boundary term vanishes if $\partial k/\partial z=0$.

Substituting (15)-(16) to (7) yields

[TABLE]

Differentiating with respect to $z$ gives

[TABLE]

In practice, the gradients cannot be evaluated right at the boundary because $\Psi$ can be changing rapidly. We evaluate the boundary conditions at distance $\Delta z/2$ from the boundary where $\Delta z$ is the vertical element size. This choice is consistent with most finite volume implementations, and results in stable behavior. The $\Psi$ boundary conditions then become

[TABLE]

2.1.2 Boundary conditions for the momentum equation

Regarding the momentum equation (see Kärnä et al. 2018), we impose the log-layer bottom friction condition on the bottom boundary:

[TABLE]

where $C_{d}$ is the drag coefficient, $z_{0}$ is the bottom roughness length, $z_{b}$ is the $z$-coordinate at the middle of the bottom most element, and $\bm{u}_{b}=\bm{u}(z_{b})$.

On the free surface, wind stress ($\bm{\tau}_{w}$) is imposed:

[TABLE]

In applications with atmospheric wind forcing, the stress is computed using the formulation by Large and Pond (1981), and a constant air density, $\rho_{air}=1.22\ \text{kg}\ \text{m}^{-3}$.

For tracers, impermeable boundary conditions, with zero advective and diffusive fluxes, are imposed at the surface and bottom boundaries.

2.2 Steady state solution

In a steady state, and assuming homogeneous turbulence, the source terms of the governing equations (1)-(2) balance each other resulting in

[TABLE]

It is convenient to define non-dimensional shear and buoyancy frequencies (Burchard and Bolding, 2001),

[TABLE]

Using the gradient Richardson number,

[TABLE]

the equilibrium condition (29)-(30) can be expressed as (Umlauf and Burchard, 2003)

[TABLE]

where $c_{\mu}$ and $c^{\prime}_{\mu}$ are the stability functions, and $R_{i}^{st}$ denotes the steady state gradient Richardson number (commonly set to value 0.25).

2.3 Length scale limitation

Galperin et al. (1988) suggested to limit $l$ as

[TABLE]

Umlauf and Burchard (2005) show that the limiting factor $c_{lim}$ is a function of the steady-state Richardson number:

[TABLE]

The length scale limitation can be formulated as a limit on $\Psi$ (Warner et al., 2005)

[TABLE]

Note that for negative $n$, (37) imposes a lower limit on $\Psi$. In this work we limit $\Psi$ only, i.e., no limit is applied on $l$ or $\varepsilon$.

3 Temporal discretization

This section outlines the temporal discretization before introducing the DG discretization in Section 4.

3.1 Thetis coastal ocean model

The turbulence closure model was implemented in the Thetis three-dimensional circulation model (Kärnä et al., 2018). Thetis solves the hydrostatic equations with a semi-implicit DG finite element method. The equations are solved in a time-dependent 3D mesh, unstructured in the horizontal direction. The mesh moves in the vertical direction to track the free surface undulation; the mesh movement is implemented with the Arbitrary Lagrangian–Eulerian (ALE) method. A second-order split-implicit time integration scheme is used to advance the equations in time. Advection of 3D fields is solved with an explicit Strong Stability-Preserving (SSP) scheme. Vertical diffusion is treated implicitly. The model formulation is described in detail in Kärnä et al. (2018).

For 3D advection, Thetis uses a two-stage, second-order SSP Runge-Kutta scheme, SSPRK(2,2) (Shu and Osher, 1988). The SSP property ensures that the solution is non-oscillatory as long as the CFL condition is satisfied. With DG spatial discretization, this implies that the element mean value is positive definite; overshoots can still arise within an element. To this end, a vertex-based slope limiter (Kuzmin, 2010) is used to redistribute mass in the element if local overshoots are present. The slope limiter is applied as a post-processing step after solving the equations. Kärnä et al. (2018) demonstrate that the model is mass conservative (typical error being $\mathcal{O}(10^{-12})$) and non-oscillatory (largest overshoots are of order $\mathcal{O}(10^{-5})$ in extreme cases).

In the GLS implementation, the turbulent quantities are advanced in time similarly to other 3D fields: Horizontal and vertical advection is solved with the same SSPRK ALE scheme. Vertical dynamics of $k$ and $\Psi$ are solved implicitly.

Thetis source code is freely available. We have also archived the exact source code version used to to produce the results in this paper (see C).

3.2 Implicit and explicit terms

The governing equations (1)-(2) are advanced in time with a fractional step method. We first separate the explicit advection terms, the remaining implicit terms:

[TABLE]

where $A_{i},\ i=\{k,\Psi\},$ denote the horizontal and vertical advection terms, $D_{i}$ the vertical diffusion terms, and $S_{i}$ the source terms:

[TABLE]

Let $\Delta t$ denote the time step. Given the state variables at time level $n$, the time update of $k$ can be divided into an explicit and implicit steps ($\Psi$ is treated analogously):

[TABLE]

The explicit advection update (42) is solved with the ALE SSPRK method as described in Kärnä et al. (2018). The implicit update (43) is presented below. Note that $D_{k}$ and $S_{k}$ have been linearized with respect to the prognostic variable: only $k$ is taken at time level $n+1$.

To ensure mass conservation and positivity preservation of the whole scheme, both stages (42) and (43) must satisfy these properties. As stated above, the 3D advection step is conservative and positivity preserving. The implicit scheme, however, must be designed with care: The diffusion operator may introduce oscillations if viscosity (or time step) is very large. Moreover, the non-linear source terms ($S_{k}$ and $S_{\Psi}$) may cause negative values unless special treatment is applied.

3.3 Implicit solver for the GLS equations

The prognostic variables $k$ and $\Psi$ are positive by definition. The discretized system must maintain positivity of these variables at all times to ensure physically sound and numerically stable solution.

As the vertical problem is generally stiff, it is crucial to use an L-stable implicit method to avoid spurious oscillations (Alexander, 1977; Hairer and Wanner, 1996). L-stability means that the amplification factor of the scheme tends to zero as the time step tends to infinity, effectively dampening unresolved high-frequency oscillations. It is worth noting that the commonly-used Crank-Nicolson method is not L-stable and may lead to erratic and unphysical solutions. The simplest L-stable choice is the Backward Euler method (43). It is used in Kärnä et al. (2018) and also employed in this work. Alternatively, a more accurate multi-stage, L-stable Diagonally Implicit Runge Kutta scheme could be used (e.g., Ascher et al., 1997). Practical evaluations, however, showed that using a two-stage second order method did not affect the results significantly but did increase the cost of the implicit solve by a factor of two.

In Thetis, the domain is decomposed into parallel sub-regions in the horizontal direction, i.e., each vertical column of elements lies within a single process. Because the implicit solve only contains vertical fluxes, the implicit solve is local and does not involve parallel communication.

3.4 Patankar treatment of source terms

The non-linear source terms must be treated carefully because they can trigger instabilities and/or generate negative values. Positivity preservation is based on the Patankar treatment where all sink terms (i.e., $S_{i}<0$) are treated implicitly (Patankar, 1980; Burchard et al., 2003) to ensure that they cannot overshoot to produce negative quantities during temporal integration. A similar approach is used in GOTM, for example.

In the $k$ equation, the shear production term, $P$, and TKE dissipation rate, $\varepsilon$, are always positive, resulting in a source and sink term, respectively. The sign of the buoyancy production term, $B$, on the other hand, depends on the sign of $N^{2}$. We therefore split the term in two parts: $B=B_{+}+B_{-}$ with $B_{+}\geq 0$ and $B_{-}\leq 0$:

[TABLE]

In the $\Psi$ equation, we have $c_{1}P\geq 0$ and $c_{2}\varepsilon\geq 0$. In this case the sign of the $B$ term, however, also depends on the parameter $c_{3}$: $c_{3}^{+}>0$ is used when $B\geq 0$, and $c_{3}^{-}$ when $B<0$ ($c_{3}^{-}$ can be negative). Denoting $c_{3}B=c_{3}^{+}B_{+}+c_{3}^{-}B_{-}$, we define:

[TABLE]

We can now split the source terms in (43) to explicit and implicit parts,

[TABLE]

Note that in $S_{k}^{-}$ we have employed the Patankar treatment, i.e., scaled the sink term with the ratio $k^{n+1}/k^{*}$ to render it implicit (here $k^{*}$ is the “old” value due to the fractional stepping).

The update for $\Psi$ is analogous with the source terms:

[TABLE]

3.5 Setting minimum values for $k$ and $\varepsilon$

The presented solver ensures that $k$ and $\Psi$ remain positive throughout the simulation. However, due to round-off errors small negative values can appear which lead into undefined values (e.g., due to the square root in (12)). The turbulent quantities are, therefore, cropped to a small positive minimum value; this is a common practice in several implementations (e.g., Umlauf et al., 2005; Warner et al., 2005). The minimum values for $k$ and $\Psi$ are listed in Table 1. In addition we impose minimum values $l_{min}=1.0\times 10^{-12}\ \text{m}$ and $\nu_{min}=\nu^{\prime}_{min}=1.0\times 10^{-8}\ \text{m}^{2}\ \text{s}^{-1}$. It should be noted that the user can modify these limits and also introduce additional background viscosity/diffusivity if needed. These limits are imposed after every update.

4 Spatial discretization

4.1 Mathematical notation

The domain $\Omega$ is divided into three-dimensional elements $e\in\mathcal{P}$. The mesh is generated by extruding a two-dimensional surface mesh over the vertical dimension. The surface mesh consist of either triangles or quads, resulting in triangular prisms, or hexahedral 3D elements, respectively.

We denote discontinuous Galerkin function spaces of degree $p$ by $\text{P}^{\text{DG}}_{p}$. A function $f\in\text{P}^{\text{DG}}_{p}$ is a polynomial of degree $p$ in each element and discontinuous at the element interfaces. On the extruded mesh, we define function spaces as a Cartesian product of the horizontal (2D) and vertical (1D) elements: $\text{P}^{\text{DG}}_{p}\times\text{P}^{\text{DG}}_{q}$; a function in this space is therefore a polynomial of degree $p$ in the horizontal direction and degree $q$ in the vertical direction.

The set of horizontal and vertical interfaces are denoted by $\mathcal{I}_{h}$ and $\mathcal{I}_{v}$, respectively (Figure 1 b). The outward unit normal vector is denoted by $\mathbf{n}=(n_{x},n_{y},n_{z})$; its restriction on the horizontal and vertical interfaces are denoted by $\mathbf{n}_{h}$, and $\mathbf{n}_{v}$, respectively. The vertical facets $\mathcal{I}_{v}$ are strictly vertical, implying $\mathbf{n}_{v}=(n_{x},n_{y},0)$.

In the weak forms we use the following notation for volume and interface integrals

[TABLE]

In interface terms we additionally use the average $\{\!\!\{\cdot\}\!\!\}$ and jump $[\![\cdot]\!]$ operators for arbitrary scalar $a$ and vector $\mathbf{u}$ fields:

[TABLE]

where the superscripts ’$+$’ and ’$-$’ arbitrarily label the values on either side of the interface. Note that for the outward normal vectors it holds $\mathbf{n}^{+}=-\mathbf{n}^{-}$.

4.2 Choosing function spaces

The prognostic variables of the GLS system (1)-(2) are $k$ and $\Psi$. Diagnostic variables include the dissipation rate $\varepsilon$, turbulent length scale $l$, vertical eddy viscosity $\nu$, and diffusivity $\nu^{\prime}$. Diagnostic variables originating from the 3D circulation model are the vertical shear frequency $M$ and buoyancy frequency $N$. The choice of the function spaces where these variables reside is crucial for numerical stability and accuracy.

The production terms $P$ and $B$ are of crucial importance as they can excite spurious modes in the numerical solution. Because $P=\nu M^{2}$, the function spaces of $\nu$ and $M$ determine the polynomial degree of $P$: for example, if $\nu$ and $M$ are linear, $P$ is nominally cubic (but could be projected to a lower function space).

For numerical stability it is natural to choose function spaces that can satisfy the steady state condition $P+B-\varepsilon=0$ point-wise (i.e., in a strong sense instead of an integral, weak sense). This implies that $P$, $B$, and $\varepsilon$ should belong to the same function space. As we can have $\Psi=\varepsilon$ it is convenient to choose that $k$ and $\Psi$ also belong to the same space. This is also desirable because it ensures that the source terms can be directly mapped to the space of the prognostic variables, and cannot thus excite spurious modes.

Thetis uses linear discontinuous Galerkin elements for 3D tracers and all velocity components. That is, the variables belong to a function space, $\mathcal{W}_{1}=\text{P}^{\text{DG}}_{1}\times\text{P}^{\text{DG}}_{1}$. Consequently, $\partial u/\partial z$ belongs to $\text{P}^{\text{DG}}_{0}$ space in the vertical direction. Noting that a product of two functions in $\text{P}^{\text{DG}}_{0}$ still belongs to the same space, $\text{P}^{\text{DG}}_{0}$ is a natural choice for the shear frequency $M$ as well. The same reasoning holds also for the tracer gradients and squared buoyancy frequency $N^{2}$ if one assumes a linear equation of state.

The above requirements can be satisfied by choosing that the turbulent quantities belong to the degree zero DG space, $\mathcal{W}_{0}=\text{P}^{\text{DG}}_{0}\times\text{P}^{\text{DG}}_{0}$. With this choice we in fact have $f,g\in\mathcal{W}_{0}\Rightarrow fg,f/g,f\sqrt{g}\in\mathcal{W}_{0}$ which implies that the stability functions and diagnostic variables can be computed point-wise. The function spaces are summarized in Table 2.

It would be desirable to use the same $\mathcal{W}_{1}$ space for $k$ and $\Psi$ as well, but (at least in the case of $k-\varepsilon$ model) this choice violates the $P=\varepsilon$ steady state condition and leads to unstable solutions. Finding higher order function spaces for the turbulent quantities remains an open question.

4.3 Discontinuous Galerkin discretization of the GLS equations

Let $\phi\in\mathcal{W}$ be a test function. The weak formulations are derived by taking the continuous equations (1-2), multiplying them by $\phi$, and integrating over the domain $\Omega$. The weak formulation then reads: find $k,\Psi\in\mathcal{W}$ such that

[TABLE]

where $\mathcal{A}_{h,i}$ and $\mathcal{A}_{v,i}$, $i=k,\Psi$ denote the horizontal and vertical advection terms, while $\mathcal{D}_{i}$ and $\mathcal{S}_{i}$ are the diffusion and source terms, respectively.

The advection terms are discretized with a standard first order DG upwind technique as implemented in the hydrodynamical model (Kärnä et al., 2018). The diffusion terms are treated with Symmetric Interior Penalty Galerkin method (SIPG; Epshteyn and Rivière, 2007) with a penalty factor $\sigma=1/L$ where $L$ is the vertical element size (Kärnä et al., 2018).

The source terms are split to sources and sinks as shown in Section 3.4:

[TABLE]

Let $\phi_{i}$ be the basis of the function space $\mathcal{W}$. We require that the solution belongs to this space, i.e., $k,\Psi\in\mathcal{W}$, and use the finite element approximations (e.g., $k=\sum_{i}k_{i}\phi_{i}$). Using the basis functions as test functions, we can write (53-54) in a bilinear form:

[TABLE]

where [\mathbf{M}]_{i,j}=\Big{\langle}\phi_{i}\phi_{j}\Big{\rangle}_{\!\!\Omega} is the mass matrix, and the other terms are also bilinear, e.g, $\mathbf{D}_{k}=\mathbf{D}_{k}(\phi_{i},\phi_{j})$.

4.4 Computing shear and buoyancy frequencies

Computing the vertical shear and buoyancy frequencies accurately is crucial for maintaining numerical stability (Kärnä et al., 2012). We compute the vertical gradients weakly. Denoting the shear frequency vector by $\bm{M}=\frac{\partial\bm{u}}{\partial z}$ we define

[TABLE]

Analogously, we can compute the $\frac{\partial\rho^{\prime}}{\partial z}=R_{z}$

[TABLE]

4.5 Time integration

For convenience we re-write the $k$ and $\Psi$ weak forms as

[TABLE]

where $F_{k}$ and $F_{\Psi}$ contain the advection terms whilst $G_{k}$ and $G_{\Psi}$ contain the vertical diffusion and source terms.

The solution can then be split in explicit advection stage and implicit vertical dynamics stage. The explicit stage is solved with SSPRK(2,2) scheme (Kärnä et al., 2018):

[TABLE]

where $w_{m}$ denotes the vertical mesh velocity and $\tilde{F}$ denotes the ALE weak forms of $F_{k}$ terms (see Kärnä et al. (2018)). The advection stage for $\Psi$ is derived analogously.

Using a linearized Backward Euler update, the implicit stage reads:

[TABLE]

The full time integration scheme is summarized in Algorithm 1.

5 Results

We first verify the GLS implementation with a series of standard turbulence closure test cases, followed by a realistic simulation of the Columbia River plume. All the tests have been carried out with the $k-\varepsilon$, $k-\omega$, and the $gen$ models as described in Table 1.

5.1 Bottom friction

The first test is an open channel with bottom friction. The fluid is initially at rest, forced only by a constant free surface slope $\partial\eta/\partial x=-1\times 10^{-5}$. In the absence of rotation, the flow converges to a steady state solution where the pressure gradient is balanced by the bottom friction.

Near the bed the flow velocity follows the well-known logarithmic profile, which can be expressed as (e.g., Hanert et al., 2007)

[TABLE]

where $\bm{u}^{*}_{b}$ is the bottom friction velocity and $z_{0}^{b}=1.5\times 10^{-3}\ \text{m}$ is the bottom roughness length. The bottom friction velocity can be solved from the momentum balance:

[TABLE]

The eddy viscosity follows a parabolic profile

[TABLE]

The free flow was simulated in a 5.0 km wide square domain with constant 15 m depth. Horizontal mesh resolution was 2.5 km. Periodic boundary conditions were used in the $x$ direction. Initially the flow is at rest. A steady state solution is reached after roughly 12 h of simulation. All simulations were carried out with 25 s time step.

First, we experimented the three closure models, GLS, $k-\varepsilon$, and $k-\omega$ with Canuto A stability functions. Vertical profiles of velocity, TKE, $\varepsilon$, $l$, and viscosity are shown in Figure 2. In all simulations 250 vertical levels were used, resulting in 6 cm vertical resolution. The parabolic Courant number, $\nu\Delta t/(\Delta z)^{2}$, was approximately 300. The velocity profile is close to the analytical logarithmic solution in all cases. TKE shows the expected nearly linear profile. In terms of eddy viscosity, the $k-\omega$ model matches best to the analytical solution. GLS and $k-\varepsilon$ models tend to overestimate viscosity in the upper water column, while underestimating it in the middle. This behavior is in-line with other results in the literature (Warner et al., 2005). One possible cause for this deviation is the surface boundary condition for $\Psi$, which affects both $\Psi$ and $l$ near the free surface (Figure 2 d).

Second, the $k-\omega$ model was run with different vertical grids consisting of 5, 25, and 250 levels (Figure 3). The corresponding parabolic Courant numbers are 0.1, 3, and 300. The velocity profiles are close to the analytical solution for all the resolutions. As resolution is increased viscosity converges to a parabola close to the analytical solution. This test demonstrates that the effective friction felt by the water column does not depend strongly on the vertical resolution. In addition, it shows that the numerical solver remains stable even with small vertical elements (6 cm).

Third, we experimented with the different stability functions. We ran the $k-\omega$ model with Canuto A, B and Cheng stability functions (Figure 4). In general the choice of the stability function does not have a major impact on the results. The results of Canuto A and Cheng stability functions are nearly identical. Canuto B functions, on the other hand, yields slightly lower TKE, but the difference in viscosity and velocity is negligible. The effect of stability functions was similar with other closure models as well (not shown).

5.2 Wind-driven entrainment

The next test examines mixed layer deepening due to surface stress, based on the laboratory experiment originally conducted by Kato and Phillips (1969). Initially the water column is at rest and linearly stratified. A constant surface stress is applied at the free surface. Stress induced mixing leads to a homogeneous surface layer that grows deeper in time. The depth of the mixed layer follows the empirical formula suggested by Price (1979):

[TABLE]

where $u^{*}_{s}$ is the surface friction velocity and $N_{0}$ is the initial spatially invariable buoyancy frequency. Here values $u^{*}_{s}=0.01\ \text{m}\ \text{s}^{-1}$ and $N_{0}=0.01\ \text{s}^{-1}$ are used. $N_{0}$ is prescribed by imposing a suitable linear salinity field and using a linear equation of state.

The mixed layer deepening was simulated in 5 km wide square domain with constant 50 m depth and 2.5 km mesh resolution. All simulations were carried out with 30 s time step unless otherwise mentioned.

Six different mesh resolutions, $\Delta z=0.2,0.5,1.0,2.0,3.125$, and $5.0\ \text{m}$, were experimented with, resulting in 250, 100, 50, 25, 16, and 10 vertical levels, respectively. The corresponding parabolic Courant number ranges from 7 to 0.03.

The evolution of the mixed layer depth is presented in Figure 5 for $k-\varepsilon$ model and Canuto A stability functions. All resolutions result in mixed layer deepening that is in good agreement with the empirical formula (72). With low vertical resolution, the mixed layer depth advances in steps as mixing penetrates new elements.

Next we tested the three different turbulence closure models (with Canuto A stability functions) on the medium resolution $\Delta z=0.5\ \text{m}$ mesh (Figure 6). All models reproduce a realistic mixed layer depth. The $k-\varepsilon$ and GLS models are close to the empirical mixed layer depth. The $k-\omega$ model tends to under estimate the mixed layer depth, especially in the beginning of the simulation, but the difference is small. These results are in good agreement with previous studies (e.g., Warner et al. 2005; Kärnä et al. 2012).

We also carried out experiments with different stability functions (not shown); in all cases, the mixed layer depth behaved correctly, its variance being similar to what is seen in Figure 6.

The above tests were carried out with a short time step (30 s) close to typical values in regional applications. To test how the model can handle stiff problems we increased the time step in the fine grid ($\Delta z=0.2\ \text{m}$) case up to $\Delta t=1200\ \text{s}$ resulting in parabolic Courant number of 300 (Figure 7). In all the presented cases, the model remains stable and produces a deepening mixed layer. With long time steps ($>300\ \text{s}$), however, the deepening of the mixed layer is underestimated especially during the first 10 h of the simulation. These results suggests that the scheme is indeed stable for stiff problems but starts to dissipate rapid processes (which occur mainly on the onset of the mixed layer in the beginning). The observed behavior is consistent with results by Reffray et al. (2015) where the NEMO 1D model also began to underestimate the mixed layer depth for very long time steps.

5.3 Idealised estuary

Estuarine circulation is an essential feature in coastal domains. It is dominated by the interplay of buoyancy driven stratification and vertical mixing. Buoyancy input from a river tends to tilt isopycnals and increase stratification in the estuary. On the other hand, bottom friction and up-estuary currents tilt the isopycnals in the opposite direction causing unstable stratification in the bottom layer. This generates vigorous mixing at the bed which, under suitable conditions, can mix the entire water column. Under sufficiently strong stratification, however, turbulence cannot penetrate the pycnocline resulting in a well mixed bottom layer that is de-coupled from the fresher surface layer. The competition between these two tendencies, and temporal evolution of the turbulence fields, cannot be modeled without a dynamic two-equation turbulence closure model.

We test the model’s capability of representing estuarine circulation with the idealized estuary test case by Warner et al. (2005). This test case has been used to validate other models as well, e.g., SLIM (Kärnä et al., 2012), and SELFE (Lopez and Baptista, 2017).

The idealized estuary test case does not have an analytical solution. Furthermore, our tests suggest that the numerical solution is highly sensitive to the open boundary conditions which vary in each model due to its discretization. As such, carrying out a direct comparison of the modeled fields is not beneficial. Nevertheless, this is a useful dynamical test to verify that the coupled model can generate realistic estuarine circulation features; It also reveals possible numerical instabilities associated with strong horizontal gradients, advection, and tidal currents.

The domain is a rectangular channel 100 km long and 1 km wide. Depth varies linearly from 10 m at the deep (ocean) boundary to 5 m in the shallow (river) end. Horizontal mesh resolution is 500 m; 20 equally distributed vertical levels are used.

At the river boundary, a constant discharge of $F_{river}=400\ \text{m}^{3}\ \text{s}^{-1}$ is imposed, corresponding to a 8 $\text{cm}\ \text{s}^{-1}$ seaward velocity. Water elevation is unprescribed, and water salinity is set to zero. At the ocean boundary, a sinusoidal tidal current is imposed with $0.4\ \text{m}\ \text{s}^{-1}$ amplitude and period $T_{tide}=12\ \text{h}$. In addition, we superimpose the (seaward) river flux, $F_{river}$, at ocean boundary as well to ensure tidally averaged volume conservation. Salinity of in-flowing ocean water is set to is 30 psu. Throughout the simulation temperature is kept at constant $10\ ^{\circ}\text{C}$.

Initially water is at rest. Salinity varies from 30 to 0 psu along the channel between 30 and 80 km from the ocean boundary (Figure 8 a). The ocean and river boundary conditions are ramped up for one hour. During the simulation estuarine circulation quickly develops, driving freshwater seaward in the surface layer. A thick salt wedge forms, oscillating with the tides. The model is spun up for 15 days after which the solution has reached a quasi-periodic state where each tidal cycle is nearly identical.

Figure 8, shows the instantaneous salinity distribution after 16 days of simulation for $k-\varepsilon$ (b), $k-\omega$ (c) and GLS models (d). In all simulations the Canuto A stability functions were used. The results show a thick salt wedge, and a fairly thin, seaward flowing surface layer. Salinity reaches about 80 km from the estuary mouth. The $k-\omega$ model shows lower stratification and shorter salinity intrusion indicating stronger effective mixing. The GLS model results in the largest salinity intrusion.

Vertical profiles are presented in Figure 9 for the three closures. The profiles indicate that the well-mixed bottom layer fills most of the water column. In the upper water column stratification is very strong which limits mixing and results in almost turbulence-free surface layer. Salinity profiles in Figure 9 (a) show that stratification is indeed weaker in the $k-\omega$ model. It also generates higher TKE and viscosity in the bottom layer suggesting a larger overall (tidally-averaged) mixing. The $k-\varepsilon$ and GLS model behave quite similarly. GLS, however, generates lower TKE and viscosity which is consistent with the larger salinity intrusion length (Figure 8).

Qualitatively the presented results are similar to those in Warner et al. (2005), Kärnä et al. (2012) and Lopez and Baptista (2017). Some differences do exist in salt intrusion and stratification, for example, which could be attributed to differences in the ocean boundary conditions and, more generally, the model’s discretization. It should also be noted that the behavior of the three turbulence closure models strongly depends on the chosen parameters.

5.4 Application to Columbia River plume

The Columbia River is a major freshwater source in the west coast of the United States (Figure 10). With an annual mean discharge of 5500 $\text{m}^{3}\ \text{s}^{-1}$ it is the second largest river in continental USA (Chawla et al., 2008). Maximal daily tidal range varies between from less than 2.0 to 3.6 m at the river mouth (Kärnä et al., 2015). The river plume is a major feature in the Pacific Northwest Coast, affecting stratification, currents, and nutrients in the coastal waters (Hickey and Banas, 2003).

The river plume can be divided into different water masses. Horner-Devine et al. (2009) define four water masses: the source water, followed by the tidal, re-circulating, and far-field plume. The tidal plume (with 6 to 12 h time scale) consists of a freshwater lens being emitted from the river mouth every tidal cycle. Once released, it rapidly spreads out and becomes thinner, its depth ranging from 6 to 3 m. Salinity is below 21 psu. The tidal plume eventually merges with the re-circulating plume.

The re-circulating plume forms an anticyclonic bulge at the vicinity of the river mouth. In the cross-shore direction, the bulge is roughly 40 km in diameter. The re-circulating plume consists of freshwater volume equivalent to 3–4 days of river discharge; salinity is between 21 and 24 psu.

The far-field plume lies beyond the re-circulating plume. It is the zone where final mixing of the plume and ambient ocean waters takes place, unaffected by the momentum of the river discharge. Typically, the far-field plume forms a large coastal current, driven by buoyancy and Earth’s rotation, that can extend hundreds of kilometers North from the river mouth. The evolution of the far-field plume, however, is strongly affected by winds and coastal currents.

The typical northward propagation of the river plume can be arrested by northerly, upwelling-favorable winds. Upwelling generates an off-shore current in the surface layer, eroding the river plume from the coast. The plume detaches from the coast, spreads out and becomes thinner. Under sufficiently strong and prevailing winds, the plume begins to travel southwest. Northerly winds are more common in the Pacific Northwest during summer months, which is why the southwest-oriented plume is often called the summer plume. Upwelling-favorable winds, therefore, can reverse the orientation of the plume.

Southerly, downwelling-favorable winds, on the other hand, enhance the northward coastal current. Under downwelling conditions, the surface layer flows onshore and pushes the plume against the coast, making it narrower and thicker. Consequently, the coastal current becomes more pronounced: it becomes narrower and propagates faster towards the North.

The evolution of the river plume was simulated with Thetis. The computational domain spans from roughly 40 ∘N to 50 ∘N in the zonal direction (Figure 10). The western boundary is located roughly 100 km west of the coast. The Columbia River is included in the domain, excluding the shallow lateral bays. The river boundary is located at Beaver Army, 86 km upstream from the mouth. The horizontal mesh resolution (triangle edge length) ranges from roughly 900 m in the estuary to 15 km at the open ocean boundaries. In the vertical direction, 20 terrain-following layers were used. The layer thickness ranges from roughly 1 m near the surface to a maximum of 180 m at the bottom. The entire vertical grid adjusts to the free surface similarly to $z^{*}$ grid used in structured grid models (Adcroft and Campin, 2004).

We used a composite bathymetry data (Kärnä and Baptista, 2016a) that was interpolated on the model grid (Figure 10). In addition, the bathymetry was smoothed with a Laplacian diffusion operator to limit strong gradients which could introduce internal pressure errors. Minimum water depth was set to 3.5 m. Wetting and drying was neglected.

Atmospheric pressure and 10 m wind were obtained from the North American Mesoscale Model (NAM) model. Wind stress was computed with the formulation by Large and Pond (1981). At the ocean boundary temperature, salinity and sub-tidal velocity was imposed from the global Navy Coastal Ocean Model (NCOM) model (Barron et al., 2006). In addition, tidal elevation and velocities were imposed from the TPXO v9.1 model (Egbert and Erofeeva, 2002). At the river boundary discharge and water temperature are imposed using US Geological Survey (USGS) gauge data. Salinity is set to zero.

The simulation covers a time period from May 1 2006 to July 2. The first month is used to spin up the simulation; the analysis period is from June 1 onward. This spin-up time is sufficient for the estuary–plume system as the estuary has a residence time of a few days (Kärnä and Baptista, 2016b) and the far-field plume develops in a couple of weeks.

It is worth noting that the Thetis model setup has not been fully calibrated for the system and the results are therefore preliminary. An exhaustive model calibration and skill analysis are out of the scope of the present article and will be addressed in the future.

The simulations were carried out on a Linux cluster with 16-core Intel Xeon E5620 processors and Mellanox Infiniband interconnect. 96 MPI processes were used. The simulation had $2\times 10^{6}$ tracer degrees of freedom (DOF) and used 15 s time step; the full 90-day simulation took 7.2 days.

Figure 11 shows the atmospheric and fluvial forcing conditions for the analysis period. Winds alternate between up- (northerly) and downwelling-favorable direction (Figure 11 a). Tides exhibit a typical spring-neap variation to the system, tidal range varying between 1.6 and 3.4 m; River discharge is roughly 10000 $\text{m}^{3}\ \text{s}^{-1}$ during the simulation (Figure 11 b and c).

Simulated salinity field is compared against an observational data set from the RISE project (River Influences on Shelf Ecosystems, e.g., MacCready et al. 2009). MacCready et al. (2009) and Liu et al. (2009) use data from the 2004 summer measurement campaign. It should be noted that for the 2006 campaign, used in the present analysis, the north and south moorings were moved further out from the river mouth by some 80 km (see Figure 10 for the mooring locations). Therefore, the 2006 data set is well-suited for evaluating the behavior of the coastal current and the far-field plume under changing wind conditions.

Figure 12 shows time series of observed surface salinity at stations rino, rice, ogi01, and riso. The rice station (Figure 12 b) is located in the tidal plume recording the passing of each tidal freshwater lens. The model replicates sub-tidal and tidal evolution variability of the near field plume well. Note that exact replication of the tidal signal is difficult as salinity gradients are very large. Measured salinity therefore strongly depends on the passing of the plume fronts, which is sensitive to atmospheric forcing, mesh resolution, mixing, and bathymetry features, for instance.

The ogi01 station (Figure 12 c) is located some 25 km southeast from the river mouth. It captures the re-circulating plume when it is pushed southward from the river mouth. Typically this occurs during southerly winds. During the simulation period, two of such events are recorded, on June 6–10 and June 18–20. The model captures the southward transport of the plume relatively well, especially during the latter event.

The northern station, rino, is located in the far-field plume region, 24 km off the coast (Figure 12 a). The narrow coastal current, however, is often situated closer to the coast. The plume passes through rino station only in cases where the winds erode the plume from the coast. Such events are seen in the observations on June 6 and 18. The model captures those events well, although the arrival of the plume is delayed by roughly one day. Salinity in Thetis also tends to be higher compared to the observations.

The southern station, riso, captures the wind-driven southward traversal of the plume. The observations show such an event between June 10 and 21. Thetis simulation shows the passage of the plume, albeit the onset is delayed by roughly 3 days, and surface salinity does not drop as much as in the observations.

Figure 13 shows the evolution of the surface salinity in the Thetis simulation. Prior to the first event (panels a and b), the bulge of the re-circulating plume is clearly visible just north of the river mouth. The far-field plume is quite narrow and fresh; it travels northward by the coast and does not reach the rino station. On the onset of northerly winds (panels c and d), the far-field plume detaches from the coast, spreads out, and passes through the station. As the winds prevail, the plume spreads out more and starts to travel southwest (panels e - h). Eventually, a new coastal plume starts to develop. Note that in response to the winds, the re-circulating plume also spreads out, dilutes and begins to meander.

6 Discussion

The idealized 1D and estuary test cases show that the GLS implementation remains stable and produces realistic bottom and surface layer dynamics as well as estuarine circulation. Overall the model’s performance is in good agreement with results presented in the literature (e.g., Warner et al. 2005; Kärnä et al. 2012).

All the tested closure models and stability functions yield similar results, although some differences do appear. For example, the differences in the bottom friction test case between different closures (Figure 2) are mostly due to numerics related to the surface boundary conditions. The commonly-used choice of $k-\varepsilon$ and Canuto A stability functions, however, appears to perform well in all test cases. The model performance is also strongly dependent on the chosen parameters which can be tuned in applications. For example, vertical profiles in the bottom friction test are affected by the choice of $c_{\mu}^{0}$ and $\sigma_{\Psi}$ parameters, among others. Similarly, the length of salinity intrusion in the idealized estuary test case (Figure 8) depends on the magnitude of total mixing which is partially controlled by the minimum value of TKE, $k_{min}$. Evidently, the parameters could be tuned further to reduce the gap between different closure models, for example. However, such an investigation of the parameter space is out of the scope of the present study.

The Columbia River plume simulation indicates that Thetis can represent tidal, highly dynamic river plumes. The plume features numerous eddies and responds quickly to changing tidal and wind conditions. Overall the simulation is in good agreement with the observations and known behavior of the Columbia plume under up- and downwelling conditions (e.g., Hickey and Banas, 2003). Qualitatively, Thetis results compare well with ROMS and SLIM (Vallaeys et al., 2018). Thetis surface fields show eddies and filaments, and strong frontal features are retained over several tidal cycles. This suggests that numerical mixing is low, as Kärnä et al. (2018) demonstrated with idealised test cases.

The presented model uses $\text{P}^{\text{DG}}_{0}$ function space for the turbulent quantities while the mean flow values are in $\text{P}^{\text{DG}}_{1}$ space. This choice resembles the staggered discretization used in finite volume models (e.g., in GOTM) where the mean flow and turbulent variables are defined at cell centers and cell interfaces, respectively. As the turbulent quantities reside in a lower-degree function space, it could impact the accuracy of the coupled model. However, as the turbulence closure model only affects the mean flow model through the eddy viscosity and diffusivity fields, the impact is likely to be small. Deriving a higher-order discretization for the turbulent quantities remains a topic for future research. Similarly, ensuring strict energy conservation between the turbulent and mean flow fields will be addressed in the future.

7 Conclusions

We have presented a DG finite element discretization for the GLS turbulence closure model and a positivity-preserving coupled time integration scheme. The turbulence closure model was implemented in the Thetis circulation model.

The model was tested with a series of idealized test cases. The results verify that bottom boundary layer, surface mixing, and estuarine circulation processes are simulated correctly. All the three closures ($k-\varepsilon$, $k-\omega$, and GLS) and stability functions (Canuto A, B, and Cheng) yield similar performance.

Finally, the model was tested in a realistic Columbia River plume application. The results indicate that the coupled model performs well in complex, highly dynamic plume applications. The surface salinity values are in good agreement with observations and the plume responds to changing wind forcing. The plume shows strong frontal features which suggests low numerical diffusion. Thetis, therefore, can retain fronts between water masses and capture the observed plume traversals.

There is an ever-increasing need for accurate, less dissipative coastal ocean models that can simulate buoyancy-driven coastal flows with a realistic representation of frontal features and vertical mixing. As such, the presented coupled circulation-turbulence model constitutes an important step towards next-generation ocean modeling.

Acknowledgments

I am grateful to anonymous reviewers whose insightful comments helped to improve the manuscript. Lawrence Mitchell (Durham University), Stephan Kramer and David Ham (Imperial College London) are thanked for their contributions to Thetis and continuous support with Firedrake. Antonio Baptista (Oregon Health & Science University) is acknowledged for his contribution on the Columbia River plume application.

Appendix A Stability functions

All the non-equilibrium stability functions are defined by a set of canonical parameters (Table 1 in Umlauf and Burchard, 2005). The functions can, however, be expressed as rational functions of $\alpha_{M}$ and $\alpha_{N}$ with parameters $n_{i},d_{i},n_{bi}$ (Umlauf and Burchard, 2005),

[TABLE]

Canuto et al. (2001) introduced two sets of full non-equilibrium stability functions, usually referred to as Canuto A and B. The stability region of Canuto B functions is somewhat larger than that of version A, suggesting that it is more robust in practice (Umlauf and Burchard, 2005). Another stability function, based on a very similar formulation, is given by Cheng et al. (2002). All the three stability function families have been implemented using the formulation (73)-(74).

A.1 Limiting $\alpha_{M}$ and $\alpha_{N}$

The stability functions are only defined for a certain values of $\alpha_{N}$ and $\alpha_{M}$ due to numerical and physical reasons. First, the shear frequency $\alpha_{M}$ is positive by definition. To ensure its positivity $\alpha_{N}$ must be limited from below. The condition $\alpha_{M}=0$ corresponds to the equilibrium relation $B=\varepsilon$, for which it holds (eq. A.23 in Umlauf and Burchard, 2005)

[TABLE]

This equilibrium condition corresponds to unstably stratified, convective regime where $\alpha_{N}<0$ and $B>0$. Substituting the stability function $c_{\mu}^{\prime}$ into (75) results in

[TABLE]

Typical value of $\alpha_{N}^{\text{min}}$ is around -3.

An alternative constrain to $\alpha_{N}$ is to require monotonicity (eq. 47 in Umlauf and Burchard, 2005):

[TABLE]

which yields more strict limit. In this work (76) is used as it was sufficient in all the tested cases.

Cropping $\alpha_{N}$ at the minimum value may lead to numerical oscillations as the gradient of the stability function is large for negative $\alpha_{N}$. Therefore a smoother transition has been proposed (eq. (19) in Burchard et al., 1999):

[TABLE]

The smoothness of the transition is controlled by $\alpha_{c}$, defined in range $\alpha_{N}^{\text{min}}<\alpha_{c}<0$. Here a value $\alpha_{c}=-1.2$ is used.

In addition, the shear frequency, $\alpha_{M}$, must be limited above to ensure physically sound behavior and numerical stability (eq. 44 in Umlauf and Burchard, 2005; Burchard and Deleersnijder, 2001):

[TABLE]

Note that $\alpha_{M}^{\text{max}}$ depends on $\alpha_{N}$.

The limit (79) is not always sufficient to ensure stability of the numerical scheme; Deleersnijder et al. (2008) present a more generic stability constraint which guarantees that the vertical gradient of the turbulent quantities remains bounded. In this work, (79) is used as it has been sufficient for stability: We first set $\alpha_{N}$ to $\tilde{\alpha}_{N}$ according to (78), and then apply the $\alpha_{M}^{\text{max}}$ limit.

Appendix B Choosing empirical parameters

B.1 Computing $c_{3}^{-}$

The parameter $c_{3}^{-}$ controls the buoyancy production of $\Psi$ under stable stratification. Its value can be computed from (34):

[TABLE]

The value of $c_{3}^{-}$ thus depends on both the closure (parameters $c_{i}$) and the stability functions. To evaluate (80) one first needs to evaluate the stability functions under steady state. The steady state condition (29) can also be written as (eq. (49) in Umlauf and Burchard, 2005)

[TABLE]

from which one obtains a non-linear equation

[TABLE]

From the above equation $\alpha_{M}$ can be solved either analytically or numerically. In practice, the analytical approach has been more robust and it is used in this work; the values for the three closures and the Canuto A stability functions are shown in Table 1.

B.2 Parameter $c_{\mu}^{0}$

Parameter $c_{\mu}^{0}$ stands for the value of the stability function in unstratified equilibrium where $P=\varepsilon$.

From (81) and (33) it follows that in equilibrium the stability functions only depend on $R_{i}^{st}$: $c_{\mu}=c_{\mu}(R_{i}^{st}),c^{\prime}_{\mu}=c^{\prime}_{\mu}(R_{i}^{st})$. Its value can therefore be computed for each stability function (see eq. A.22 in Umlauf and Burchard, 2005).

B.3 Parameter $\sigma_{\Psi}$

A relation for the Schmidt number, $\sigma_{\Psi}$, and the von Karman constant $\kappa$ can be derived for the logarithmic boundary layer assumption (eq. 14 in Umlauf and Burchard, 2003):

[TABLE]

In this work we compute $\sigma_{\Psi}$ for each closure with $\kappa=0.4$.

Appendix C Source code availability

The source code used to perform the experiments in this papers is publicly available. Firedrake, and its components, may be obtained from https://www.firedrakeproject.org/ (last access: 1 February 2020); Thetis from http://thetisproject.org/ (last access: 1 February 2020). For reproducibility, we also cite archives of the exact software versions used to produce the results in this paper. All major Firedrake components and the Thetis source code have been archived on Zenodo (zenodo/Firedrake, 2019; zenodo/Thetis, 2019).