Numerics

For the spatial discretisation, the water column is divided into $N_i$ layers of not necessarily equal thickness $h_i$,

$$h_i=(\gamma_i-\gamma_{i-1})D, \qquad i=1,\dots,N_i \; , \quad$$ (3)

with nondimensional interfaces $\gamma_i$ with $\gamma_0=-1$, $\gamma_{i-1}< \gamma_i$ and $\gamma_{N_i}=0$, see Burchard and Petersen (1997).

The discrete values for the mean flow quantities $U$, $V$, $\Theta$, and $S$ represent interval means and are therefore located at the centres of the intervals, and the turbulent quantities like $k$, $L$, $\epsilon$, $\nu_t$, $\nu'_t$, $N$, $P$, $G$, $c_{\mu}$, and $c_{\mu}'$ are positioned at the interfaces of the intervals (see section 4.7). The indexing is such, that the interface above an interval has the same index as the interval itself. This means that mean flow quantities range from $i=1,..,N_i$ while turbulent quantities range from $i=0,..,N_i$ (see figure 1).

Figure 1: Spatial organisation and indexing of the numerical grid.

The staggering of the grid allows for a straight-forward discretisation of the vertical fluxes of momentum and tracers without averaging. However, for the vertical fluxes of e.g. $k$ and $\epsilon$, averaging of the eddy diffusivities is necessary. This is only problematic for the fluxes near the surface and the bottom, where viscosities at the boundaries have to be considered for the averaging. These can however be derived from the law of the wall.

Figure 2: Temporal organisation and indexing of the numerical grid. Here, a time stepping slightly more implicit than the Crank and Nicolson (1947) scheme with $\sigma =0.6$ is shown.

The time stepping is equidistant, based on two time levels and not limited by Courant numbers, because of the absence of advection and an implicit treatment of vertical diffusion, see figure 2. In the following, the discretisation of a simple diffusion equation,

$$\dfrac{\partial {X}}{\partial {t}} - \dfrac{\partial {}}{\partial {z}} \left( \nu \dfrac{\partial {X}}{\partial {z}} \right) = 0 \; , \quad$$ (4)

will be illustrated for Neumann-type boundary conditions

$$\nu \dfrac{\partial {X}}{\partial {z}} = F_s z=\zeta,\qquad$$ (5)

and

$$\nu \dfrac{\partial {X}}{\partial {z}} = F_b z=-H.\qquad$$ (6)

The semi-implicit discretisation of (4) can then be written as

$$\displaystyle \frac{X^{n+1}_{N_i}-X^n_{N_i}}{\Delta t} -\frac{F_s... ...gma}_{N_i-1}}{0.5(h^{n+1}_{N_i}+h^{n+1}_{N_i-1})}}{h^{n +1}_{N_i}} = \; , \quad$$ (7)

$$\displaystyle \frac{X^{n+1}_i-X^n_i}{\Delta t} -\frac{\nu^n_i\fra... ...i}-X^{n+\sigma}_{i-1}}{0.5(h^{n+1}_i+h^{n+1}_{i-1})}}{h^{n +1}_i} =0 \; , \quad$$ (8)

$$\displaystyle \frac{X^{n+1}_1-X^n_1}{\Delta t} -\frac{\nu^n_1\fra... ...2}-X^{n+\sigma}_{1}}{0.5(h^{n+1}_{2}+h^{n+1}_1)} -F_b}{h^{n+1}_1} =0 \; , \quad$$ (9)

for $1<i<N_i$. Here, the semi-implicit time level is defined by

$$X^{n+\sigma}=\sigma X^{n+1}+(1-\sigma)X^n.$$ (10)

Thus, for $\sigma=0$, a fully explicit, for $\sigma=1$ a fully implicit, and for $\sigma=0.5$ the Crank and Nicolson (1947) second-order scheme are obtained. Figure 2 shows an example for $\sigma =0.6$. It should be noted that often a time stepping is preferable which is slightly more implicit than the Crank and Nicolson (1947) scheme in order to obtain asymptotic stability. The resulting linear system of equations (7) - (9) with tri-diagonal matrix structure is solved by means of the simplified Gaussian elimination.

With the same strategy, a very similar system of equations can be derived for variables located at the interfaces of the grid cells, i.e. variables describing turbulence.