In search of improving the numerical accuracy of the k−ɛ model by a transformation to the k−τ model

Ocean Modelling 104 (2016) 129–142

Contents lists available at ScienceDirect

Ocean Modelling journal homepage: www.elsevier.com/locate/ocemod

Virtual Special Issue Coastal ocean modelling

In search of improving the numerical accuracy of the k − ε model by a transformation to the k − τ model Yoeri M. Dijkstra a,b,∗, Rob E. Uittenbogaard a, Jan A.Th.M. van Kester a, Julie D. Pietrzak b a b

Deltares, Boussinesqweg 1, NL-2629 HV Delft, The Netherlands Section of Environmental Fluid Mechanics, Delft University of Technology, P.O. Box, 5048, NL-2600 GA, Delft, The Netherlands

a r t i c l e

i n f o

Article history: Received 21 July 2015 Revised 14 May 2016 Accepted 31 May 2016 Available online 31 May 2016 Keywords: Turbulence modelling Numerical accuracy k − ε model k − τ model Boundary layers

a b s t r a c t This study presents a detailed comparison between the k − ε and k − τ turbulence models. It is demonstrated that the numerical accuracy of the k − ε turbulence model can be improved in geophysical and environmental high Reynolds number boundary layer flows. This is achieved by transforming the k − ε model to the k − τ model, so that both models use the same physical parametrisation. The models therefore only differ in numerical aspects. A comparison between the two models is carried out using four idealised one-dimensional vertical (1DV) test cases. The advantage of a 1DV model is that it is feasible to carry out convergence tests with grids containing 5 to several thousands of vertical layers. It is shown hat the k − τ model is more accurate than the k − ε model in stratified and non-stratified boundary layer flows for grid resolutions between 10 and 100 layers. The k − τ model also shows a more monotonous convergence behaviour than the k − ε model. The price for the improved accuracy is about 20% more computational time for the k − τ model, which is due to additional terms in the model equations. The improved performance of the k − τ model is explained by the linearity of τ in the boundary layer and the better defined boundary condition. © 2016 Published by Elsevier Ltd.

1. Introduction Resolving the physical processes that occur in the boundary layer is of great importance to the modelling of flows in the coastal ocean, estuaries and rivers. Many of these physical processes are associated with turbulent mixing and need to be parametrised in large-scale three-dimensional models. Mixing processes in these models are generally parametrised by the use of two-equation turbulence models, popular examples of which are the k − ε model (Jones and Launder, 1972; Rodi, 1993), the buoyancy extended version of the k − ω model (Wilcox, 1988; Umlauf et al., 2003) and the k − kl model (Mellor and Yamada, 1982), but many more exist. All of these models parametrise the physics of turbulence using similar, but not identical methods. In addition, these models differ in terms of the numerical choices made, which in turn affect the numerical accuracy. It is unclear which differences in model accuracy are due to differences in the physical parametrisations and which are due to the numerical accuracy. Moreover, such a distinc∗

Corresponding author. E-mail address: [email protected], [email protected] (Y.M. Dijkstra).

http://dx.doi.org/10.1016/j.ocemod.2016.05.010 1463-5003/© 2016 Published by Elsevier Ltd.

tion is not easily made. The question that then arises is whether those turbulence models that are most popular for their physical accuracy also possess the most favourable numerical properties. More specifically, the question arises whether these models have the highest possible accuracy at resolutions typically used in 3D modelling. This question is addressed here by comparing two different turbulence models at a wide range of different grid resolutions up to a resolution where the numerical error almost vanishes. In this paper, a detailed comparison of the k − ε and k − τ models is carried out. Here the effect of the grid resolution on the numerical accuracy of these models is explored. Significantly, we show here that it is possible to improve upon the numerical properties of the k − ε model in frictional boundary layers by transforming it to the k − τ model. The k − τ model is obtained from the k − ε model by using the definition τ = k/ε . The new variable τ represents a typical timescale of turbulent eddies. The k − τ model has been derived earlier by Speziale et al. (1992) and was tested by Thangam et al. (1992). These tests focussed on non-stratified boundary layers in industrial and aerospace applications. The model has otherwise been largely disregarded. According to Hanjalic´ (1994), one reason for

130

Y.M. Dijkstra et al. / Ocean Modelling 104 (2016) 129–142

Fig. 1. Schematised representation of the vertical profiles of ε and τ in a steady, uniform flow that is dominated by bed friction. The profile of ε is hyperbolic and its value goes to infinity at the bed, while the profile of τ is linear near the bed and finite at the bed itself.

this is that it does not bring new insights into the physics of turbulence in comparison with the k − ε model. However, approaching the model from the perspective of numerical accuracy, it will be shown here that the k − τ model has several superior properties over the k − ε model in boundary layer flows. The k − ε model has two particular properties in boundary layers that are disadvantages from a numerical perspective. First, the profile of ε is hyperbolic, see Fig. 1. The gradients of such a hyperbolic profile are hard to approximate on a coarse grid. Second, the numerical approximation of the boundary value of ε is severely inaccurate (Burchard and Petersen, 1999). In contrast, the profile of τ is linear in the boundary layer, so that the profile is exactly represented in any first-order or higher-order numerical scheme. Additionally, the boundary condition is much better behaved. It is generally difficult to conclude that one turbulence model is better than another. This is because the parametrisation of turbulence is highly complex; it involves a large number of assumptions, which differ slightly between different models. This makes it difficult to say which set of assumptions is the most accurate for a general field of application. As a result, a large number of papers exist comparing different two-equation turbulence models for different cases. An excellent review is provided by Warner et al. (2005), concluding that there is no clear theoretical preference for any particular two-equation model. Moreover, numerous comparisons of these models for a wide range of environmental and geophysical cases have not resulted in the identification of a single best model (Martin, 1985; Burchard et al., 1998; Umlauf and Burchard, 2003; Wijesekera et al., 2003; Warner et al., 2005; Ilicak et al., 2008). The interpretation of these comparisons is further complicated by the fact that the accuracy of turbulence modelling is not only affected by the model equations, but also by the numerical implementation (Patankar, 1980; Stelling, 1995; Ezer et al., 2002; Ezer and Mellor, 2004) and the way in which other model components and the computational grid are designed (Warner et al., 2005; Platzek et al., 2014). In order to gain a better understanding of the nature of the differences between the turbulence models, it is useful to make a clear distinction between the physical and numerical aspects. Restricting our attention to high Reynolds number flows, most twoequation models have a similar equation structure. However, the models differ in the way the coefficients are calibrated and the way in which the diffusion of turbulence is parametrised. Therefore, the models represent slightly different physical parametrisations. Most model comparisons focus on these different parametrisations. In so doing, the importance of numerical differences between the models remains unclear. Here we isolate the numerical differences between the turbulence models by using a transformation from the k − ε model to the k − τ model. The two models thus implement the same physical parametrisation. Consequently, differences in the results can be traced back to the numerical implementation.

A comparison between the k − ε and k − τ models is carried out using four test cases. These test cases represent types of flows commonly encountered in environmental and geophysical modelling. The first case is an open channel flow case, which has been tested by Burchard et al. (1998) and Warner et al. (2005). This case represents the class of non-stratified boundary layer flows and is used to quantify the effects of differences between the profiles of ε and τ , as shown in Fig. 1. The second test case is an idealised estuarine flow, such as tested by Nunez Vaz and Simpson (1994), Monismith and Fong (1996) and Burchard and Hetland (2010). This case represents stratified boundary layer flows. It indicates to what extent the advantages of the k − τ model in the boundary layer remain relevant if the flow involves stratification. The final two cases focus on flows dominated by buoyancy processes: the mixing of unstable stratification (Burchard and Petersen, 1999) and the simulation of the temperature in an idealised lake model. Boundary friction is not important in the latter two test cases. Mixing in three-dimensional models is dominated by processes acting in the vertical dimension. A one-dimensional vertical (1DV) model is therefore used to isolate these processes. The lower computational costs of the 1DV model compared to 3D models makes it possible to carry out test cases on high resolution grids, of up to several thousands of vertical layers. At these resolutions the numerical error should become negligibly small. The numerical accuracy of the models is tested at a number of resolutions, ranging from 5 to 10,0 0 0 layers in the first case and 5 to 20 0 0 layers in the other cases. In this way the convergence behaviour of the turbulence model with increasing resolution is explored. The 1DV model closely corresponds, apart from some details, to the vertical grid design and k − ε model implementation of Delft 3D-FLOW (Uittenbogaard et al., 1992; Deltares, 2014), which is a full-featured 3D flow and transport model, specialised for coastal, estuarine and river applications. The k − τ model is currently not included in Delft 3D-FLOW. A new k − τ routine for Delft 3D-FLOW has therefore been developed in this study. We will provide the details of the numerical implementation of both the Delft 3D-FLOW version of the k − ε model and the new implementation of the k − τ model in Section 2. We will then present the results of the four case studies in Section 3. A discussion of the results is presented in Section 4.

