A Reasoned Comparison between Two Hydrodynamic Models: Delft3D-Flow and ROMS (Regional Oceanic Modelling System)

Abstract

Useful information, such as water levels, currents, salinity and temperature dynamics in water bodies, are obtained through numerical models in order to pursue scientific research or consultancy. Model validation dates back long ago, since such models started to be developed in the 1960s. Despite their usefulness and reliability in complex situations, some issues related to well-known benchmarks are still present. This work aims to analyse in detail the behaviour of the velocity profile, vertical eddy viscosity and tangential stresses at the bed in two cases of free surface flows; namely: uniform flow in an inclined rectangular channel and a wind-induced circulation in a closed basin. Computational results strongly depend on the turbulence closure model employed and a reasoned comparison is necessary to highlight possible improvements of these models. The strong differences that arise are deeply analysed in this work.

1. Introduction

Environmental Science and Engineering dealing with riverine and coastal areas [1] are in need of reliable hydrodynamic models [2,3,4]. Such numerical models are of fundamental importance in order to provide some vital information. Currents, water levels, turbulence intensity and bottom stresses are only some of the variables required in order to carry out technical and scientific analyses that meet the stringent requirements that either public stakeholders or the scientific community require in order to accept the results of a study [5,6,7].

Several open-source and community based suites are now available in order to carry out such analyses. Among the vast options it is worth mentioning Delft3D-Flow [8], regional oceanic modelling system (ROMS), Princeton ocean model (POM) [9], shallow water hydrodynamic finite element model (SHYFEM) [10] and nucleus for European modelling of the ocean (NEMO) [11]. In particular, the present study is focused on the analysis of the performance of two models of the latter list, i.e., Delft3D-Flow and ROMS, them being two widely used models, particularly the former in the engineering community, and the latter in oceanography. However, several conclusions drawn in the present work might be applied to models based on similar approaches, especially regarding the turbulence closures.

During the last few decades, great efforts have been spent to provide such numerical models with turbulence closures [12] and grid generators [13] that allow for customising the numerical models for the various needs.

Both models have been extensively used in ocean and coastal areas. A vast literature is available on the implementation, validation and use of such numerical models. As far as Delft3D is concerned, it is possible to cite [14,15] as seminal works that established the reliability of this code. Applications of Delft3D can be found, for example, in [16], wherein tidal dynamics in a schematised mangrove creek catchment were modelled, and in [17], wherein water renewal was studied inside a Spanish estuary. Several other implementation can be found in literature; for example, in the open sea [18] or in estuaries [19]. Analogously, ROMS was effectively validated, for example, in [20], wherein the three dimensional terrain following grid was applied; in [21], where turbulence models were developed; and [22], wherein a high-resolution numerical model of the circulation around the Corsica Channel was implemented. Other applications of ROMS can be found, for example, in [23], in which the numerical simulations of the circulation in the Hudson River estuary are compared with an extensive set of time series and spatially resolved measurements, and in [24], in which the model was implemented to evaluate the near-shore circulation in Monterey Bay, California.

This work aims at critically comparing the results of Delft3D-Flow and ROMS over two well-known benchmarks. The first one is an open channel flow; its velocity profiles are carefully compared with the logarithmic theoretical one. Such a case can well model a tidal channel or an estuarine area. Theoretical bed shear stresses and model outputs are compared evaluating the relative error in order to assess the reliability of turbulence models adopted. The second benchmark considers the circulation inside a wind driven closed basin whose velocity profiles are analogous to Couette–Poiseuille turbulent profiles. Such velocity profiles are critically analysed considering the ratio between surface and bottom stresses and depth-averaged velocity.

The paper proceeds with Section 2 where the numerical models are introduced and the benchmarks described. Then, Section 3 describes the results obtained by comparing Delft3D-Flow, ROMS and their turbulence models over the chosen benchmarks. Eventually, the conclusions are drawn in Section 4.

2. Material and Methods

2.1. Model Description

In the present section we introduce the main characteristics of the modules under study, starting from the momentum equations, boundary conditions and turbulence closures. Owing to the strong similarities between Delft3D and ROMS, some of the following parts will be valid for both models. On the contrary, the differences in the formulations will be clearly stated and commented.

2.1.1. Primitive Equations

Both Delft3D and ROMS are based on the solution of the Reynolds averaged Navier Stokes equations (RANS) [25], which in their hydrostatic primitive form read:

$$\nabla ·\mathit{U}=0$$

(1)

$$\frac{\partial \mathit{U}}{\partial t}+(\mathit{U}·\nabla )\mathit{U}=-\frac{1}{{\rho }_0}\nabla P+\mathit{g}+\nu {\nabla }^2\mathit{U}+\nabla ·{\mathit{T}}^{\mathit{R}},$$

(2)

where $\mathit{U}=(U,V,W)$ stands for the mean velocity vector, P is the pressure, ${\rho }_0$ is the reference density (kept constant in this work), $\mathit{g}$ is the gravity acceleration vector $\mathit{g}=(0,0,-g)$, $\nu$ is the kinematic viscosity and ${\mathit{T}}^{\mathit{R}}$ is the Reynolds stress tensor defined as follows:

$${\mathit{T}}^{\mathit{R}}=-{\rho }_0\left(\begin{array}{ccc}\overline{uu}& \overline{uv}& \overline{uw}\\ \overline{vu}& \overline{vv}& \overline{vw}\\ \overline{wu}& \overline{wv}& \overline{ww}\end{array}\right),$$

