Generate a two-equation model

INTERFACE:

    subroutine generate_model

DESCRIPTION:

Computes the parameters of an instance of the `generic' two-equation model according to the specifications set in gotmturb.nml. This model solves (150) for the $k$ and (166) for the generic length-scale defined in section 4.7.29 together with an Algebraic Stress Model. For several simple turbulent flows, analytical solutions of this models exist and can be used to calibrate the model coefficients. The method is described in great detail in Umlauf and Burchard (2003). Also users that are not interested in the generic part of GOTM should have a look in this section, because results derived here are referenced in later parts of the manual.

After the call to generate_model(), all parameters of the generic two-equation model are known. The user has full control over specific properties of the resulting model (see section 4.7.29).

In the following sections, the effects of model parameters on the behaviour of two-equation models in specific situations are briefly reviewed. For a more in-depth discussion, see Umlauf and Burchard (2003).

The logarithmic boundary layer

In the logarithmic boundary layer one has $P=\epsilon$ and $k\propto u^2_*$ by defintion. Under these conditions it is easy to show that a solution of (150) is

$$k = \dfrac{u_*^2}{(c_\mu^0)^2} \; , \quad$$ (101)

and a solution of (166) can only be obtained if the condition

$$\sigma_\psi = \dfrac{n^2 \kappa^2}{(c_\mu^0)^2 (c_{\psi 2}-c_{\psi 1})}$$ (102)

is satisfied. (101) can be conveniently used to obtain boundary conditions for $k$, whereas (102) yields for example the value for the turbulent Schmidt-number $\sigma_\psi$ as a function of the von Kármán constant (provided, of course, that the other constants are known). The value of the von Kármán constant is usually assumed to be $\kappa \approx 0.4$.

Decay of homogeneous turbulence

Another example of a simple but fundamental turbulent situation is the temporal decay of isotropic, homogeneous turbulence (approximated by the spatial decay of turbulence behind grids in laboratory settings). At large times, $t$, data from many experiments are well described by a power law of the form

$$\frac{k}{k_0} = A \left( \frac{t}{\tau_0} \right)^d \; , \quad$$ (103)

with constant $A$ and initial values of the kinetic energy, $k_0$, and the eddy turnover time, $\tau_0$. The decay rates, $d$, have been thoroughly documented. Experiments (Bradshaw (1975), Townsend (1976), Domaradzki and Mellor (1984), Mohamed and Larue (1990)) suggest that $d$ is in the range $-1.3 < d < -1$. DNS, generally conducted at low Reynolds numbers, produce consistently higher values. For example, Briggs et al. (1996) obtain a value near $-1.5$ from their DNS.

In homogeneous decaying turbulence, (150) and (166) reduce to a balance between the rate and dissipation terms, respectively. The coupled system of ordinary differential equations can be solved for given initial values $k_0$ and $\psi_0$ (see e.g. Wilcox (1998)). The solution can be shown to reduce to (103) at large times. Then, the decay exponent, $d$, is determined by

$$d = - \frac{2n}{2m+n-2c_{\psi 2}} \; , \quad$$ (104)

and thus depends only on the exponents $m$ and $n$ defined in (165) and the model constant $c_{\psi 2}$. For given exponents $m$ and $n$, the experimental values of $d$ can be used to derive the value of the model constant $c_{\psi 2}$. Note, that the predicted decay rate, $d$, is completely independent of the ASM (or the stability functions in other words).

Homogeneous turbulent shear-flows

A natural extension of decaying homogeneous turbulence is the inclusion of a homogeneous shear and an aligned homogeneous stratification. Since turbulence is still assumed to be homogeneous, the divergence of any turbulent transport term vanishes and the interplay between the stabilizing effects of stratification and the destabilizing action of shear can be isolated. Thus, it is not surprising that this interesting special case of turbulence has been explored extensively by laboratory experiments (Tavoularis and Corrsin (1981a); Tavoularis and Corrsin (1981b), Tavoularis and Karnik (1989), Rohr et al. (1988)), by Direct Numerical Simulation (Gerz et al. (1989), Holt et al. (1991), Jacobitz et al. (1997), Shih et al. (2000)) and by Large-Eddy Simulation (Kaltenbach et al. (1994)). That flows of this kind are crucial also in many oceanographic flows has been pointed out by Baumert and Peters (2000).

In the context of the generic two-equation model, this turbulent flow is mathematically established by neglecting the turbulent transport terms and the advective part of the material time derivative. Then, (150) and (166) reduce to a set of ordinary differential equations. Using the chain rule of differentiation, the relation

$$\frac{1}{l} \dfrac{\text{d} {l}}{\text{d} {t}} = \frac{1}{n} \fra... ...si}}{\text{d} {t}} - \frac{m}{n} \frac{1}{k} \dfrac{\text{d} {k}}{\text{d} {t}}$$ (105)

for the mixing length, $l$, follows immediately from (165). With (105), the generic model expressed by (150) and (166) can be used to derive an evolution equation for the integral length scale, $l$,

$$\begin{array}{rcl} \dfrac{1}{l} \dfrac{\text{d} {l}}{\text{d}... ...c_{\psi_3} - \dfrac{m}{n} \right) G \right) \quad . \end{array}$$ (106)

Tennekes (1989) derived an equation similar to (106), however only for the special case of the $k$-$\epsilon$ model applied to unstratified flows, and stated that `on dimensional grounds, $l$ cannot depend upon the shear because the shear is homogeneous and cannot impose a length scale'. This argument requires immediately