2. Model 2.1. Model equations The 1DV model solves the equations for mass and momentum conservation for flow in one direction. The momentum equation uses an eddy viscosity to parametrise the effects of turbulent mixing. These equations are simplified further by using the hydrostatic pressure assumption and by neglecting the Coriolis terms. The resulting momentum equations can be formulated as

  1 ∂p ∂u ∂ ∂u =− + , (ν + νt ) ∂t ρ0 ∂ x ∂ z ∂z ∂p = −ρ g. ∂z

(1)

(2)

In these equations the horizontal velocity is denoted by u in the Cartesian coordinates (x, z), with the z-axis pointing upwards. Pressure is denoted by p, ρ is the density, ρ 0 is a reference density, ν is the molecular viscosity of water and ν t is the eddy viscosity. Boundary conditions to these equations are prescribed at the bottom and free surface. These boundary conditions are given by

Y.M. Dijkstra et al. / Ocean Modelling 104 (2016) 129–142

The k − ε model in 1DV reads

Table 1 Parameter values used in the k − ε and k − τ models.

cε 1 , cε 2 , cε 3 , cε 3 , cμ

cτ 1 cτ 2 cτ 3 , stable stratification cτ 3 , unstable stratification

σρ σk σ ε, σ τ

k−ε

k−τ

1.44 1.92 0 1.44 0.09 0.7 1.0 1.3

−0.44 −0.92 1 −0.44 0.09 0.7 1.0 1.3

131

  ∂k ∂  ν  ∂k νt 2 ν+ t = + νt M2 − N − 

ε ,   σρ ∂t ∂ z σk ∂ z  

Ek 

 Pk

(8)

−Bk

Dk

  ∂ε ∂  ν  ∂ε ε ε νt 2 ε2 = + cε1 νt M2 − cε3 N − cε 2 . ν+ t ∂t ∂z σε ∂ z k k σ

  ρ   k

 

  P E ε

Dε

−Bε

ε

(9)

quadratic-slip formulations, which for the u-equation read

 ∂ u  νt  ∂z

= u∗,b |u∗,b |,

(3)

z=−H

νt

 ∂ u  = u∗,s |u∗,s |, ∂ z z=ζ

(4)

where u∗, b and u∗, s are the friction velocities at the bed and surface respectively. Assuming a rough boundary in a turbulent flow, these friction velocities are derived from the assumption that the velocity profile is logarithmic. Numerically, this means that the velocity point closest to the bed is used to compute u∗, b . As this velocity point is half a grid cell ( 12 z) from the boundary it follows



that u∗,b = u−H+ 1 z κ / ln 1 + 2

1 z 2 z0



, where κ is the Von Kármán

coefficient and z0 is the roughness height. The bed level is denoted by z = −H. The surface friction velocity follows from wind or wave impact, which is not taken into account in this study, i.e. u∗,s = 0. The surface level in general is denoted by z = ζ . As our experiments all use a rigid lid, we set ζ = 0. The 1DV model additionally solves an equation to model the vertical distribution of either temperature or salinity, depending on the application. Both are modelled by the advection–diffusion equation

  ∂c ∂c ∂ ∂c +u = D + D , ( t) ∂t ∂x ∂z ∂z

k2

ε

, or

νt = cμ kτ .

 ∂ u 2

 2 2 g ∂ρ 2 + ∂v ∂z ∂ z , N is the buoyancy frequency N = − ρ0 ∂ z and cε1 , cε2 and cε3 are coefficients, see Table 1. The symbols under the terms in the equations denote a short-hand notation, where D denotes the diffusion of turbulence, P is the shear production, B is the buoyancy production or dissipation and E is the viscous dissipation. The k − ε model can be transformed to the k − τ model by using the definition τ = εk . This yields the following expression for the k − τ model

  k ∂k ∂  ν  ∂k νt 2 = + νt M2 − N − , ν+ t   σρ ∂t ∂ z σk ∂ z τ 

 

 

Pk E

(6)

(7)

The coefficient cμ in this equation is treated as a constant, see Table 1.

−Bk

Dk

(10)

k

  ∂τ ∂  νt  ∂τ τ τ νt 2 = + cτ 1 νt M2 − cτ 3 N − cτ 2 ν+ ∂t ∂z στ ∂ z k k σ

  ρ 

 

  Eτ P τ

Dτ

−Bτ

 ν  ∂τ ∂ k 1 ν  ∂τ ∂τ ν+ t ν+ t −2 στ ∂ z ∂ z τ σ ∂z ∂z  

  τ

2 + k

(5)

where c can denote either temperature or salinity and D and Dt are the molecular and turbulent diffusivity. The term u ∂∂ xc on the lefthand side denotes advection of c for a prescribed and vertically and temporally constant horizontal gradient of c. The advection– diffusion equation has no-flux boundary conditions at the surface and at the bed. The calculated temperature and salinity are converted to a density using a linear equation of state in the case of estuarine flow and unstable stratification (case 2 and 3) and using the UNESCO 1981 equation of state in the lake case (case 4). The resulting density profile affects the flow in two ways: via the turbulence model and via the baroclinic pressure in the momentum equations. The final equations that are solved by the 1DV model are those for the two-equation turbulence model. In this study, either the k − ε model or the k − τ model is used as the turbulence model. The eddy viscosity follows from the turbulence parameters k, ε and τ according to

νt = cμ

The σ k , σ ε and σ ρ are constant Prandtl–Schmidt numbers, M2 =



Dkτ

−

τ ∂ k ∂z 



 ∂k  1 1 ν − . σε σk t ∂ z 

Dτ τ

(11)

Dkk

The structure of the first four terms on the right-hand side of the τ -equation is the same as in the ε-equation. Eq. (11) contains three additional terms which contain gradients of k or τ . These terms are called diffusive-type terms as they are derived from diffusive term in the k − ε model. These terms are denoted in short by Dkτ , Dτ τ and Dkk . Additionally, the values of the coefficients differ from the k − ε model, see Table 1. These additional diffusive terms and different coefficients are a consequence of the transformation. In general, any transformation of the k − ε model to other variables will yield a model with a similar structure with the possible addition of several diffusive-type terms and different values of coefficients. Similarly one could derive a k − τ model from, for example, the k − ω model or the k − kl model. Such versions of the k − τ model would also differ from Eq. (11) in their diffusive-type terms and the value of the coefficients. The values of the coefficients used in this research are listed in Table 1 and correspond to the formulation of Launder and Spalding (1974), Rodi (1993) and Delft 3D-FLOW (Uittenbogaard et al., 1992; Deltares, 2014). The coefficients cτ 1 , cτ 2 , cτ 3 for the k − τ model follow from the transformation and satisfy the relation cτ i = 1 − cεi for i = 1, 2, 3.

132

Y.M. Dijkstra et al. / Ocean Modelling 104 (2016) 129–142

diffusion of turbulent kinetic energy k at grid interface i. Its numerical implementation uses a fully implicit method.

Dk,i =

  n+1 νn k − kni +1 ν + t,i+1/2 i+1 σk

zi+1   n+1  n n+1 νt,i k − k −1/2 i i−1 − ν+ . σk

z i 1

zi+1/2

(16)

The eddy viscosity νt,i+1/2 in this expression is defined as νt,i+1/2 =  1 2 νt,i+1 + νt,i . The other terms in the k-equation are implemented according to n Pk,i = νt,i Mi2 ,

Fig. 2. Numbering of the staggered 1DV grid.

The k − ε model uses a Dirichlet boundary condition for k and a Neumann boundary condition for ε , following the recommendation of Burchard and Petersen (1999). These boundary conditions read

Bk,i =

⎧ n ν kn+1 ⎪ ⎨− σt,iρ Ni2 ki ni if Ni2 > 0; stable stratification ⎪ ⎩− νt,in N2 σρ

k|z=−H =

u2∗b

1/2 cμ

,

 ∂ε  |u |3 = − ∗b2 .  ∂ z z=−H κ z0

(12)

Ek,i = (13)

Similar boundary conditions are applied at the surface. The profile of τ does not go to infinity near a frictional boundary (see Fig. 1). We therefore formulate a Dirichlet boundary condition, by transforming the Dirichlet condition for ε to τ . The boundary condition for τ reads

τ |z=−H =

κ z0 . cμ |u∗b | 1/2

(14)

(17)

(18)

if Ni2 < 0; unstable stratification