(3)

where u, v and w are the components of the turbulent fluctuating velocity vector $\mathit{u}$. Reynolds stress tensor is modelled adopting the Boussinesq approach:

$$\overline{uw}=-{\nu }_{tv}\frac{\partial U}{\partial z},\overline{vw}=-{\nu }_{tv}\frac{\partial V}{\partial z},\overline{uv}=-{\nu }_{th}\left(\frac{\partial U}{\partial y}+\frac{\partial V}{\partial x}\right)$$

(4)

where ${\nu }_{tv}$ is the vertical eddy viscosity and ${\nu }_{th}$ is the horizontal eddy viscosity. In Delft3D and ROMS, ${\nu }_{th}$ can be set as a function of the grid resolutions and ${\nu }_{tv}$ can be calculated by the turbulence model adopted.

2.1.2. Boundary Conditions

Boundary conditions are needed to solve Equations (1) and (2). These conditions are set at the free surface and at the bottom and involve both velocities and tangential stresses. Both ROMS and Delft3D apply the vertical boundary conditions at the bottom reference level $z_{ob}$ where velocity is zero. Such a height is evaluated as $z_{ob}=z_r/30$, where $z_r$ is the Nikuradse roughness length scale. Both the numerical models use the reference level to compute the drag coefficient $C_{3D}$. In ROMS, it is defined as

$$C_{{3D}_R}={\left(\frac{\kappa }{log\left(\frac{\Delta z_b}{2z_{ob}}\right)}\right)}^2$$

(5)

which is slightly different from the one of Delft3D

$$C_{{3D}_D}={\left(\frac{\kappa }{log(1+\frac{\Delta z_b}{2z_{ob}})}\right)}^2.$$

(6)

Such a difference is present since ROMS was developed for ocean circulations characterised by very deep bathymetry, whereas Delft3D was initially intended for shallow flows of coastal areas. The former case is, therefore, insensible to the summation present in Equation (6), since for very deep flows $1+\frac{\Delta z_b}{2z_{ob}}\approx \frac{\Delta z_b}{2z_{ob}}$.

The drag coefficient $C_{3D}$ is the constant of proportionality between the bed velocity ${\mathit{U}}_{\mathit{b}}=(U_b,V_b)$ and the bottom tangential stress ${\mathit{\tau }}_{\mathit{b}}=({\tau }_b^x,{\tau }_b^y)$ using the typical quadratic closure:

$$({\tau }_b^x,{\tau }_b^y)=({\rho }_0{C}_{3D}{U}_b\sqrt{U_b^2+V_b^2},{\rho }_0{C}_{3D}{V}_b\sqrt{U_b^2+V_b^2}).$$

(7)

The models use ${\mathit{\tau }}_{\mathit{b}}$ to compute the bed friction velocity vector ${\mathit{u}}_{\mathit{b}}^{*}=(u_b^{*},v_b^{*})$ through the relation:

$${\mathit{u}}_{\mathit{b}}^{*}=\sqrt{\frac{{\mathit{\tau }}_{\mathit{b}}}{{\rho }_0}}.$$

(8)

Bottom boundary conditions for velocity ${\mathit{U}}_{\mathit{b}}=(U_b,V_b)$ are computed by assuming a logarithmic velocity profile in the lowest computational cell:

$${\mathit{U}}_{\mathit{b}}=\frac{{\mathit{u}}_{\mathit{b}}^{*}}{\kappa }log\left(1+\frac{\Delta z_b}{2z_{ob}}\right)$$

(9)

where $\kappa$ is the von Kármán’s constant equal to $0.41$, $\Delta z_b$ is the distance to the computational grid point closest to the bed and $z_{ob}$ is the reference level where velocity goes to zero. A typical drag coefficient $C_{3D}$ is adopted to relate tangential stresses at the bed with velocities at the first bottom layer.

Boundary conditions at the bottom and at the free surface involve the vertical eddy viscosity ${\nu }_v$ and relate the derivatives of the plane components of velocity with the respective tangential stresses. Vertical velocity is instead obtained from continuity. At the free surface, i.e., for $z=\zeta (x,y,t)$, being $\zeta (x,y,t)$ the free surface elevation above reference plane, the following conditions hold:

$$\frac{\partial U}{\partial z}=\frac{{\tau }_s^x(x,y,t)}{{\rho }_0{\nu }_v},\frac{\partial V}{\partial z}=\frac{{\tau }_s^y(x,y,t)}{{\rho }_0{\nu }_v},W=\frac{\partial \zeta }{\partial t}.$$

(10)

At the bottom, for $z=-h(x,y)$, where $h(x,y)$ is the depth below the reference plane, the conditions are:

$$\frac{\partial U}{\partial z}=\frac{{\tau }_b^x(x,y,t)}{{\rho }_0{\nu }_v},\frac{\partial V}{\partial z}=\frac{{\tau }_b^y(x,y,t)}{{\rho }_0{\nu }_v},W=(U,V)·\nabla h.$$

(11)

2.1.3. Turbulence Closures in ROMS

In ROMS the classical two-equation models are present through the generic length scale (GLS) approach [26]. The GLS is a two-equation model that takes advantage of the similarities in other two-equation formulations. The first equation in GLS model is the standard equation for the turbulent kinetic energy k, whereas the second equation is for a generic variable $\psi$. Hence, the equations are:

$$\frac{\partial k}{\partial t}+U_i\frac{\partial k}{\partial x_i}=\frac{\partial }{\partial z}\left(\frac{{\nu }_{tv}}{{\sigma }_k}\frac{\partial k}{\partial z}\right)+P+B-ϵ$$

(12)

$$\frac{\partial \psi }{\partial t}+U_i\frac{\partial \psi }{\partial x_i}=\frac{\partial }{\partial z}\left(\frac{{\nu }_{tv}}{{\sigma }_{\psi }}\frac{\partial \psi }{\partial z}\right)+\frac{\psi }{k}(c_{1}P+c_{3}B-c_2ϵ),$$

(13)

where ${\sigma }_k$ and ${\sigma }_{\psi }$ are the turbulence Schmidt numbers respectively for k and $\psi$ defined as the ratio between the vertical eddy viscosity and the vertical eddy diffusivity. $c_1$, $c_3$ and $c_2$ are coefficients chosen to be consistent with von Kármán’s constant and with experimental observations. The generic parameter $\psi$ can be related to the other turbulent quantities as $\omega$ (vorticity), $ϵ$ and $kl$ through the following relation, where l is the length scale of characteristic vortices:

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

(14)

P represents production by shear:

$$P=-\overline{uw}\frac{\partial U}{\partial z}-\overline{vw}\frac{\partial V}{\partial z}={\nu }_{tv}{M}^2,$$

(15)

where

$$M^2={\left(\frac{\partial U}{\partial z}\right)}^2+{\left(\frac{\partial V}{\partial z}\right)}^2.$$

(16)

B is the buoyancy term:

$$B=-\frac{g}{{\rho }_0}\overline{\rho w}={\nu }_{th}{N}^2,$$

(17)

where N is the Brunt–Vaisala frequency:

$$N^2=-\frac{g}{{\rho }_0}\frac{\partial \rho }{\partial z}$$

(18)

The dissipation rate $ϵ$ is modelled according to this dimensional relation:

$$ϵ={\left(c_{\mu }^0\right)}^{3+p/n}{k}^{3/2+m/n}{\psi }^{-1/n},$$

(19)

where $c_{\mu }^0$ is the stability coefficient based on experimental data for unstratified channel flow with a log-layer solution.

The explicit formulation for three closures, $k-kl$, $k-\omega$ and $k-ϵ$, are described with specific combinations of coefficients as reported in

Table 1. Coefficient combination used to set turbulence closures through GLS approach.

Closure	p	m	n
$k-kl$/$MY25$	0.0	1.0	1.0
$k-\omega$	$-1.0$	0.5	$-1.0$
$k-ϵ$	3.0	1.5	$-1.0$

.

ROMS also provides the Mellor–Yamada level 2.5 scheme (MY25) [27], which is commonly used in oceanographic models. It can be selected with this parameter combination: $p=0$, $m=1$ and $n=1$ (same as $k-kl$). It differs from the classical $k-kl$ in the vertical eddy diffusivity coefficient $D_Q$ for transport of k and $kl$, and in the use of specific wall functions $F_{wall}$. The equations for k and $kl$ in MY25 are:

$$\frac{\partial k}{\partial t}+U_i\frac{\partial k}{\partial x_i}=\frac{\partial }{\partial z}\left(D_Q\frac{\partial k}{\partial z}\right)+P+B-ϵ$$

(20)

$$\frac{\partial \left(kl\right)}{\partial t}+U_i\frac{\partial \left(kl\right)}{\partial x_i}=\frac{\partial }{\partial z}\left(D_Q\frac{\partial \left(kl\right)}{\partial z}\right)+l(c_{1}P+c_{3}B-c_2ϵF_{wall}),$$

(21)

where in the standard implementation $D_Q=\sqrt{2k}l{S}_q$ with $S_q=0.2$ and $F_{wall}$ is the defined as follows:

$$F_{wall}=\left(1+E_2{\left(\frac{l}{k}\frac{d_b+d_s}{d_b{d}_s}\right)}^2\right),$$

(22)

where $E_2=1.33$. The parameters $d_b$ and $d_s$ are the distances to the bottom and free surface, respectively. The other parameters are described above.

The boundary conditions for k are applied in flux form assuming local, steady state, no-gradient conditions with the equilibrium layer hypothesis, that provide $P=ϵ$, to yield no-flux conditions respectively for $z=\zeta (x,y,t)$ at the free surface and for $z=-h(x,y)$ at the bed [12]:

$${\left(\frac{{\nu }_{tv}}{{\sigma }_k}\frac{\partial k}{\partial z}\right)}_b=0,{\left(\frac{{\nu }_{tv}}{{\sigma }_k}\frac{\partial k}{\partial z}\right)}_s=0.$$

(23)

The boundary conditions for the generic length scale $\psi$ are applied in flux form, again, respectively, for $z=\zeta (x,y,t)$ and for $z=-h(x,y)$ as follows:

$${\left(\frac{{\nu }_{tv}}{{\sigma }_{\psi }}\frac{\partial \psi }{\partial z}\right)}_b=-n\frac{{\nu }_{tv}}{{\sigma }_{\psi }}{\left(c_{\mu }^0\right)}^p{k}^{m}kn{z}_{ob}^{n-1},{\left(\frac{{\nu }_{tv}}{{\sigma }_{\psi }}\frac{\partial \psi }{\partial z}\right)}_s=-n\frac{{\nu }_{tv}}{{\sigma }_{\psi }}{\left(c_{\mu }^0\right)}^p{k}^{m}kn{z}_{os}^{n-1}.$$