$$c_{\psi 1} = m \; , \quad$$ (107)

which is used in the subroutine to determine the model parameter $c_{\psi 1}$. A more detailed discussion of this method is given in Umlauf and Burchard (2003).

Shear-free turbulence, wave-breaking

The first step in understanding the behaviour of two-equation models in the surface layer affected by breaking gravity waves is the investigation of a special case, in which turbulence decays spatially away from a planar source without mean shear. Turbulence generated by an oscillating grid in a water tank has been used in various laboratory settings to study the spatial decay of velocity fluctuations in this basic turbulent flow, where turbulent transport and dissipation balance exactly. For a summary of these results, see Umlauf et al. (2003).

All grid stirring experiments confirmed a power law for the decay of $k$ and a linear increase of the length scale, $l$, according to

$$k = K (z+z_0)^\alpha \; , \quad l = L(z+z_0) \; , \quad$$ (108)

where $K$, $L$, and $z_0$ are constants, and the source of turbulence has been assumed to be at $z=0$. In these experiments, $z_0 = l / L$ at $z=0$ is not related to any kind of surface roughness length. Rather, it is connected to the length scale of injected turbulence, uniquely determined by the spectral properties of turbulence at the source. Experiments suggest that the decay rate for the turbulent kinetic energy is likely to be in the range $-3<\alpha<-2$. The value of $L$, i.e. the slope of the turbulent length scale, $l$, was found to be consistently smaller than in wall-bounded shear flows, $L < \kappa \approx 0.4$, see Umlauf et al. (2003).

In stationary, shear-free, unstratified turbulence, the generic model simplifies to a balance between the turbulent transport terms and the dissipative terms in (150) and (166). Using the definition of $\psi$ in (165) and the scaling for the rate of dissipation, (153), the transport and dissipation of $k$ and $\psi$ are balanced according to

$$\begin{array}{rcl} \dfrac{\text{d} }{\text{d} {z}} \left( \df... ...} (c_\mu^0)^{p+3} k^{m+\frac{1}{2}} l^{n-1} \quad . \end{array}$$ (109)

Note, that in shear-free turbulence, the shear-number defined in (45) is $\alpha_M=0$ by definition, and stability functions always reduce to a constant which is, however, different from the constant $c_\mu^0$ approached in the logarithmic boundary layer, see section 4.7.13.

For the solution of this non-linear system , we inserted the expressions (108) in (109). From (153) and (44), power-laws follow then also for $\epsilon = E(z+z_0)^\beta$ and $\nu_t = N (z+z_0)^\gamma$.

Inserting (108) into (109)$_1$ yields the equation

$$(\alpha L)^2 = \frac{2}{3} (c_\mu^0)^2 R \, \sigma_k^\psi \; , \quad$$ (110)

where the constant ratio $R=c_\mu^0/c_\mu$ follows uniquely from the respective ASM. The power-law (108) can also be inserted in (109)$_2$ to yield

$$\left( \alpha m + n \right) \left( \left( \dfrac{1}{2} + m \right... ...a + n \right) L^2 = \left(c_\mu^0 \right)^2 R \, \sigma_\psi c_{\psi 2} \quad .$$ (111)

We note that with the help of (104) and (107), the relation (102) can be rewritten as

$$\sigma_\psi = \frac{2 \kappa^2 d}{(c_\mu^0)^2 (d+2)} n \quad .$$ (112)

Expressing now $\sigma_\psi$ with (112) and $c_{\psi 2}$ with the help of (104) on the right hand side of (111), an equation expressing the exponent $m$ in terms of $n$ (or vice-versa) can be obtained. The result for $n$ can be written as

$$\begin{array}{rcl} n &=& - \dfrac{1}{4 (2+d) (\kappa^2 R - L^... ...+ (2+d) (1+4m) \alpha L^2 \Big)^2 } \right) \quad . \end{array}$$ (113)