i

⎧ kn+1 ⎪ ⎨2εin ki ni − εin , k − ε model (19)

⎪ ⎩ kni +1n ,

k − τ model.

τi

The terms Bk, i for stable stratification and Ek, i can act as a sink to the equations and could potentially cause k to become negative. However, the turbulence equations guarantee that k, ε and τ remain positive at all times (Mohammadi and Pironneau, 1994). It is important to make sure that this positivity is also guaranteed in the numerical implementation. To this end, a method by Patankar (1980) is used. Applied to the k-equation, this method uses prodkni +1 kni

ucts of the type

in sink terms to guarantee that k remains

2.2. Discretisation

positive.

The numerical implementation of the k − ε model is close to the implementation of the vertical dimension in Delft 3D-FLOW. The only exception is that we will use a Neumann boundary condition at the bed, while Delft 3D-FLOW uses an adapted Dirichlet condition (Van Kester, 1994). The k − τ model is currently not implemented in Delft 3D-FLOW. We have therefore developed an implementation of the k − τ model based on similar numerical operators. The model uses a staggered grid, such as is drawn schematically in Fig. 2. The velocity and pressure are defined in the cell centres (dots in the figure), while the turbulence parameters k, ε , τ and ν t are defined at the cell interfaces (horizontal lines in the figure). The cell interfaces are numbered i = 0, . . . , I, with interface 0 equal to the bed and I at the surface. The cell centres are numbered i = 1/2, . . . , I − 1/2. The model solves the equations for k and ε or τ independently and does not include an iterative procedure to deal with the non-linear coupling of the equations. As long as the time step is small relative to the time-scale of the flow, the change of turbulence between two time levels is fairly small, therefore justifying the non-iterative procedure. The implementation of the three equations for k, ε and τ are provided below.

2.2.2. Discretisation of the ε -equation The equation for ε has a similar structure to that of k and can be written as

2.2.1. Discretisation of the k-equation The equation for k is integrated in time by one-step procedure. The resulting equation can be written in short as

kni +1 − kni

t

= Dk,i + Pk,i + Bk,i − Ek,i .

(15)

The terms on the right-hand side of the equation correspond to the short-hand notation used in Eq. (8). The term Dk, i denotes vertical

εin+1 − εin = Dε,i + Pε,i + Bε,i − Eε,i ,

t

(20)

where the terms on the right-hand side correspond to the shorthand notation used in Eq. (9). The diffusive term Dε, i is implemented in the same way as the term Dk, i above. The implementation of the other terms is done according to the following schemes:

Pεn,i = cε1 cμ kni Mi2 ,

Bnε,i = −cε3

Eεn,i = cε2

cμ

σρ

Ni2

⎧ n+1 ⎨kni εεi n if cε3 Ni2 > 0 i

⎩

2εin εin+1 kni

(21)

kni

− cε 2

if

(εin )2 kni

cε3 Ni2 .

(22)

<0

(23)

The potential sink terms again use positivity corrections as described above. The term Ek, i uses a Newton method to deal with the non-linearity. The Neumann boundary condition (13) for ε is implemented at half a grid cell above the bed. The discretised boundary condition reads

ε1 − ε0 |u∗b |3 =−

z 1 κ (z0 + 12 z1 )2

(24)

Y.M. Dijkstra et al. / Ocean Modelling 104 (2016) 129–142

2.3. Discretisation of the τ -equation The equation for τ reads

τ

n+1 i

− τin

t

= Dτ ,i + Pτ ,i + Bτ ,i − Eτ ,i + Dkτ ,i + Dτ τ ,i + Dkk,i ,

(25)

where the right-hand side terms correspond to the short-hand notation used in Eq. (11). The diffusive term Dτ , i again uses the same discretisation as the term Dk, i . The local production terms read

Pτn,i = cτ 1 cμ Mi2 τin τin+1 ,

 n

Bτ ,i =

−cτ 3 σμρ Ni2 τin τin+1 c −cτ 3 σμρ Ni2 (τin )2 c

(26)

if cτ 3 Ni2 > 0 if cτ 3 Ni2 < 0

(27)

and Eτ , i is equal to the constant coefficient cτ 2 . The additional diffusive type terms Dkτ and Dτ τ are implemented together according to a first-order upwind scheme

Dnkτ ,i

n

− Dτ τ ,i = −2

where

wtur

cμ

στ

 τ n+1 −τ n+1 i+1

wtur

i

zi+1 +1 τin+1 −τin−1

z i

if wtur < 0 if wtur > 0

,

n n  ∂τ  ∂ k  =− τ −k , ∂ z i ∂ z i  τ n + τin−1 kni − kni−1 1 τin+1 + τin kni+1 − kni =− + i 2 2

zi+1 2

z i  kni+1 + kni τin+1 − τin kni + kni−1 τin − τin−1 − − . 2

zi+1 2

z i

(28)



133

are dominated by boundary friction. The last two cases are mixing by unstable stratification and temperature modelling in an idealised lake model. These cases are dominated by buoyancy-induced mixing. They are included to test the performance of the k − τ model in cases where boundary friction is absent or weak. To compare the k − ε and k − τ models the sensitivity of the results to grid resolution is investigated for each of the four test cases, and the convergence behaviour is studies. This is done for a wide range of grid resolutions. In all cases an equidistant grid is used. The simulations additionally use the same time step for all resolutions. This chosen time step is the result of many tests with many time steps and is chosen such that numerical errors resulting from the time stepping procedure are small. The effect of the time step size is also studied separately for a particular vertical grid resolution. Whenever possible, the results are discussed qualitatively by explaining how these results follow from the turbulence model equations. A quantitative discussion of the results then follows by comparing the results at low resolutions, with the converged result at a high resolution. By using a 1DV model it is feasible to carry out tests for up to several thousands of layers, which is generally sufficient to obtain reasonably converged results. The results are also compared against theoretical results if these are available. In all cases, the model results are tuned to exactly fit the theoretical relations at a certain point in time by adjusting calibration parameters. The comparison with theory is therefore mainly to verify that the model is applied in a realistic parameter space. 3.1. Case 1: open channel flow

(29)

This upwind implementation is inspired by the observation that wtur has the dimension of a velocity. The first-order upwind method is preferred over a higher-order method, because it is more stable. A second-order central method has also been tested, but this proved to be an important source for unstable results, especially in cases involving strong vertical density gradients. The molecular viscosity is neglected in the implementation of the Dkτ term, because it is likely to cause model instabilities. The molecular eddy viscosity appears for example in a product ν 1k ∂∂ kz .

If at some point k tends to zero, then ∂∂ kz will also tend to zero, because of the positivity of k. However, the numerical approximation of ∂∂ kz is not necessarily small if k is small, so that this product can grow to infinity. The relative unimportance of the molecular viscosity compared to the typical magnitude of the eddy viscosity justifies neglecting this term. The term Dkk contains a similar numerically unstable product 1k ∂∂ kz . It has been chosen to neglect the Dkk term in the k − τ model, because of the stability problems. Neglecting the Dkk term creates a difference between the k − ε and k − τ models. This difference is only small for the present test cases, indicating that the effect of the Dkk term is rather small. Although this term may be significant, it is expected to be relatively small in general. This is because of the factor σε−1 − σk−1 , which equals −0.23 here. This factor is relatively small compared to 1, so that it contributes to Dkk being typically small compared to the other terms. In the results section we will verify the importance of neglecting this term by repeating the simulation using σε = σk , for which the Dkk cancels from the equation. 3. Test cases & results Four test cases are used to compare the k − ε and k − τ models. The test case were selected as they represent different types of flows. The first two cases are of a non-stratified open channel flow and an idealised estuarine flow. These cases represent flows which

An open channel flow is a non-stratified, steady, uniform flow which is driven by a constant pressure gradient and affected by bed friction. This case has also been used by Burchard et al. (1998) and Warner et al. (2005) to compare turbulence models which each implemented a somewhat different physics. Their analyses did not include tests at multiple resolutions, so that the convergence behaviour of the turbulence models is unknown. The analytical estimate of the depth-averaged velocity u in a sloping open channel flow with prescribed pressure gradient ∂ζ ∂ x is given by Chézy’s law



uChézy =

gh ∂ζ cf ∂x

1 / 2

,

(30)

where cf is a dimensionless friction parameter. This friction parameter is converted to a roughness height z0 by assuming an approximate logarithmic law-of-the-wall type velocity profile. This approximation is consistent with the logarithmic velocity profile as long as the roughness height is much smaller (i.e. more than a factor 100 smaller) than the depth. The conversion relation reads





/2 z0 = h exp −1 − κ c−1 , f

(31)