(24)

2.1.4. Turbulence Closures in Delft3D

Delft3D provides three possible models of turbulence closures. The algebraic closure is a zero-equation model and assumes a logarithmic velocity profile. This leads to a linear relation between the turbulent kinetic energy at the bed and the turbulent kinetic energy at the free surface:

$$k=\frac{1}{\sqrt{c_{\mu }^0}}\left({\left(u_b^{*}\right)}^2\left(1-\frac{z+H}{H}\right)+{\left(u_s^{*}\right)}^2\left(\frac{z+H}{H}\right)\right),$$

(25)

where $H=h(x,y)+\zeta (x,y,t)$ is the total depth; and $c_{\mu }^0$ is the stability coefficient which was discussed earlier and is assumed equal to $0.09$. $u_b^{*}$ and $u_s^{*}$ are the friction velocities at the bottom and at the free surface respectively. In this case the vertical eddy viscosity ${\nu }_{tv}$ is defined as follows:

$${\nu }_{tv}={\left(c_{\mu }^0\right)}^{1/4}L\sqrt{k},$$

(26)

where L the mixing length defined by Prandtl [28] as follows

$$L=\kappa (z+H)\sqrt{1-\frac{z+H}{H}}.$$

(27)

Then, Delft3D implements the $k-L$ closure that is a one-equation model, based on the same definition of L (Equation (27)) and on the transport equation for the turbulent kinetic energy k (Equation (12)). The third model implemented is the $k-ϵ$ closure. This is the usual two-equation model [29] based on the solution of the transport equation for the turbulent kinetic energy k and of the transport equation of the dissipation rate $ϵ$

$$\frac{\partial k}{\partial t}+U_i\frac{\partial k}{\partial x_i}=\frac{\partial }{\partial z}\left(\frac{{\nu }_{tv}}{{\sigma }_k}\frac{\partial k}{\partial z}\right)+P+B-ϵ$$

(28)

$$\frac{\partial ϵ}{\partial t}+U_i\frac{\partial ϵ}{\partial x_i}=\frac{\partial }{\partial z}\left(\frac{{\nu }_{tv}}{{\sigma }_{ϵ}}\frac{\partial ϵ}{\partial z}\right)+\frac{ϵ}{k}\left(c_{1}P+c_{3}B-c_2ϵ\right).$$

(29)

In this case the vertical eddy viscosity ${\nu }_{tv}$ is defined as follows:

$${\nu }_{tv}=c_{\mu }^0\frac{k^2}{ϵ}.$$

(30)

The boundary conditions for k and $ϵ$ are applied in Dirichlet form assuming that the hypothesis of equilibrium $P=ϵ$ is valid, for which [30]:

$$k_b=\frac{{u_b^{*}}^2}{\sqrt{c_{\mu }^0}},k_s=\frac{{u_s^{*}}^2}{\sqrt{c_{\mu }^0}},{ϵ}_b=\frac{{u_b^{*}}^3}{k{z}_{ob}},{ϵ}_s=\frac{{u_s^{*}}^3}{\frac{1}{2}k\Delta z_s}.$$

(31)

2.2. Description of the Test Cases

We tested both Delft3D and ROMS against two different well-known benchmarks. The first one is a fundamental study case where the solution of the vertical profile of the longitudinal velocity and the relative shear stress is well known. This problem is characterised by one turbulence scale that is the friction velocity at the bottom $u_b^{*}$. By modelling this benchmark we aimed at testing the different turbulence closures mentioned before, since turbulence properties are well-known. Then, it was possible to detect the influence of the different types of bottom boundary conditions relative to the velocities and to the turbulent quantities. The second one is a wind driven circulation in a closed basin, a well studied problem (see for example [31,32]) characterised by two turbulence scales that are the friction velocities at the bottom $u_b^{*}$ and at the free surface $u_s^{*}$. In addition, this benchmark is used to test and compare the different turbulence closures. However, in this case the turbulence dynamics are not simple to describe. Indeed, thanks to the particular geometry, in a depth-varying cross-section vertical profiles of velocities are different moving along the transverse direction, reproducing, in part, some of the possible solutions that can be generated from a flow between two moving plates (Couette–Poiseuille type).

The analogy between the wind driven circulation in closed basins and the turbulent Couette–Poiseuille flows has been already discussed in [31], and a detail analysis of the characteristics of these flows can be found in the cited study. In the present context, the wind-driven circulation can be a useful benchmark for the model at hand, owing to the variety of vertical profiles of longitudinal velocity that can be encountered in a simple geometry, such as a closed basin.

2.2.1. Open Channel Flow

The uniform open channel flow is a free-surface channel flow characterised by a constant slope and a constant discharge. The vertical velocity profile is the well-known logarithmic profile [33]. Figure 1 shows a sketch of the free-surface channel flow we modelled. The bottom of the channel is depicted in red while the surface in blue. In steady conditions these surfaces should be parallel to each other. The profiles of velocity, tangential stresses and turbulent eddy viscosity are shown as a function of the vertical axis z.

We define a non-dimensional conductance coefficient C as the ratio between the depth-averaged velocity $\overline{U}$ and the bed friction velocity $u_b^{*}$, which can be related to the flow depth and the sediment diameter $d_s$ through the Keulegan [33] equation for fully rough turbulent flows:

$$C=\frac{\overline{U}}{u_b^{*}}=\frac{1}{\kappa }log\left(\frac{11.09R_i}{2.5d_s}\right),$$