After assigning appropriate values for the von Kármán constant, $\kappa$, the decay coefficient of homogeneous turbulence, $d$, the spatial decay rate, $\alpha$, and the slope, $L$, an infinite number of pairs of $m$ and $n$ satisfying (113) can be derived. Each corresponds to a different two-equation model. Some example are given in table 5 (see Umlauf and Burchard (2003)).

Table 5: Some parameter sets for the generic model with $\kappa = 0.4$, $d = -1.2$, $(c_\mu ^0)^2=0.3$, $c_{\psi _1}=m$ and obeying the log-layer compatibility relation (112).

$\alpha$	$L$	$m$	$n$	$c_{\psi 2}$	$\sigma_k^\psi$	$\sigma_\psi$
$-2.0$	$0.20$	$1.00$	$-0.67$	$1.22$	$0.80$	$1.07$
$-2.0$	$0.20$	$2.00$	$-1.09$	$2.36$	$0.80$	$1.75$
$-2.5$	$0.20$	$1.00$	$-1.05$	$1.35$	$1.25$	$1.68$
$-2.5$	$0.20$	$2.00$	$-1.74$	$2.58$	$1.25$	$2.78$

Even though each line in this table represents a different two-equation model with completely different model constants, each of the two groups of models (with $\alpha=-2.0$ and $\alpha=-2.5$, respectively) performs completely identical in all situations discussed until here. Thus, the generic model allows for the formulation of groups of two-equation models with fully controlled properties from the outset. As discussed by Umlauf and Burchard (2003), one more constraint is necessary to obtain the final values of all parameters, including the exponents $m$ and $n$. These authors suggested that the first line in table 5 yields a model with excellent properties in all flows they considered.

Mixed layer deepending

The correct prediction of mixed layer deepening into a stratified fluid due to a wind stress at the surface is one of the most crucial requirements for an oceanic turbulence model. This situation has been frequently interpreted by analogy with the classical experiment of Kato and Phillips (1969) and its re-interpretation by Price (1979), in which the entrainment in a linearly stratified fluid subject to a constant surface stress was investigated. The results of this experiment have been used by numerous authors to calibrate their turbulence models.

In particular, it has been shown by Burchard and Bolding (2001) for the $k$-$\epsilon$ model of Rodi (1987), by Burchard (2001b) for the $q^2 l$ model of Mellor and Yamada (1982), and by Umlauf et al. (2003) for the $k$-$\omega$ model of Wilcox (1988) that, remarkably, the mixed layer depth predicted by these models depends almost exclusively on the value of the Richardson number, $Ri=N^2/M^2$, computed in a homogeneous, stratified shear-flow in steady-state. This value is usually referred to as the steady-state Richardson number, $Ri_{st}$ (Rohr et al. (1988), Kaltenbach et al. (1994), Jacobitz et al. (1997), Shih et al. (2000)).

Umlauf et al. (2003) showed that in the context of models considered in GOTM, the steady-state Richardson number is determined by the relation

$$Ri_{st}=\dfrac{c_\mu}{{c_\mu}'} \dfrac{c_{\psi 2} - c_{\psi 1}}{c_{\psi 2} - c_{\psi 3}} \quad .$$ (114)

Since it is well-known that, with the equilibrium assumption $P+G=\epsilon$, stability functions reduce to functions of $Ri$ only (Mellor and Yamada (1974), Galperin et al. (1988)), (114) is a non-linear equation for the model constant $c_{\psi 3}$ for given $Ri_{st}$. Note, that the structure parameters, $m$ and $n$, do not appear in (114). This implies that the type of the two-equation model is irrelevant for the prediction of the mixed layer depth, as long as (114) is fulfilled for identical $Ri_{st}$. Numerical examples with very different values of $m$ and $n$ confirmed indeed that the mixed layer depth only depends on $Ri_{st}$. The experiment of Kato and Phillips (1969) could almost perfectly be reproduced, provided the parameter $c_{\psi 3}$ was chosen to correspond to $Ri_{st}\approx0.25$, see Umlauf et al. (2003).

Note, that in instable situations, a different value of the parameter $c_{\psi 3}$ needs to be used. This does not cause a discontinuity in the model because the buoyancy term in (166) is zero at the transition. An evaluation of the length-scale equations in convective flows, however, is intimately related to the third-order modelling of the triple correlation terms, a topic outside the scope of this documentation.

USES:

   IMPLICIT NONE

REVISION HISTORY:

    Original author(s): Lars Umlauf