where κ is the Von Kármán constant. We will express the results in terms of a dimensionless average velocity q∗ which is defined as the average calculated velocity divided by the Chézy velocity, so that q∗ = 1 according to Chézy’s law. Although it is common practice to switch between cf and z0 by assuming a logarithmic profile, measurements indicate that this is only true by approximation and that this leads to an underestimate of the velocity (Ueda et al., 1977; Nezu and Rodi, 1986). A better estimate can be obtained by adding Coles’ wake law to the Chézy estimate (Coles, 1956). The dimensionless discharge can then be approximated by

q∗ = 1 +

c1f /2

κ

.

(32)

Nezu and Rodi (1986) advice a value of = 0.2. For this case we −5 −3 m (i.e. c = 2.8 · 10−3 ), use ∂ζ f ∂ x = 5 · 10 , h = 5 m and z0 = 1 · 10

134

Y.M. Dijkstra et al. / Ocean Modelling 104 (2016) 129–142

Open channel flow − layer convergence 1.1 Chézy’s law Wake law Π = 0.2 k−τ k−ε

1.08

q*

1.06 1.04 1.02 1 0.98

5

10

20

50

100 500 1000 Number of layers

1500

2000

5000

10000

Fig. 3. Dimensionless velocity plotted against the number of layers in the case of open channel flow (case 1). Note that the scale on the horizontal axis is non-linear.

which should result in an average velocity of just under 1 m/s and a q∗ of 1.025 assuming = 0.2. We will plot both the result from Chézy’s law and the wake law for = 0.2 to indicate an area where the numerical model results are expected based on theory. 3.1.1. Results The results are presented in Fig. 3. This figure shows the results of simulations with both turbulence models at a range of grid resolutions. The figure shows that the results of the k − τ model vary only by 1% between 5 and 10,0 0 0 layers, while the results of the k − ε model vary by 3%. In addition to this, the results of the k − ε model up to 50 layers seem to suggest a trend with q∗ ≈ 1.05. The result only starts to converge for resolutions above 50 layers and the result has converged within 1% after 500 layers. Such a resolution is often not feasible for large-scale 3D simulations. Additionally, the convergence behaviour of the k − τ model is not monotonous. However, the variation of q∗ with the grid resolution is smaller. It should be noted that the k − ε and k − τ models converge to slightly different values of q∗ , but this difference is less than 0.1% for 10,0 0 0 layers. This becomes more clear when considering the differences in the eddy viscosity predicted by both models. This shows differences between the models of the order of 10−4 m2 /s, which is small compared to the typical eddy viscosity of 10−2 m2 /s, but still considerably larger than the molecular viscosity of 10−6 m2 /s. These differences are caused by neglecting the Dkk term in the τ -equation (11). This has been verified by repeating the experiment for σε = σk , which makes the Dkk term cancel from the τ equation (also see end of Section 2.3). The difference between the eddy viscosity in both models reduces to the order of 10−6 m2 /s at 10,0 0 0 layers. The results are examined more closely by considering the profiles of k, τ , ε , ν t and u at 10 and 10,0 0 0 layers in Fig. 4. The upper-left panel shows that both turbulence models represent the profile of k almost as well at 10 layers, as at 10,0 0 0 layers. The profile of k in both models is also virtually identical for 10,0 0 0 layers. The upper-right panel shows the profiles of ε and τ , where ε has been rewritten in terms of τ in order to compare the two models. The profile of τ resulting from the k − τ model using 10 layers is almost equal to the converged result in the lower half of the water column. This is due to the linearity of the profile, the gradient of which is represented exactly by any first-order or higher-order linear numerical difference operator. In comparison, the profile of ε in the lower part of the water column is hyperbolic and less accurately approximated. When ε is rewritten to τ it shows that this value is too large at the bed and too small just above the bed. The opposite result is found in the upper half of the water column. The strongly curved profile of τ in this case is significantly underes-

timated at 10 layers, while the value of ε is fairly constant and therefore captured more accurately. This difference in τ near the surface remains, even for 10,0 0 0 layers. Similar results are found in the profile of the eddy viscosity (lower-left panel) at 10 layers, with the k − τ model being more accurate in the lower half of the water column and the k − ε model being more accurate in the upper half of the water column. The results at 10,0 0 0 layers are virtually the same in the lower half of the water column, but show some remaining differences in the upper half. Although no clear preference for either of the models follows from the eddy viscosity profiles, the velocity profile obtained by using the k − τ model with 10 layers is significantly better than that obtained by using the k − ε model at the same resolution. These velocity profiles, shown in the lower-right panel of Fig. 4, can be explained by considering the relation between the velocity and eddy viscosity. This relation reads R = νt ∂∂uz , where R is the Reynolds stress. The Reynolds stress computed from the model is free of numerical errors. This is because the Reynolds stress is solved for directly in the model by using the conservative form of the momentum equations and because the Reynolds stress is linear in the vertical direction. As a result, an error in the computation of the eddy viscosity will therefore directly result in an error in the velocity gradient. Because the velocity gradient is largest near the bed, a small error in the eddy viscosity near the bed will cause a significant error in the velocity gradient. Conversely, because the velocity gradient is small in the upper half of the water column, even a significant error in the eddy viscosity will only cause a small error in the small velocity gradient there. As a consequence it is more important to have a turbulence model that is accurate near the frictional bottom boundary than near the surface. This explains why the velocity profile obtained by using 10 layers is more accurately modelled when using the k − τ model, than when using the k − ε model. The velocity profiles from both models at 10,0 0 0 layers show hardly any difference. 3.2. Case 2: estuarine flow Next, we will consider an idealised estuarine flow. The velocity consists of a single M2 tidal component and a prescribed, vertically uniform horizontal salinity gradient. The simulation starts with a vertically uniform salinity distribution. As the tidal flow velocity develops, it deforms the vertical salinity profile, so that fresh water tends to move over salt water during ebb and vice versa during flood, see Fig. 5. This process is called strain-induced periodic stratification (SIPS, see e.g. Simpson et al., 1990). As a result of the density straining, the water column tends to be stratified during ebb and well-mixed during flood. This induces a significant variation in turbulent mixing over the tidal cycle due to buoyancy effects.

Y.M. Dijkstra et al. / Ocean Modelling 104 (2016) 129–142

τ and ε rewritten to τ 5

4

4

3

3

z (m)

z (m)

TKE 5

2 1 0

2 1

0

0.002

0.004

0.006

0.008

0

0.01

0

20

60

80

100

1

1.2

Velocity

5

5

4

4

3

3

z (m)

z (m)

40 τ (s)

2 2

k (m /s ) Eddy viscosity

2 1 0

135

k−τ − 10,000 layers k−τ − 10 layers k−ε − 10000 layers k−ε − 10 layers

2 1

0

0.005

0.01

0.015

0.02

0.025

2

ν (m /s)

0 0.2

0.4

0.6 0.8 u (m/s)

t

Fig. 4. Vertical profiles of turbulent kinetic energy (TKE), τ (note that for comparison ε has been rewritten in terms of τ ), eddy viscosity and velocity in the open channel flow case (case 1) with 10 layers compared to 10,0 0 0 layers.

Fig. 5. Schematised side view of estuarine flow, where blue indicates salt water and white indicates fresh water. The 1DV model simulates a single water column in the middle of this domain. Left: the ebb tidal flow in combination with the horizontal salinity gradient strains fresh water over salt water. Middle: the flood tidal flow tends to strain salt water over fresh water, but turbulent mixing prevents the formation of strong unstable density gradients. Right: The density gradient induces a sub-tidal exchange flow which is typically directed landwards near the bed and seawards near the surface. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Similar model set-ups have been used by Nunez Vaz and Simpson (1994), Monismith and Fong (1996) and Burchard and Hetland (2010). The depth-averaged velocity is prescribed according to u = U sin(ωt ), where ω = 1.4 · 10−4 1/s is the angular frequency corresponding to the M2 tide and U = 1 m/s is the amplitude. The model uses a fixed horizontal buoyancy gradient ∂∂ bx = 10−5 1/s2 , ρ −ρ

where the buoyancy is defined as b = −g ρ 0 . The density in this 0 expression is linked to the salinity through a linear equation of state: ρ = ρ0 (1 + β s ), with β = 7.6 · 10−4 1/psu. The depth is 8 m and the roughness height is z0 = 1.5 · 10−3 m, which corresponds approximately to c f = 2.8 · 10−3 according to Expression (31). The model is run for 4 tidal cycles, where the fourth cycle is used to compute the results. These results differ only marginally from those of the third tidal cycle, so that the results are for a fully developed periodic situation.