(32)

where $R_i$ is the hydraulic radius, i.e., $R_i=\Omega /(W+2H)$, W being the width of the channel and H the water depth, and $d_s$ is the average diameter of the sediments. In this expression $2.5d_s$ is considered as the Nikuradse roughness length scale $z_r$. The choice of this expression is justified by the fact that it derives directly from the rough wall law and it is adjusted for rectangular riverbeds. The expected water depth can be evaluated thanks to the uniformity of the flow by solving iteratively the following equation:

$$Q=\Omega C\sqrt{g{R}_i{i}_f},$$

(33)

where Q is the flow, $\Omega$ the cross section of the channel, g the acceleration of gravity and $i_f$ the slope.

We simulate such a flow over a 1000 m long and 50 m wide channel where a constant discharge of 300 m³/s flows. The grid is made of constant square cells of 5 m sides. The velocity profile and all the results are taken at the centre of the channel in order to discard possible effects of the boundaries. The channel presents a constant slope of 0.1%. Delft3D presents two open boundaries conditions of total discharge at both ends of the channel. The total discharge is modelled in order to increase linearly from a zero value to 300 m³/s in one hour. After that, steady conditions of flows are reached. Analogously for ROMS, conditions of a constant current parallel to the bed are chosen since conditions of total discharge are not provided as a default option in the code. The current velocity is chosen as 2.47 m/s in order to allow a constant flow of 300 m³/s in steady conditions and to carry out a fruitful comparison with Delft3D. The simulation parameters are summarised in

Table 2. Parameter for the open channel flow test case.

Model Parameter	Variable	Value
Length, Width	L, W	1000, 50 m
Bed Slope	$i_f$	0.001
Number of grid spacing	$Lm$, $Mm$, $Nm$	200, 10, 40
Bottom roughness	$z_r$	0.03 m
Time step	$\Delta t$	1 s
Total simulation time	$t_{end}$	86,400 s
Uniform water depth	H	2.428 m
Hydraulic Radius	$R_i$	2.213 m
Flow, Velocity	Q, U	300 m3/s, 2.47 m/s

.

Results of the benchmark are analysed in Section 3 and compared with the theoretical solution expected, which is the usual logarithmic profile:

$$\frac{U}{u_b^{*}}=\frac{1}{\kappa }log\left(\frac{z}{z_r}\right)+8.5.$$

(34)

The velocity U goes to zero when $z=z_r/30$. Consequently, tangential stresses and vertical eddy viscosity were expected to have a linear and a parabolic profile, respectively [34].

2.2.2. Wind Driven Trapezoidal Closed Basin

The wind driven circulation in a closed basin consists of a squared closed basin with a trapezoidal cross-section where the motion is induced only by a steady wind blowing in the x-direction.

Figure 2 shows the geometry of the basin in blue. The basin is 10 km × 10 km with a maximum depth of 15 m in the flat central region of 5 km width. Lateral banks have a water depth that linearly decreases from 5 m to 15 m along 2.5 km. The numerical grid is made of square cells of 100 m × 100 m side and the vertical is discretised with 40 equidistant layers. Layer b) of Figure 2 shows vectors representing the horizontal velocity field induced by wind, and the free surface is depicted in red in Layer c). The simulation parameters are summarised in

Table 3. Model parameters for the wind driven trapezoidal closed basin test case.

Model Parameter	Variable	Value
Length, Width	L, W	10,000, 10,000 m
Depth	H	from $-5$ m to $-15$ m
Number of grid spacing	$Lm$, $Mm$, $Nm$	100, 100, 40
Bottom roughness	$z_r$	0.03 m
Time step	$\Delta t$	1 s
Total simulation time	$t_{end}$	86,400 s
Wind Driven Surface Stress	${\tau }_x^s$	0.325 N/m2

.

In a depth-varying cross-section, vertical profiles of velocities are quite different moving along the transverse direction. Such velocity profiles were studied in detail by [31,32,35,36]. In the shallow regions close to the lateral banks, flow velocity is always positive (i.e., in the wind direction), while a reverse flow occurs in the central region. Zero-mean flow only appears at two symmetrical locations in the cross-section, corresponding to the centre of the two counter-rotating cells of the depth-averaged circulation pattern, visible in Figure 2.

A condition that can be safely assumed to be valid in any region of the basin is a linear distribution of the shear stress along the vertical, with the only exception of the area close to the leeward and windward banks, where the flow is strongly three-dimensional. The shear stress can be thought to be due to the superposition of two contributions: a constant one due to the wind stress acting on the free surface and a linear one that can be related to the free surface slope (the wind set-up).

Since the set-up is fairly constant along the transverse coordinate, so is the vertical gradient of the shear stresses. Then, as sketched in Figure 3a, the vertical distribution of shear stresses is always positive (no reversal) in the region close to the lateral banks. This profile is associated with point A depicted in Figure 2. The velocity profile is also shown in red and can be associated to that of a Couette–Poiseuille flow. A negative value of the bed shear stress (associated to flow reversal) is expected in the central region, Figure 2c and point C. The flow reversal is evinced by the velocity profile in red: the maximum velocity is located at the surface and shows a positive direction. At the middle of the water column there is the flow reversal. The peculiar situation of zero bed shear stress, point B of Figure 2b, marks the boundary between the two regions. In this case, the velocity profile tends to the bottom with a zero-gradient.