The turbulence models will be compared by the magnitude of the exchange flow (see Fig. 5). The exchange flow is a residual flow with depth-averaged value of zero, seaward flow near the surface and landward flow near the bed. This exchange flow is considered to be important for the modelling and understanding of salt and sediment transport (e.g. Burchard and Hetland, 2010). The exchange flow can be separated into a part induced by the horizontal density gradient (gravitational circulation, see e.g. Hansen and Rattray, 1965) and a part due to tidal correlations of the eddy viscosity and velocity (also known as straining circulation, see e.g. Jay and Musiak, 1994; Burchard and Hetland, 2010). The gravitational circulation is sensitive to the magnitude of the time-averaged eddy viscosity, while the straining circulation depends strongly on the temporal variations of the eddy viscosity. The exchange flow magnitude is therefore a good means to evaluate the accuracy of the numerical implementation of the turbulence models. The magnitude of the exchange flow is measured by taking the depth-average

136

Y.M. Dijkstra et al. / Ocean Modelling 104 (2016) 129–142

Depth−averaged absolute exchange flow − time step convergence using 20 layers 0.07 k−τ k−ε

|u| (m/s)

0.068 0.066 0.064 0.062 0.06

5

10

20

50

100 500 Number of layers

1000

1500

2000

Fig. 6. The absolute depth-averaged exchange flow in the estuarine case (case 2) plotted against the grid resolution for a time step of 1 s. Note that the scale on the horizontal axis is non-linear.

of the absolute velocity;

uexchange

1 = H



0 −H

|uexchange | dz.

(33)

3.2.1. Results Fig. 6 shows the exchange flow magnitude against a range of grid resolutions for a fixed time step of 1 s. The k − ε model shows a better result than the k − τ model at 5 layers. Here a result is considered better than another when it is closer to the converged result. However, the k − ε model then shows a worsening of the result as the number of layers increases to 20. The k − τ model on the other hand converges almost monotonously, with its result already being better than that of the k − ε model at 10 layers. The difference between the model becomes significant at 20 layers, where the k − τ model differs by 5% from the converged result and the k − ε model differs by 9%. At 50 layers the difference between the models has increased to a 2% deviation of the k − τ model and an 8% deviation of the k − ε model. The k − ε model results only start to converge at a high rate for more than 100 layers. It should be noted that the models converge to results which differ by approximately 1% of the exchange flow magnitude when using 20 0 0 layers. One of the reasons for this difference is the neglected Dkk term. Setting σε = στ = 1 to eliminate this term indeed reduces the difference to 0.5%. The average difference in the eddy viscosity between the two models at 20 0 0 layers is of the order of 10−5 m2 /s, with a maximum difference in time and space of 10−4 m2 /s. These difference are quite small compared to the typical eddy viscosity magnitude of 10−2 m2 /s, but large compared to the molecular viscosity. The most important source of differences in the eddy viscosity is a couple of cells difference in the position of density gradients. One possible reason for the differences is the discrete switch in the value of cε3 and cτ 3 from stable to unstable stratification. Such a switch can amplify the smallest of numerical differences to a significant difference. Another possible reason is the minimum value of 10−32 for k, ε and τ , which is set to prevent divisions by zero. Such a truncation can affect the results especially around sharp density gradients, where the eddy viscosity is close to zero. The better performance of the k − τ model at most low resolutions is again explained from the linear profile of τ near the bed. Fig. 7 shows time series of the depth-averaged shear production and buoyancy production of turbulence in the k − τ model with 100 layers and a time step of 1 s. The depth-averaged shear production is at least one order of magnitude larger than the depthaveraged buoyancy production. This dominance of the shear production comes from its dominant role in the lower half of the water column, while shear and buoyancy are of an equal order of magnitude in the upper half of the water column. The turbulence near the bed therefore behaves similar to a non-stratified flow case

and the buoyancy mainly affects the upper half of the water column. Concluding, the conclusions from the open channel flow case, regarding the higher accuracy of the k − τ model in the lower half of the water column, also apply to the present case. We will proceed by investigating the role of the time step in the accuracy of the turbulence model. This is done by fixing the vertical resolution at 20 layers, a resolution that is commonly used in coastal ocean modelling, and by varying the time step. The results are presented in Fig. 8. Neither model shows an instability even for time steps of a minute or more, although the results do become inaccurate for such time steps. Overall, both models show a similar convergence behaviour and error for time steps smaller than 20 s. The results of both models have converged to within 1% at a time step of 10 s and within 0.1% at a time step of 1 s.

3.3. Case 3: mixing by unstable stratification The previous cases show the performance of both turbulence models in flows dominated by a frictional boundary layer. The next cases present results of cases that are dominated by buoyancy effects. The first is of mixing by unstable stratification, also called free convection, which is relevant to many coastal ocean applications. Examples are surface cooling in oceans and lakes and shearing of salt water over fresh water in coastal seas and estuaries. The experiments by Deardorff et al. (1969) and Willis and Deardorff (1974) provide a simple case for testing the ability of a turbulence model to simulate such free convection. These experiments have also been repeated in a numerical 1DV model by Burchard and Petersen (1999). It has been remarked that the two-equation turbulence models, such as the ones used here, do not provide an accurate representation of free convection (e.g. Burchard and Petersen, 1999). The reason is that these models do not contain the physics of convection. Convection is a fundamentally non-hydrostatic process that can only be parametrised with hydrostatic models. This parametrisation is incorporated in the two-equation turbulence model and uses a combination of a diffusion process and local buoyancy production to simulate the effect of free convection. The numerical study by Burchard and Petersen (1999) shows that the k − ε model is nonetheless capable of reproducing the measured density profiles in a simple case of free convection. We will repeat experiment A of Deardorff et al. (1969) in the numerical model. The experiment starts with a tank of water with a surface temperature Ts and a bottom temperature Tb,stable < Ts so that a stationary stably stratified state is reached. Then, at t = 0, the temperature under the bottom of the tank is changed to Tb,unstable > Ts by letting hot water flow under the tank. This

Y.M. Dijkstra et al. / Ocean Modelling 104 (2016) 129–142 −4

x 10

Depth−averaged shear production 1

137

−5 x 10 Depth−averaged buoyancy production

0

k

B (m2/s3)

Pk (m2/s3)

0.5 2

1

−0.5

0

0

0.2

0.4

0.6

0.8

−1

1

0

0.2

0.4

0.6

t/T

0.8

1

t/T

Fig. 7. Time series of depth-averaged shear production (left) and buoyancy production (right) in the estuarine flow case (case 2) with the k − τ model at 100 layers and a time step of 1 s. The horizontal axis shows one tidal cycle with tidal slack at t/T = 0 and t/T = 0.5, peak flood at t/T = 0.25 and peak ebb at t/T = 0.75.

Depth−averaged absolute exchange flow − time step convergence using 20 layers 0.08 0.075 |u| (m/s)

0.07 0.065 0.06 k−τ k−ε

0.055 0.05 120

60

30

20

10 Time step (s)

5

1

0.5

0.1

Fig. 8. The absolute depth-averaged exchange flow in the estuarine case (case 2) plotted against the time step using 20 grid cells. Note that the scale on the horizontal axis is non-linear.

forces convective mixing in the water column and a mixed-layer starts to develop. The experiment uses a depth of 0.36 m and temperatures Ts = 39◦ C, Tb,stable = 21.5◦ C and Tb, unstable = 43.4◦ C. The mixed-layer depth is defined as the level of the upper side of the thermocline. Numerically this is the lowermost cell interface with a TKE level below 10−8 m2 /s2 . The diffusion of heat through the bottom of the tank is modelled numerically by adding two artificial cells below the bed level. A heat diffusion equation is solved in these cells to compute the diffusion through the bottom of the tank. The thickness of the cells and rate of diffusion are used as calibration parameters to approximately match the model results to the measurements. The diffusion coefficient in the artificial cells is 7.5 · 10−8 m2 /s and the thickness of each cell is 0.001 m. 3.3.1. Results Fig. 9 shows the evolution of the mixed-layer depth in time using 20 and 100 layers and a time step of 0.5 s. The figure shows that the resolution affects the results significantly, especially at the start of the simulation. When using 20 layers or less, both turbulence models need some time before the mixed layer starts to develop. This delay disappears when using a resolution of 100 layers. The reason for the delay in turbulence production at low resolutions is best understood from the buoyancy production term: Bk = σνρt ρg ∂ρ ∂ z . In this expression, the eddy viscosity is almost zero initially. Additionally, ∂ρ ∂ z is small, because z is large on coarse grids. The buoyancy production is therefore small and turbulence production takes time to initiate. The delay of the mixed-layer development disappears when the resolution is increased. The

increased resolution allows for a large unstable density gradient to develop near the bottom of the tank. This high gradient accelerates the production, so that the mixed-layer develops faster. At 100 layers, both models reproduce the theoretical relation that is suggested by Deardorff et al. (1969) well, see Fig. 9. The k − τ model is more accurate than the k − ε model, but the results do highly depend on the way the diffusion of heat through the bed is calibrated. Several model experiments have shown that a different calibration of the heat diffusion through the bed would lead to a closer match between the k − ε model and the theoretical result. At the same time it would decrease the match between the k − τ model and the theoretical result. It can therefore be concluded that the k − τ and k − ε behave equivalently well at a resolution of 100 layers. When the resolution is increased further to 500 layers, the k − ε model again displays a delay in the mixed-layer development. The most important reason for this delay is the use of a Neumann boundary condition for ε at the bed. This boundary condition cannot accurately approximate the value of ε at the bed. As a result, the value of the eddy viscosity is very small at the bed, so that the production of turbulence is retarded. This has been verified by testing the k − ε model with a Dirichlet boundary condition, which indeed shows no delay in turbulence production for high resolutions. Similarly, the k − τ shows no delay of production at high resolutions, because this model uses a Dirichlet boundary condition for τ . It should be noted that the described delay is an artefact of this particular model case that will rarely be encountered in more realistic modelling cases. This is because it requires a flow to be fully at rest initially. The Neumann and Dirichlet boundary

138

Y.M. Dijkstra et al. / Ocean Modelling 104 (2016) 129–142

Fig. 9. Deepening of the mixed layer in time in case 3. Left: using 20 layers. Right: using 100 layers.

Deardorff experiment A after 800 seconds − layer convergence −0.05

Dm (m)

−0.1 −0.15 −0.2 −0.25

k−τ k−ε − Neumann BC k−ε − Dirichlet BC

−0.3 −0.35 5

10

20

50

100 500 Number of layers

1000

1500

2000

Fig. 10. Depth of the mixed layer after 800 s plotted against the grid resolution using a time step of 0.5 s in the free convection case (case 3). The error bars denote the cell thickness. Note that the scale on the horizontal axis is non-linear.

conditions are technically inappropriate for such conditions. This problem is resolved by any disturbance in the form of a weak background flow. Indeed, a background flow of 1 mm/s or more is sufficient to start mixing, even when only using 5 layers. The convergence of the mixed-layer depth after 800 s with increasing grid resolution is plotted in Fig. 10. The figure displays the results of the k − τ model and k − ε model with both Neumann and Dirichlet boundary condition (see e.g. Burchard and Petersen, 1999). It is shown that, if only 5 or 10 layers are used, turbulence has not started to develop after 800 s, because of the reasons explained above. At resolutions of 500 layers or more, the results of the Neumann boundary condition of the k − ε model starts retarding turbulence production, leading to a reduction of the mixed layer depth. The k − ε model with Dirichlet boundary condition does not suffer from such problems and produces a fairly constant result for all resolutions. The total variation of the result is 15% in the k − ε model with Neumann boundary condition, 4% in the k − ε model with Dirichlet boundary condition and 8% in the k − τ model between 20 and 20 0 0 layers. The boundary condition of the k − ε model thus determines a large part of the different converged results between both turbulence models at 20 0 0 layers. The remaining difference is small; with the eddy viscosity between the models differing a maximally 3 · 10−5 m2 /s at t = 800 s using 2000 layers and a time step of 0.5 s. The remaining difference can be explained from the Dkk term and the neglected effect of the molecular viscosity in the Dkτ and Dτ τ terms in the k − τ model. If the experiments are repeated for cancelling Dkk term and zero molecular viscosity, the maximum eddy viscosity difference reduces to 7 · 10−6 m2 /s at t = 800 s using 20 0 0 layers and a time step of 0.5 s. The results still seem to be converging further for increasing resolution.

The time step also has a clear effect on the development of the mixed-layer depth. Fig. 11 shows the mixed-layer depth simulated using several time steps and using 20 0 0 grid cells. Both models under-predict the mixed-layer depth for large time steps. The most likely explanation for this is that non-linear effects in the model cannot be properly evaluated if the time step is too large. Both turbulence models have been linearised for the numerical implementation and need a sufficiently small time step to yield an accurate approximation of the non-linear equations. The results show that the k − ε model is the less accurate of the two models for small time steps and has a more erratic convergence behaviour than the k − τ model. It is difficult to analyse the reasons why the k − ε and k − τ models show a different response to large time steps. This depends strongly on the way the k, τ and ε equations are coupled in the numerical implementation. For time steps smaller than 2 s, both models show a fairly converged result.

3.4. Case 4: mixing processes in a lake The final case is of mixing processes in an idealised lake model. This 1DV schematisation of a water column in a lake is forced externally by prescribed atmospheric conditions. The turbulence models calculate the mixing processes related to stable and unstable buoyancy gradients and a weak wind induced flow. The basis of this study is the study by Uittenbogaard and AparicioMedrano (2011), who modelled the temperature development in a 1DV model of Lake Vlietland, the Netherlands, based on measured wind, air temperature and cloud cover data. These data are incorporated into the 1DV model using the water-atmosphere model of Delft 3D. This model determines several heat flux components between the water column and the atmosphere. We refer readers to Deltares (2014) for details on this model.

Y.M. Dijkstra et al. / Ocean Modelling 104 (2016) 129–142

139

Deardorff experiment A after 800 seconds − time step convergence using 2000 layers −0.05

D (m) m

−0.1 −0.15 −0.2 −0.25 k−τ k−ε

−0.3 −0.35 50

20

10

5

2 Time step (s)

1

0.5

0.2

0.1

Fig. 11. Depth of the mixed layer after 800 s plotted against the time step using 20 0 0 grid cells in the free convection case (case 3). The error bars denote the cell thickness. Note that the scale on the horizontal axis is non-linear.

Simulated temperature (Celsius) over a year − k−τ model

z (m)

30

20

20 10

10 0

50

100

150 200 250 300 Time (days) Simulated temperature (Celsius) over a year − k−ε model

350

z (m)

30

0

20

20 10

10 0

50

100

150 200 Time (days)

250

300

350

0

Fig. 12. Modelled temperature in the water column over the span of a year in the idealised lake case (case 4). Top: results when using the k − τ model. Bottom: results when using the k − ε model. Both models use 20 0 0 layers and a time step of 0.1 min.

In their study, Uittenbogaard and Aparicio-Medrano (2011) compared the modelled temperature at 10.00 AM every day for a year to measurements. They show that the 1DV model can model the temperature evolution reasonably well. It captures the main seasonal thermocline evolution and several smaller mixing events well. In this study, we do not compare the results against their measurements, as these are too coarse in space to differentiate clearly between the model results at very high resolution. Instead we compare the results of all tested vertical resolutions against the reference result using the k − ε model with 20 0 0 layers and a 0.1 min time step. Our goal is not to test the validity of the model results in practice, but to investigate the sensitivity to vertical grid resolution and time step. The 1DV model uses a depth of 36 m. The only other parameters of relevance are the Stanton and Dalton numbers that are calibration parameters of the water-atmosphere model. These parameters are set to values 1.3 · 10−3 and 1.0 · 10−3 , corresponding to the study by Uittenbogaard and Aparicio-Medrano (2011).

3.4.1. Results A typical result of the modelled temperature over a year is shown in Fig. 12, which shows the computed temperature using 20 0 0 layers and a time step of 0.1 min both models. The lake starts to stratify in late April and develops the main thermocline at a depth of about 15 m. Over the summer, surface temperatures reach up to 25°, while the deepest water remains at around 7°. The thermocline disappears again in late autumn due to a combination

of autumn storms and lower air temperatures. The results of the two models show only minor differences. The convergence of both turbulence models is measured by comparing the computed temperatures using different grid resolutions against the k − ε model with 20 0 0 layers. We compare the results of a simulation using n layers (n ≤ 20 0 0) to this reference case by interpolating the reference case results onto the n-layer grid. We then calculate the root-mean square (RMS) of the difference averaged over space and time. The resulting RMS errors are presented in Fig. 13. Both models show a similar, and fairly small, degree of variation of the error with the grid resolution between 50 and 20 0 0 layers. The k − τ model results show a large degree of variation for very low resolutions of 5 to 20 layers. However, the k − τ and k − ε models show a similar degree of variation for resolutions exceeding 20 layers. The k − τ model at 50 layers has an RMS error that differs 0.07° from the its 20 0 0 layer result. The k − ε model at 50 layers has a 0.12° RMS error. It is difficult to conclude from these numbers whether one model is less dependent on the grid resolution than the other, as this depends on the exact error measure that is used. The final convergence of the models on the detailed scale is slow; the RMS error shows that the results have not fully converged when using 20 0 0 layers. The models converge slightly different results. Note here that the k − ε model converges to a zero error by definition, as the reference case uses the k − ε model. The simulation has been repeated using σε = σk so that the Dkk term cancels from the equation. This yields a similar convergence behaviour for both models and the RMS error of the k − τ model at