The importance of this benchmark relies on the possibility of spanning a great number of depth averaged velocities and tangential stresses at the same time, moving along the central cross-section. As a result, the comparison between the turbulence models will be much more accurate than a single point-to-point difference since a wider range of results will be taken into consideration.

3. Results and Discussion

Both Delft3D and ROMS have been tested over the two benchmarks we introduced in Section 2. We report the results with particular attention to the behaviour of several turbulence models at the bottom layers.

3.1. Open Channel Flow

The case of uniform open channel flow was carried out assuming a bed roughness $z_r$ = 0.03 m. Employing Equations (32) and (33) it is possible to compute the water depth H and the friction velocity at the bottom $u_b^{*}$ for a constant flow of 300 m³/s, respectively, to be equal to 2.428 m and 0.147 m/s. The conditions are those of a rough surface since the rough Reynolds number $R{e}_{*}=u_b^{*}{z}_r/\nu =4420$ is much greater than 100 [37]. Figure 4 shows the comparison among the non dimensional vertical profiles of the longitudinal velocity, eddy viscosity and turbulent kinetic energy obtained with the available turbulence closures of Delft3D and ROMS. In particular, Delft3D adopts $k-ϵ$, $k-L$ and algebraic closure models. Conversely, ROMS adopts $k-ϵ$, $k-\omega$, $k-kl$ and $MY25$ turbulence models. Figure 4a,b show on a semi-log plot, the vertical profile of the longitudinal velocity made non dimensional with the bed friction velocity, together with the theoretical log law of Equation (34). Each profile obtained from the turbulence models is made non-dimensional either using the model bed shear stress outputs since $u_b^{*}=\sqrt{\frac{{\tau }_b}{{\rho }_o}}$ or the theoretical $u_b^{*}$ = 0.147 m/s. Such results are shown in panels (a) and (b) respectively. The choice of normalising the velocity profile with theoretical and output friction velocities is due to the fact that the curves should collapse over the rough wall law. In particular, comparing panels (a) and (b), it is possible to appreciate that the spread of the profiles is greater for high values of $z/z_r$ adopting the output friction velocity. Results reported in Figure 4b show an increasing overlapping. Therefore, it is possible to argue that, although the profiles follow the expected slope in the semi-log plot, the output bottom stress from the model does not conform to the profiles obtained. In addition, our attention is drawn from the points at the bed since they do not follow the rough wall law. Such an issue is independent from the normalisation adopted since both panels show the same shortcomings. Boundary conditions at bottom do influence the bed friction velocity, and consequently, the velocity at the bottom layer. It is the evaluation of the bed friction velocity that must be modified adequately in order to produce a velocity at the bottom that is in agreement with the rough wall law. Another issue consists in the differences between the bed friction velocities generated from one turbulence model to another. As a result, whether we normalise with the theoretical bed friction velocity at the bottom or not, the curves tend to disperse around the wall law. This is the sign of possible inconsistencies between the various turbulence models adopted by Delft3D and ROMS.

Figure 4c reports the vertical eddy viscosity normalised with theoretical friction velocity and the total depth for the several turbulence models adopted. The analytic solution was obtained assuming homogeneous turbulence, velocities V and W equal to zero and a logarithmic profile of velocity. As a result, the vertical eddy viscosity follows the well known parabolic profile that reads ([30]):

$${\nu }_{tv}=\kappa u_b^{*}z\left(1-\frac{z}{H}\right).$$

(35)

It is possible to see that there are high discrepancies in the central region, especially for the two-equation models, whereas the algebraic and the $k-L$ model of Delft3D seem to show a better agreement with the expected one. Note that the best performance among the two equation models is achieved by the $k-\omega$ formulation of ROMS. It is interesting to note that the numerical profiles show significant differences not only in the central part of the channel, but also in the vicinity of the bottom and the free surface. Note that the slope of the turbulent viscosity profile near the walls also influences the calculation of the tangential stresses at the boundaries, as discussed at the end of this section. In addition, the large differences along the depth may also have consequences in applications related to mass transport. In fact, models such as Delft3D and ROMS are able to solve water quality problems through the solution of the advection–diffusion equation, whose mass transport estimates strongly depend on the diffusive and dispersive transport coefficients. As far as the coefficients of turbulent diffusion along the vertical are concerned, these models assume that the mass transport and the momentum transport are substantially the same, based on the so-called Reynolds analogy, which implies that the number of turbulent Schmidt is actually unitary [38]. This implies that an inaccurate estimation of the turbulent viscosity along the vertical leads inevitably to an inaccurate estimation of the vertical turbulent diffusion coefficients.

Panel (d) shows the turbulent kinetic energy k normalised with the squared theoretical friction velocity $u_b^{*2}$. A physical inconsistency is here quite evident. At the bed, k should be zero since velocities and turbulent fluctuations should disappear due to the no-slip conditions. The condition according to which $k/u_b^{*2}$ is equal to $1/\sqrt{c_{\mu }^0}$ is valid in correspondence of the equilibrium layer that is set in $z=30\nu /u_b^{*}$ starting from the bed [30]. The consequence of this physical inconsistency relies in the fact that the production term P and the dissipation rate $ϵ$ are equal starting from $z=0$, instead of the right depth, which is $z=30\nu /u_b^{*}$.