140

Y.M. Dijkstra et al. / Ocean Modelling 104 (2016) 129–142

RMS Error of temperature compared to reference solution

Error (deg. Celsius)

2 k−τ k−ε

1.5

1

0.5

0

5

10

20

50

100 500 Number of layers

1000

1500

2000

Error (deg. Celsius)

Fig. 13. Root-mean square average difference between the temperature computed using different grid resolutions, compared to the results when using the k − ε model with 20 0 0 layers in the idealised lake case (case 4). All simulations use a time step of 0.1 min. Note that the scale on the horizontal axis is non-linear.

RMS Error of temperature compared to reference solution − time step convergence using 100 layers 2.5 k−τ k−ε 2 1.5 1 0.5 0 30

10

5

2

1 Time step (min)

0.5

0.2

0.1

0.05

Fig. 14. Root-mean square average difference between the temperature computed using different time steps, compared to the results when using the k − ε model with 20 0 0 layers in the idealised lake case (case 4). All simulations use a grid resolution of 100 layers. Note that the scale on the horizontal axis is non-linear.

20 0 0 layers relative to the k − ε model at 20 0 0 layers decreases to 0.04°. Visually, there are no major differences to the temperature profiles presented in Fig. 12. This means that the Dkk term, in this case, does not result in a significant qualitative difference in the temperature profile. Nevertheless it is quantitatively relevant on the level of small details, as eliminating it reduces the RMS error by approximately a factor 8. Even though the RMS difference between the models is only 0.04°, the average spatial and temporal difference in eddy viscosity is still 10−5 m2 /s. The maximum eddy viscosity difference is even 10−3 m2 /s. These differences originate from very small differences between the model results, where sharp density gradients have a different position by only a couple of grid cells. The convergence of the models with the time step and a fixed resolution of 100 layers is plotted in Fig. 14. The convergence is monotonous and quite similar for both models. At 100 layers, the results seem fairly well converged already for a time step of 1 min. Some small improvements can be seen when decreasing the time step further to 0.2 min. It should be noted that although a time step of 1 min might be sufficient for 100 layers, this might not be sufficient for the convergence of the results at resolutions of 10 0 0 to 20 0 0 layers. It should also be noted that the requirements on the time step and grid resolution for reaching converged results are high compared to typical numbers used in practice. For example the original study by Uittenbogaard and Aparicio-Medrano (2011) uses 200 layers and a time step of 5 min. 3.5. Computation time Finally we present the differences in computation time between the k − ε and k − τ models. It is expected that the k − τ model takes more computation time per run, because of the additional diffusion terms that need to be computed. Apart from these diffu-

sion terms, the numerical implementation of the models is similar. The computation time is determined by recording the total time consumed by the k − ε and k − τ routines in the estuarine flow case (Section 3.2). The computation time is thus purely the time consumed by the turbulence model. The results are for grid resolutions between 5 and 20 0 0 layers and a fixed time step of 1 s. The simulations over all grid resolutions are repeated three times in order to estimate both the mean and standard deviation of the computation time. The computation time of the k − τ model, relative to the k − ε model is 1.19 with a standard deviation of 0.03, i.e. the k − τ model takes typically 19% ± 3% more computation time than the k − ε model. 4. Discussion The k − ε turbulence model is one of the most popular twoequation turbulence models in coastal ocean models. Many studies have shown that the model is able to represent turbulent mixing fairly well in a wide range of applications. However, the numerical properties of the k − ε model in boundary layers are not ideal. The most important reason for this is related to the hyperbolic shape of ε in non-stratified shear flows, as has been illustrated in Fig. 1. By transforming the k − ε model to the k − τ model one obtains a model that has the same physical properties, but more favourable numerical properties in boundary layer flows. One can construct an infinite number of turbulence models with different length-scale parameters of the form km ln , with arbitrary m and n (Launder and Spalding, 1974; Umlauf and Burchard, 2003). Models with different choices of length-scales, such as ε ∼ k3/2 l −1 , ω ∼ k1/2 l −1 and kl are typically accompanied by choices relating to the calibration of the model coefficients and parametrisation of the diffusion of turbulence. As a result, different turbulence models generally have different physical and numerical

Y.M. Dijkstra et al. / Ocean Modelling 104 (2016) 129–142

properties, which both affect the results to a similar extent. Consequently, most comparisons of turbulence models in literature show a combination of the physical and numerical differences, without distinguishing between the two. Here we take a different approach to the length-scale variable, by using it to change only the numerical properties of a turbulence model. The transformation of the length-scale variable used here adds a number of diffusive-type terms to the equation, so that it represents the same physical schematisation as the original lengthscale variable. The comparison of different turbulence now allows us to highlight the numerical differences between the models. The use of transformations thus supports the conclusions by, amongst others, Launder and Spalding (1974) and Speziale et al. (1992) that the choice of length-scale variable is only of a limited physical significance. Instead of transforming ε to τ , one could transform it to any general length-scale variable km ln . We have chosen τ = k−1/2 l 1 here, because of its linearity in the frictional boundary layer, but in fact any length-scale using n = 1 has this property. The reason for choosing τ specifically is that the profile of τ keeps its linearity for a large part of the water column, while variables like kl (Mellor and Yamada, 1982) or kτ (Zeierman and Wolfshtein, 1986) have a more parabolic shape. The equations for such variables have similar diffusive terms and so offer no clear advantages over τ . It may nevertheless be that the k − τ model is not the single best possible transformation. This is sufficient for our study as we only aim to show that a transformation of the length-scale variable can have beneficial effects on the numerical accuracy and convergence behaviour. The results of this research show that the convergence behaviour of a turbulence model is a relevant subject. The results of the four 1DV cases only converge between 100 and 10 0 0 grid cells in the vertical direction, which is generally beyond the computational feasibility of complex 3D models. It is generally difficult to estimate the magnitude of the numerical error of such models. Assuming that the convergence study in 1DV at least provides an indication of the convergence in 3D, this 1DV study shows that the numerical accuracy of a turbulence model plays a significant part in the result of large-scale modelling studies. Additionally, he convergence of the models is not necessarily monotonous, so that a grid refinement will not necessarily reduce the numerical error. The linearity of τ in frictional boundary layers is especially beneficial at resolutions of 10 to 100 vertical grid cells, which are common in coastal ocean modelling. For such resolutions and in cases dominated by boundary friction, the k − τ model can achieve an accuracy that can only be attained with the k − ε model when using much more than 100 vertical grid cells. The time step has little influence on this in our test cases; the k − τ model does not require a smaller time step than the k − ε model. The price is the computation time of the k − τ model, which is about 19% ± 3% more than that of the k − ε model. From one perspective one could argue that this extra computation time is more than compensated for, because some simulations require more than a factor ten fewer grid cells to achieve the same accuracy as the k − ε model. From a different perspective one could argue that the gain in accuracy is not significant enough to justify the additional computation time. The k − τ model typically leads to a 3% to 5% more accurate approximation of the average velocity in our test cases of boundary friction dominated flows. Such an improvement may be insignificant for certain applications, but may well be important when modelling the transport of constituents, such as salt or sediments, which often scales with a power two or three of the velocity. It thus depends on the perspective and application as to whether or not it is economical to use the k − τ model instead of the k − ε model.

141

Apart from the economic considerations, it is also important to consider the predictability and stability of the turbulence models. The different boundary conditions of the models play a role in this. The Neumann boundary condition for the k − ε model is recommended, because it is more accurate than the Dirichlet boundary condition (Burchard and Petersen, 1999). However, the Neumann condition sometimes has some negative side-effects. An example is given in case 3, representing a case of convective mixing, where the Neumann boundary condition leads to problems in initiating turbulence from laminar conditions. Similar problems are not encountered with the more natural Dirichlet boundary condition of the k − τ model. Although problems with the boundary conditions of ε are only encountered in some specific cases, the natural boundary condition of τ is preferable as fewer or no problems are expected with this boundary condition. The three additional diffusive terms in the τ equation with respect to the ε equation, could pose a risk to the stability of the k − τ model. We have chosen to implement the Dkτ and Dτ τ terms in Eq. (11) according to a first-order upwind method, which leads to numerical diffusion and avoids numerical instabilities. The Dkk term has been neglected, because of the risk of numerical instabilities. Neglecting this term leads to generally small, but definite, differences between the results of the k − τ and k − ε models. As a result of the potential problems and inaccuracies caused by diffusive terms, the implementation of a transformed turbulence model is a non-trivial task. The comparison of two transformed models is certainly not trivial, as complete convergence of the models to the same result is difficult to achieve. Even after compensating for the Dkk term, the average difference between the eddy viscosity in both models is typically 10−5 m2 /s at 20 0 0 layers and small time step. Such a difference is small compared to the typical eddy viscosity, but an order of magnitude larger than the molecular viscosity. There are several possible explanations for these differences. Firstly, the results have not fully converged for 20 0 0 layers and the chosen time step. It has been demonstrated in the open channel flow case that an increase of the resolution from 20 0 0 to 10.0 0 0 layers makes a difference on this level of detail. Secondly, the buoyancy parameters cε3 and cτ 3 switch values between stable and unstable stratification in the parametrisation of Delft 3D-FLOW. Such a switch can amplify small differences between the models. Thirdly, k, τ and ε have a prescribed minimum value of 10−32 in order to avoid division by zero. It has been established that the code sometimes switches to these minima, thus creating a difference between the two models. It is unclear whether other factors such as solver precision or machine precision play a role in the differences between the model results. Summarising the results, the k − τ model as a transformation of the k − ε model does show improvements in the test cases and the metrics considered in this research. While boundary flows in general remain to be investigated, the results of our cases strongly suggest that the k − τ model shows some improvement of the numerical accuracy compared to the k − ε model in both stratified and non-stratified boundary layer flows. The 1DV cases of buoyancy dominated flows show that the k − τ and k − ε perform similarly in the absence of a frictional boundary layer, so that the k − τ model can be applied to both buoyancy dominated and boundary layer flows. Advantages of the k − τ model are the profile of τ , which is approximated more accurately in boundary friction dominated flows than the profile of ε , and the boundary condition of τ , which is better behaved. The main disadvantage of the k − τ model is the larger computational time. Given the widespread use of the k − ε model in coastal ocean models, the transformation to the k − τ model can be a valuable addition for a range of cases and applications.

142

Y.M. Dijkstra et al. / Ocean Modelling 104 (2016) 129–142

Acknowledgements This work was supported by Deltares. We would like to thank Henk Schuttelaars, Jan Mooiman and Herman Kernkamp for providing helpful ideas and comments. We also would like to thank an anonymous reviewer for many suggestions that strongly improved this work. References Burchard, H., Hetland, R.D., 2010. Quantifying the contributions of tidal straining and gravitational circulation to residual circulation in periodically stratified tidal estuaries. J. Phys. Oceanogr. 40, 1243–1262. Burchard, H., Petersen, O., 1999. Models of turbulence in the marine environment – a comparative study of two-equations turbulence models. J. Mar. Syst. 21, 29–53. Burchard, H., Petersen, O., Rippeth, T.P., 1998. Comparing the performance of the Mellor-Yamada and the k − ε two-equation turbulence models. J. Geophys. Res. 103 (C5), 10543–10554. Coles, D., 1956. The law of the wake in turbulent boundary layers. J. Fluid Mech. 1, 191–226. Deardorff, J.W., Willis, G.E., Lilly, D.K., 1969. Laboratory investigation of non-steady penetrative convection. J. Fluid Mech. 35, 7–31. Deltares (2014). Delft 3d user manual. Version 3.15.34158. Ezer, T., Arango, H.G., Shchepetkin, A.F., 2002. Development in terrain-following ocean models: intercomparisons of numerical aspects. Ocean Model. 4, 249–267. Ezer, T., Mellor, G.L., 2004. A generalized coordinate ocean model and a comparison of the bottom boundary layer dynamics in terrain-following and in z-level grids. Ocean Model. 6, 379–403. ´ K., 1994. Advanced turbulence closure models: a view of current status Hanjalic, and future prospects. Int. J. Heat Mass Transf. 15, 178–203. Hansen, D.V., Rattray, M., 1965. Gravitational circulation in straits and estuaries. J. Mar. Res. 23, 104–122. Ilicak, M., Özgökmen, T.M., Peters, H., Baumert, H.Z., Iskandarani, M., 2008. Performance of two-equation turbulence closures in three-dimensional simulations of the red sea overflow. Ocean Model. 24, 112–139. Jay, D.A., Musiak, J.D., 1994. Particle trapping in estuarine tidal flows. J. Geophys. Res. – Oceans 99, 20445–20461. Jones, W.P., Launder, B.E., 1972. The prediction of laminarization with a two-equation model of turbulence. Int. J. Heat Mass Transf. 15, 301–314. Launder, B.E., Spalding, D.B., 1974. The numerical computation of turbulent flows. Comput. Methods Appl. Mech. Eng. 3, 269–289. Martin, P.J., 1985. Simulation of the mixed layer at OWS November and Papa with several models. J. Geophys. Res. 90, 903–916. Mellor, G.L., Yamada, T., 1982. Development of a turbulence closure model for geophysical fluid problems. Rev. Geophys. Space Phys. 20, 851–875.

Mohammadi, B., Pironneau, O., 1994. Analysis of the K-epsilon Turbulence Model. Wiley. Monismith, S.G., Fong, D.A., 1996. A simple model of mixing in stratified tidal flows. J. Geophys. Res. 101, 28583–28595. Nezu, I., Rodi, W., 1986. Open-channel flow measurements with a laser doppler anemometer. J. Hydraulic Eng. 112, 335–355. Nunez Vaz, R.A., Simpson, J.H., 1994. Turbulence closure modeling of estuarine stratification. J. Geophys. Res. 99, 16143–16160. Patankar, S.V., 1980. Numerical Heat Transfer and Fluid Flow. Hemisphere Publishing Corporation. Platzek, F.W., Stelling, G.S., Jankowski, J.A., Pietrzak, J.D., 2014. Accurate vertical profiles of turbulent flow in z-layer models. Water Resources Res. 50, 2191–2211. Rodi, W., 1993. Turbulence Models and Their Application in Hydraulics - a State of the Art Review, 3rd IAHR, Delft. Balkema Simpson, J.H., Brown, J., Matthews, J., Allen, G., 1990. Tidal straining, density currents, and stirring in the control of estuarine stratification. Estuaries 13, 125–132. Speziale, C.G., Abid, R., Anderson, E.C., 1992. Critical evaluation of two-equation models for near-wall turbulence. AIAA J. 30, 324–331. Stelling, G.S., 1995. Compact differencing for stratified free surface flow. In: Advances in Hydro-science and -Engineering, 2, pp. 378–386. Thangam, A., Abid, R., Speziale, C.G., 1992. Application of a new k − τ model to near wall turbulent flows. AIAA J. 30, 552–554. Ueda, H.R., Möller, R., Komori, S., Mizushina, T., 1977. Eddy diffusivity near the free surface of open channel flow. Int. J. Heat Mass Transf. 20, 1127–1136. Uittenbogaard, R.E., Aparicio-Medrano, E., 2011. Hyrdodynamic simulations of the reduction of microcystis blooms by bubble plumes in deep lakes. Technical Report. Deltares. Uittenbogaard, R.E., Van Kester, J.A.T.M., Stelling, G.S., 1992. Implementation of three turbulence models in 3D-TRISULA for rectangular grids. Technical Report. WL | Delft Hydraulics, the Netherlands. Umlauf, L., Burchard, H., 2003. A generic length-scale equation for geophysical turbulence models. J. Mar. Res. 61, 235–265. Umlauf, L., Burchard, H., Hutter, K., 2003. Extending the k − ω turbulence model ttoward oceanic applications. Ocean Model. 5, 195–218. Van Kester, J.A.T.M., 1994. Validatie Delft3D voor menglaag proef. Technical Report. Waterloopkundig laboratorium. In Dutch Warner, J.C., Sherwood, C.R., Arango, H.G., Signell, R.P., 2005. Performance of four turbulence closure models implemented using a generic length scale method. Ocean Model. 8, 81–113. Wijesekera, H.W., Allen, J.S., Newberger, P.A., 2003. Modeling study of turbulent mixing over the continental shelf: comparison of turbulent closure schemes. J. Geophys. Res. 108, 3103. Wilcox, D.C., 1988. Reassessment of the scale-determining equation for advanced turbulence models. AIAA J. 26, 1299–1310. Willis, G.E., Deardorff, J.W., 1974. A laboratory model of the unstable planetary boundary layer. J. Atmosph. Sci. 31, 1297–1307. Zeierman, S., Wolfshtein, M., 1986. Turbulent time scale for turbulent-flow calculations. AIAA J. 24, 1606–1610.