Figure 5 shows the comparison among the vertical profiles of shear stress obtained with Delft3D and ROMS with the different turbulence closures. Note that both models do not provide the vertical profile of shear stress, but, instead, they compute only the bottom shear stress only following Equation (7), and this is the value reported in the plots as “model output.” The vertical profiles of shear stress have been computed starting from the vertical profiles of longitudinal velocity and the corresponding eddy viscosity. In the same plot of Figure 5 the theoretical distribution of shear stress for a uniform open channel flow is also reported. By inspecting Figure 5, it appears quite clearly that only few turbulent closures for both Delft3D and ROMS are able to correctly predict the vertical profile of the shear stress. In particular, the best comparison between the numerical solution against the theoretical linear profiles was obtained with the $k-\omega$ closure for ROMS and with the $k-L$ model and algebraic closures for Delft3D. The good performance of the $k-\omega$ model of ROMS was expected owing to the satisfactory comparison of the vertical velocity profile in Figure 4a,b. It is not surprising that the simplest turbulent models are most often not consistent with the whole vertical linear distribution of shear stress (compare each profile with the symbol representing the model output). This discrepancy is directly associated to the already discussed difference in the velocity derivative close to the bottom with respect to the theoretical log law of the wall; see Figure 4a,b. The effect of an incorrect vertical profile of the longitudinal velocity close to the bottom layers is reflected also into the non linear profile for the shear stress close to the bottom; see, for example, the results obtained with the ROMS $k-ϵ$. In

Table 4. Bottom tangential stresses ${\tau }_b^x$ for Delft3D and ROMS at the varying of turbulence models. The relative error refers to theoretical value equal to 21.713 N/m². More than three significant digits were employed in the calculations.

ROMS	${\mathit{\tau }}_{\mathit{b}}^{\mathit{x}}$	Relative Error (%)	Delft3D	${\mathit{\tau }}_{\mathit{b}}^{\mathit{x}}$	Relative Error (%)
$k-ϵ$	22.095	1.759	$k-ϵ$	24.471	12.702
$k-\omega$	22.411	3.215	$k-L$	21.537	0.811
$k-kl$	20.450	5.817	$Algebraic$	21.951	1.096
$MY25$	16.988	21.761

the error in the evaluation of the bottom shear stress is reported for all models; it was calculated as the relative error between the theoretical value and the model output. The relative error can reach values greater than $10\%$ depending on the model adopted. Note the $k-ϵ$ Delft3D, one of the most used models for geophysical applications, produces a difference of about $12\%$. It is worth noting that the error could be even larger if it is computed using the vertical profiles instead of the bottom shear stress provided by the model.

Finally, we stress how the uncertainties identified regarding the estimation of the bottom shear stress could have significant repercussions in many applications. Indeed, the models analysed in this study are often used for evaluations of solid transport in coastal, lagoon and estuarine environments. An inaccurate assessment of the bottom shear stress causes an incorrect estimation of the solid transport, and, ultimately, errors in the estimation of the morphodynamics of these environments. We will see in the next section that similar uncertainties are evident in the case of wind driven flows as well.

3.2. Wind Driven Circulation in a Closed Basin

The circulation in a trapezoidal closed basin was considered in order to provide an assessment of the behaviour of the numerical models corresponding to different ratios between surface and bottom stresses, as depicted in Figure 3. Particular attention is drawn by velocity profiles in points A, B and C of Figure 2 since it is possible to cast an analogy with Couette–Poiseuille turbulent flows and the present circulation.

Figure 6 shows the vertical profiles of the longitudinal velocity computed in such points. Panel (a) shows the behaviour of the profiles along the lateral banks, i.e., point A. Panel (b) shows that in point B velocities tend to zero with significantly different profiles. This happens since zero tangential stresses are reached in different areas depending on the turbulence closure model adopted. This is confirmed by Figure 7 where tangential stresses are shown for the $k-ϵ$ closures for both models (other turbulence models are omitted for the sake of brevity). As far as panel (c) of Figure 6 is concerned, it is clear that algebraic and $k-L$ model of Delft3D cannot reproduce correctly the flow inversion in point C. Indeed, this panel shows that only two-equation turbulence models can follow the strong gradients expected in this region, providing a positive velocity at the surface and a reversal in the central region.

Despite similarities in the magnitude of the bed shear stresses depicted in Figure 7, the patterns show differences that can be of importance for sediment transport. The activation of bed transport occurs when shear stresses overcome a threshold. Activation of bed transport, therefore, can differ from point to point. Therefore, it is important to observe that the different patterns of bed shear stress determine different positions of the points in which the value of the bed stress is zero. Figure 8 shows the x-component velocity profiles for every turbulence closure taken in exam in accordance of zero shear-stress at the bottom. The position in the basin of these points changes as a function of the turbulence closure model. It is evident that the behaviour of the algebraic and $k-L$ closures does not manage to capture strong gradients in the velocity profiles. Figure 8 also reports the analytical solution obtained by [32]. They introduced a bi-logarithmic model to describe the velocity profile:

$$u\left(z\right)=A{u}_s^{*}ln\left[1+\frac{z}{z_{os}}\right]+B{u}_b^{*}ln\left[1-\frac{z}{(z_{ob}+H)}\right]+C,$$

(36)

where $z_{os}$ and $z_{ob}$ are the reference levels at the free surface and at the bottom; H is the total depth; and $u_s^{*}$ and $u_b^{*}$ are the friction velocities, respectively, for the free surface and the bottom. The coefficients A, B and C are defined in function of $z_{os}$ and $z_{ob}$ as follows:

$$A=\frac{q_2}{(p_1{q}_2-q_1{p}_2)},B=-\frac{q_1}{(p_1{q}_2-q_1{p}_2)},C=0$$

(37)

$$p_1=\lambda \frac{z_{os}}{H},p_2=\lambda \frac{z_{os}}{z_{ob}}$$

(38)

$$q_1=\left(1+\frac{z_{os}}{H}\right)ln\left[1+\frac{H}{z_{os}}\right]-1,q_2=\left(\frac{z_{os}}{H}\right)ln\left[1+\frac{H}{z_{ob}}\right]-1.$$

(39)

The constant $\lambda$ represents the turbulence intensity and the values suggested are between $0.2$ and $0.5$. The mean value of $0.35$ is considered optimal and herein adopted. Attention should be drawn by this case since the assumptions behind the boundary conditions, i.e., equilibrium between production and dissipation of turbulent kinetic energy, tends to disappear. The conditions on the dissipation rate $ϵ$ start to be physically uncertain whenever the boundary shear stress tends to vanish. As shown by [31], whenever the bed shear stress tends to zero and surface and bottom friction velocities differ of one order of magnitude, there exists a behaviour called “erosion” of the logarithmic layer already observed in seminal studies, as in [39].

The results of the turbulence models are summarised in Figure 9 where the y-axis represents the ratio between the bed shear stress and surface shear stress ${\tau }_b^x/{\tau }_s^x$, whereas the x-axis represents the the ratio between the depth-averaged velocity and the surface friction velocity $U/u_s^{*}$ for the points of the central cross-section of the closed basin. Note that $u_s^{*}$ is constant for all the points of the basin and for all the simulations since the surface stress induced by wind is constant. Therefore, Figure 9a aims at highlighting the spread of results and the degree of reliability that each model can have. The spread of results is evidenced by the bunch of curves that cover a fairly wide range of values. Figure 9b is a zoom around the crossing of the x-axis by the curves. It is possible to observe that all models manage to reproduce a negative bottom stress while there exists a positive mean velocity calculated over the water column. The greater differences can be seen for high values of $U/u_s^{*}$ where the behaviour tends to be parabolic as expected by the typical quadratic relation between velocities and tensions expressed by Equation (7). Here, it is worth mentioning that this typical formula is widespread used in shallow water models that solve purely the planar components of velocity. Delft3D and ROMS could be used in this approach, too. In this case, the formula $\mathit{\tau }=C\mathit{U}\left|\mathit{U}\right|$ cannot reproduce the behaviour represented in the Figure 9b. This limitation can be overcome adopting the formulation proposed by [31].

4. Conclusions

A reasoned comparison between Delft3D-Flow and ROMS was carried out with the aim of evaluating the behaviour of such models in agreement with the turbulence models implemented, especially as far as at the bed velocities and stresses are concerned. Two benchmarks were considered: the first one was a uniform open-channel flow whose logarithmic velocity profile can be derived analytically, whereas the second one was a wind driven circulation in a closed basin induced by a steady wind. The former case was well suited to evaluate the agreement of the model with expected theoretical results, while the latter was taken into consideration since it allowed us to validate the performance of turbulence closures at the varying of the ratio between bed and surface stresses.

A degree of variability in the results can be expected, especially adopting closures usually implemented in different contexts. For example, $MY25$ model implemented in ROMS is also much used in atmospheric models such as WRF (weather research and forecasting model) [40] for stratified environments. However, concerning coastal and estuarine environments, such a model should not be used, as the first benchmark clearly demonstrates. Besides, this benchmark provides a useful indication on which model should be adopted in these areas. The profiles of longitudinal velocity obtained with algebraic and $k-L$ models of Delft3D follow the analytical solution even at lower layers close to the bed. On the other hand, $k-ϵ$ model of Delft3D is ill-behaved because does not manage to reproduce the theoretical result at the bottom. Therefore, the bed stresses computed by the models can be accepted only for algebraic and $k-L$ turbulence closures. As far as ROMS is concerned, $k-ϵ$ and $k-\omega$ models provide acceptable results in terms of bed stresses (see Table 4). Excluding the bottom layer, $k-\omega$ model follows the rough wall law up to the surface. Indeed, none of the models tested provided an acceptable velocity profile at the bottom.

The second benchmark is well-suited to evaluate the variability in the results, discarding possible influences of open boundaries. The variability is quite widespread, as depicted in Figure 9, and this must be kept in mind in case sediment transport should be evaluated. Besides, the behaviour of the models in this benchmark opposed that of the first one. Here, algebraic and $k-L$ models did not manage to reproduce the inversion of the flow in the central part of the basin, as depicted in Figure 6.

In conclusion, this work clarifies that in open channel flows, typical of riverine and estuarine areas, algebraic and one-equation models should be preferred to two-equation models. On the contrary, in the presence of strong gradients and flow reversals, two-equation models manage to provide much more realistic results despite the patterns of bed shear stresses not being in complete agreement.

Further studies will be devoted to develop better boundary conditions for two-equation models that can overcome the issues reported here. A possible solution would be to apply in $k-\omega$ models, the conditions proposed by [41] alongside a stretching of the $\sigma$-layers that emphasises the refining at the bottom layers. A useful stretching is reported in [31]. Besides, in two-dimensional (2D) shallow-water models the quadratic dependence of stresses from velocity through a conductance coefficient, $\mathit{\tau }=C\mathit{U}\left|\mathit{U}\right|$, should be used with care since it cannot reproduce the behaviour of Figure 9.