Arch Comput Methods Eng (2007) 14: 303–341

DOI 10.1007/s11831-007-9009-3

Depth Averaged Modelling of Turbulent Shallow Water Flow with

Wet-Dry Fronts
Luis Cea · Jerónimo Puertas ·
María-Elena Vázquez-Cendón

Published online: 9 August 2007

© CIMNE, Barcelona, Spain 2007

Abstract Depth averaged models are widely used in engi- Under these conditions the 3D Reynolds averaged Navier-
neering practice in order to model environmental flows in Stokes equations (3D-RANS) can be simplified in order to
river and coastal regions, as well as shallow flows in hy- obtain the depth averaged shallow water equations. Here, the
draulic structures. This paper deals with depth averaged tur- term shallow refers to a small ratio between the vertical and
bulence modelling. The most important and widely used horizontal length scales. Shallow flows appear in many en-
depth averaged turbulence models are reviewed and dis- gineering applications, mainly in river and coastal engineer-
cussed, and a depth averaged algebraic stress model is pre- ing, but also in certain hydraulic structures like open chan-
sented. A finite volume model for solving the depth aver- nels or sedimentation tanks, just to cite some examples.
aged shallow water equations coupled with several turbu- An important feature of free surface flows is that they are
lence models is described with special attention to the mod- unbounded in space, being the limits of the spatial domain
elling of wet-dry fronts. In order to asses the performance an unknown of the problem to be solved. Frequent problems
of the model, several flows are modelled and the numerical in which the limits of the fluid are unknown and unsteady in-
results are compared with experimental data. clude, among others, tidal flow in estuaries and flood propa-
gation in rivers. In those situations it is necessary to compute
a non-stationary wet-dry front, which is part of the solution
1 Introduction looked for. In addition to the former considerations, and due
to the large length scales involved in hydraulic engineering
This paper is focused on the numerical modelling of quasi- practise, the Reynolds number is usually large enough for
2D free surface turbulent flows. The two-dimensional char-
the flow to be fully turbulent. Hence, in most environmen-
acter of a free surface flow is usually enforced by a hori-
tal hydraulic engineering studies it is necessary to deal with
zontal length scale much larger than the vertical one, and
unsteady turbulent shallow water flow with wet-dry fronts,
by a velocity field quasi-homogeneous over the water depth.
and that is what this paper deals with.
The choice of a specific numerical model in hydraulic en-
L. Cea () · J. Puertas
gineering depends on the problem considered. The models
Departamento de Métodos Matemáticos y de Representación, most commonly used nowadays for the computation of shal-
Universidad de A Coruña, Campus de Elviña, 15071 A Coruña, low water flows solve either the 3D shallow water equations
Spain (3D-SWE) or the 2D shallow water equations (2D-SWE).
e-mail: [email protected]
Although its use is becoming more common due to the in-
J. Puertas
crease in computational power at the present time, the 3D-
e-mail: [email protected]
SWE are mainly used to model the flow in hydraulic struc-
M.-E. Vázquez-Cendón tures or in problems with a relatively simple geometry. In
Departamento de Matemática Aplicada, Universidad de Santiago the modelling of realistic coastal regions and estuaries the
de Compostela, Campus Sur, 15782 Santiago de Compostela,
Spain
computational cost of a 3D model becomes very expensive.
e-mail: [email protected] First, because the spatial domain is usually extensive and
304 L. Cea et al.

complex, such that a large numerical mesh is needed. Sec- a coastal estuary, in an open channel with a 90° bend and
ond, because the flow is usually unsteady, specially in en- in two designs of vertical slot fishways. The estuary con-
vironmental problems, and the presence of a non-stationary sidered has extensive flat marsh areas which flood and dry
wet-dry front limits the maximum computational time step periodically due to the tidal driven flow. This makes it pos-
in order to achieve accurate results. For these reasons depth sible to test the turbulence models in the presence of wet-
averaged models are still much more commonly used in en- dry fronts in complex two-dimensional geometries with a
gineering practise for the modelling of shallow flows. The very irregular bathimetry. The 90° bend has a simple geom-
depth averaged shallow water equations have been exten- etry, but at the same time it presents a recirculation bubble,
sively used in order to model the flow in rivers and coastal which makes it well suited for testing and comparing differ-
regions [65, 81, 83, 98], in open channels and hydraulic ent turbulence models. The flow in the fishway is highly tur-
structures [11, 36, 62, 100], as well as in flooding and drying bulent and it has strong recirculation regions, which makes it
problems [1, 13, 16, 70].
a perfect test case for the turbulence models. All the numer-
The modelling of turbulence in shallow water flows has
ical results are compared with extensive experimental data,
not been treated so profusely as in other fluid dynamics
which permits us to evaluate the degree to which the shal-
areas. Some depth averaged turbulence models have been
low water hypotheses are fulfilled in the considered flows.
proposed for the 2D-SWE. Those models derive from well
This is important because, as it is pointed out by Lloyd and
known RANS turbulence models, including somehow the
effects of bed friction in the turbulence field. Special men- Stansby [55], due to the assumptions made in the derivation
tion should be given to the depth averaged k–ε model pro- of the shallow water equations, and considering the numeri-
posed by Rastogi and Rodi in 1978 [76], which was the first cal dissipation inherent to the numerical schemes, the accu-
depth averaged two-equation eddy viscosity model, and it is racy of the results is problem dependent and usually uncer-
still the most commonly used with the 2D-SWE when tur- tain.
bulent effects are included in the computation. The rest of this paper is organised as follows. Section 2
The treatment of wet-dry fronts in the 2D-SWE has been summarises some previous theoretical, experimental and nu-
studied and improved in recent years by several researchers, merical studies about turbulent shallow water flows, as well
achieving non-diffusive and stable numerical schemes [13, as some features which distinguish them from 3D flows. The
16, 21, 34, 46, 70]. However, the modelling of turbulent derivation hypotheses made in the 2D-SWE are summarised
shallow flow with wet-dry fronts is seldom found in the lit- and discussed, in order to understand the limitations of the
erature. We have intended to give a forward step in that di- equations we are working with, and to be able to asses criti-
rection by using several turbulence models in problems with cally the numerical results obtained from them.
wet-dry fronts. Section 3 includes some notions about turbulence mod-
An unstructured finite volume model for quasi-2D turbu- elling in shallow waters. The depth averaged mixing length
lent free surface flow with wet-dry fronts is presented in this model as well as several versions of the depth averaged k–ε
paper. Several depth averaged turbulence models are consid- model are presented and discussed briefly.
ered, including a depth averaged mixing length model and In order to account for non-isotropic turbulence without
the depth averaged k–ε model of Rastogi and Rodi [76], with a significant decrease in the numerical stability, a depth av-
additional limiters to the production of turbulence. A new eraged algebraic stress model is proposed in Sect. 4 as an
depth averaged algebraic stress model is proposed and com-
extension of the 2D algebraic stress model, with additional
pared with the k–ε model.
terms which account for the production of turbulence due to
Regarding the numerical schemes, an hybrid second-
bed friction.
order/first-order scheme (first order in the water depth and
In Sect. 5, a finite volume solver for the turbulent shal-
second order in the unit discharge) is presented. The hybrid
low water equations is presented. The numerical schemes
scheme uses a second order discretisation for the two unit
discharge components, whilst keeping a first order discreti- and discretisation techniques implemented in the solver are
sation for the water depth. In such a way the numerical diffu- described in detail, and some recent advances in this field
sion is considerably reduced, without a significant reduction are discussed.
on the numerical stability of the scheme. In order to avoid Section 6 is intended to show the results obtained with
spurious oscillations of the free surface when the bathym- the model in three specific real engineering applications: the
etry is irregular, an upwind discretisation of the bed slope tidal flow in the Crouch estuary (Essex, England), the open
source term was implemented. This has proved to be more channel flow around a 90° bend, and the turbulent flow in
stable and accurate than a centred scheme [8, 38, 39]. two different designs of vertical slot fishways. The results
In order to show the capabilities of depth averaged mod- obtained with the mixing length, the k–ε and the algebraic
els in the computation of turbulent flows, the finite vol- stress models are discussed and compared with comprehen-
ume model developed has been used to model the flow in sive experimental data.
Depth Averaged Modelling of Turbulent Shallow Water Flow with Wet-Dry Fronts 305

Finally, the main conclusions and contributions of the pa- an important difference between quasi-2D shallow flows and
per are summarised in Sect. 7, and future research lines are 3D flows over 2D geometries, since the vortex stretching
proposed. process appears only in the latter.
The 2D properties are not always present in shallow wa-
ter flows. The stability of the large 2D structures depends on
2 Numerical Modelling of Shallow Flows in Hydraulic the balance between the effect of horizontal shear, which
Engineering produces 2D structures, the effect of shallowness, which
damps the 3D instabilities and promotes 2D turbulence, and
2.1 Shallow Water Flows the effect of bed friction and vertical shear, which damp
and stabilise the large 2D eddies by generating 3D turbu-
The study of turbulent shallow water flows is strongly linked lence [28, 92].
with 2D flow and 2D turbulence. Shallow water flows are Wolanski et al. [99] studied the wakes formed behind is-
those in which the vertical characteristic length scale is lands and decided to characterise the flow pattern by an is-
much smaller than the horizontal one. The water depth lim- land wake parameter P , which they defined as:
its the development of the 3D turbulent eddies in the ver-
UH2
tical direction, enhancing the 2D character of the turbulent P= , (1)
structures with a characteristic length larger than the water νt,z D
depth. The possibility of distinguishing between quasi-2D where U is the free stream velocity, H is the water depth,
and 3D turbulent structures depends on the ratio between D is the island characteristic diameter and νt,z is the vertical
the horizontal and the vertical turbulent length scales. The eddy viscosity, which can be approximated as:
2D state is approached as the water depth diminishes, but,
in any case, no matter how small the water depth is, there is νt,z = 1/6κuf H. (2)
always some interactions between the quasi-2D and the 3D
structures. From (2) and (1), we get:
There are several studies which reveal that the large hor-
6U H 2 6 U H
izontal quasi-2D coherent structures play an important role P= = . (3)
κuf H D κ uf D
in shallow water flows. Plane turbulent jets in shallow wa-
ters have been extensively studied experimentally by Dracos According to the study of Wolanski et al., when the island
et al. [35], Giger et al. [42], Thomas et al. [86] and Chen et wake parameter is very small (P  1) the flow is stable and
al. [26]. Shallow water wakes around circular cylinders and no bubble appears in the wake, for P ≈ 1 a stable turbu-
conical islands were studied by Lloyd et al. [55–57] as well lent wake is formed behind the island, and for P  1 the
as by Chen and Jirka [24, 25]. The homogeneous decaying bed friction becomes negligible and an unsteady wake ap-
turbulence produced by a grid in shallow flows was investi- pears. Further theoretical and experimental studies on the
gated by Uijttewaal and Jirka [91], who found the merging shallow water flow around islands have been made by Chen
of vortex which characterises the inverse energy transfer in and Jirka [24, 25]. In order to classify the shallow wake be-
2D turbulence. There are several studies by Uijttewaal et al. hind a cylinder they used a wake stability parameter S de-
[90, 92] and by Chu et al. [27] about the shallowness and bed fined by Ingram and Chu [49] as:
friction effects in the development of 2D turbulent structures
in free surface mixing layers. The turbulence characteris- D u2f D
S = 2cf =2 2 , (4)
tics in open shallow channels were measured and analysed H U H
by Prooijen et al. [94]. Different generation mechanisms of
where cf is the bed friction coefficient, defined as cf =
large coherent structures in shallow jets, wakes and mixing
τb u2
layers were proposed by Jirka [51]. ρU 2
= Uf2 , and τb is the bed shear stress. Chen and Jirka
All the previous theoretical and experimental studies found that above a critical vertical Reynolds number of
have revealed the critical role which large 2D turbulent Rv = UνH > 1500, the flow structures depend mainly on the
structures play in shallow water flows, as well as the impor- wake stability parameter S, being quite insensitive to the
tance of the bottom and free surface boundary conditions, horizontal Reynolds number (Rh = UνD ). This implies that
which confine the flow and modify the turbulence proper- even if shallow flows can be considered as 2D in the hor-
ties. The no-slip bottom condition enhances the 3D turbu- izontal plane, they should not be characterised by the hor-
lence production, while the slip condition at the free surface izontal Reynolds number (as it would be the case of 2D
promotes 2D turbulence [92]. The confinement of the flow flow), but the wake stability parameter should be used in-
by the bottom and free surface prevents the generation of stead. For small values of the stability parameter (S < 0.2) a
vorticity by the vortex stretching mechanism, which makes von Kármán’s vortex street appears behind the body. This is
306 L. Cea et al.

the case of a small bed friction and a large water depth. As This dimensional analysis confirms the results of Chen and
the bed friction increases and the water depth diminishes, Jirka [24, 25] regarding the importance of the wake stability
an unsteady bubble appears in the wake for values of the parameter in the analysis of shallow flows.
wake stability parameter in the range 0.2 < S < 0.5. For The results of Wolanski [99], Chen and Jirka [24, 25],
larger values of S the wake is stable and the bubble becomes and Lloyd [55, 57], show clearly that unlike in unbounded
steady. Similar results were obtained by Lloyd et al. [55, 57] flows, where the flow patterns are classified according to the
in flows around conical islands. Reynolds number, in shallow waters the flow patterns de-
The role played by the stability parameter S in shallow pend strongly on the water depth and on the bed friction.
flows can be explained by dimensional analysis of the shal- These parameters establish the differences between 2D un-
low water equations. Simple considerations about the char- bounded flow and 2D shallow water flow.
acteristic scale of each term in the x-momentum equation
gives: 2.2 RANS Models for Free Surface Flow
 
∂hU 2 ∂zs τb,i ∂ ∂U
= −gh − + (ν + νv + νh )h , The most general approach to model free surface flows, al-
∂x ∂x ρ ∂x ∂x
(5) though not the most often used, is to compute the fully 3D
U 2H gH 2 UH UH UH flow field with a specific treatment of the free surface bound-
∼ ∼ uf ∼ ν 2 ∼ uf H 2 ∼ U L 2 .
2
L L L L L ary. The main drawback to a fully 3D approach is its com-
putational cost, specially in environmental problems, where
In (5) only the terms in x-direction have been consid-
the size of the spatial domain is very large and there are
ered for simplicity in the notation. The total diffusion term
has been split in three terms which account respectively for flow patterns of different length scales involved in the flow.
viscous diffusion, turbulent diffusion due to 3D bed gener- For this reason, it is not yet efficient to use the fully 3D ap-
ated turbulence, and turbulent diffusion due to 2D horizontal proach in river engineering applications, as it is pointed out
turbulent structures. The diffusivity coefficient for each dif- in recent works by Duan [36] and Minh-Duc [61]. Still, the
fusion term is proportional to ν, uf H and U L respectively. 3D approach can be used to compute the flow in hydraulic
Considering that the convective term is of order θ (1), the structures. Olsen [66, 67] used a 3D-RANS model in order
similarity relations given by (5) can be expressed as: to compute the flow in a spillway, where the vertical ve-
locity is important, the extension of the spatial domain is
gH u2f L ν uf H small, and the geometry is simple. Wu et al. [101] used the
1∼ 2 ∼ 2 ∼ ∼ ∼ 1. (6) same approach to compute the flow and sediment transport
U U H LU U L
in open channels. Since 3D-RANS free surface models will
The following non-dimensional ratios have formed from
not be treated profusely in this paper, a brief description of
the previous dimensional analysis:
the current approaches is done in the following, in order to
gH 1 L S provide the interested reader with some basic concepts and
2
= 2, cf = , references for further reading.
U FR H 2
(7) The characteristic feature and main difficulty in 3D-
ν 1 √ H 1
= , cf = . RANS free surface models compared with models over a
U L Rh L Rv,t fixed geometry is the tracking of the unsteady free surface.
Considering usual values in hydraulic engineering, the There exist several numerical techniques to compute the free
order of magnitude of the previous ratios can be estimated surface position in 3D computations. A rather general classi-
as follows. The horizontal Reynolds number in environmen- fication distinguishes between mesh methods and meshless
tal flows is always much larger than Rh > 106 , and so the methods. Meshless methods use a Lagrangian formulation
viscous forces are negligible (except in the near wall re- in order to compute the movement of fluid particles applying
gion). The Froude number FR is usually smaller than 1 Newton’s Second Law. The most popular meshless method
and therefore, the hydrostatic pressure term is always im- nowadays is the Smoothed Particle Hydrodynamics (SPH)
portant and should always be considered. The other two method. It has the advantage of being able to treat compli-
non-dimensional ratios involve the bed shear stress. In en- cated free surface deformations, but it has problems with the
vironmental free surface flows the bed friction coefficient correct modelling of boundaries.
cf is generally in the range 0.005–0.01 [24, 51, 90], and Mesh methods can be classified in moving grid methods
the ratio HL is usually smaller than 0.1 in shallow flows.
and fixed grid methods. Moving grid methods use a La-
Then the wake stability parameter is in general larger than grangian formulation in order to move the grid nodes and
S > 0.1. The vertical turbulent Reynolds number is of the boundaries with the fluid [48, 89]. The free boundary is
order Rv,t ≈ 100 ( R1v,t ≈ 0.01) and thus, it has less impor- computed with a front tracking technique. A popular mov-
tance in the flow pattern than the wake stability parameter S. ing grid method for small deformations of the free surface is
Depth Averaged Modelling of Turbulent Shallow Water Flow with Wet-Dry Fronts 307

the Arbitrary Lagrangian-Eulerian method (ALE), in which the most popular ones, it is obvious that there are many pos-
the boundary of the grid moves with the fluid and the inner sibilities for computing the free surface in 3D models. De-
nodes of the mesh move arbitrarily in order to reduce mesh pending on the problem considered some methods are more
distortion. The advantage of this approach is that the free efficient than others. In the particular case of environmental
surface is resolved sharply, without diffusing it. The main free surface flows the most suited methods are probably the
disadvantage is the computational cost, since the nodes of height-function methods for cases in which the free surface
the mesh move at each iteration, and thus, the geometric can be defined by a single-valued function zs = f (x, y).
properties of the mesh need to be recomputed. Lagrangian This can be usually done in environmental problems, with
methods are mainly used when the movement of the free few exceptions, as for instance if we are interested in mod-
surface is small, because otherwise it is necessary to add or elling the shape of a breaking wave. In the case of a multi-
remove some nodes from the mesh in order to avoid a large valued free surface, the Volume of Fluid method in its orig-
distortion of the elements. In the latter case it is necessary to inal formulation (modelling only the fluid phase) is a sen-
remesh the whole domain, which is time consuming. sible option, and it has been reported to produce satisfac-
Fixed grid methods are more commonly used. They can tory results for modelling the flow in spillways [44, 77].
use a fully Eulerian formulation (interface capturing) or a At the present time, no published applications of the VOF
combined Eulerian-Lagrangian formulation (interface track- method for modelling the flow in rivers and coastal regions
ing). Fixed grid methods introduce a volume fraction func- are known to the authors.
tion which defines if an element is completely wet, partially In shallow water flows it is possible to simplify the 3D-
wet or completely dry. Among the many Eulerian methods RANS equations assuming an hydrostatic pressure distrib-
which have been proposed, it should be mentioned the Vol- ution. In such a case the vertical momentum equation is
ume of Fluid (VOF) method and the level set method, al- simplified to the hydrostatic pressure equation, and there-
though the latter one has some problems with mass conser- fore, only the two horizontal momentum equations need to
vation [41]. The VOF method is very popular, existing sev- be solved in a 3D mesh. The continuity equation is used
eral formulations [58, 78]. Some VOF formulations solve in order to compute the free surface level, which in turn,
only for the liquid phase, using a constant pressure bound- defines the hydrostatic pressure distribution. The numerical
ary condition at the free surface, while other VOF methods mesh is often built as a 2D horizontal mesh with several lay-
need to solve for both liquid and gas phases. Nonetheless, ers in the vertical direction, being the number of layers de-
pendent on the expected complexity of the vertical velocity
in its original formulation, proposed by Hirt and Nichols
profile. This way of building the 3D mesh is not compul-
in 1981 [44], the VOF method solves only for the liquid
sory, but it simplifies the mesh generation and also, it pro-
phase. The method solves a transport equation for a function
duces a well oriented mesh for stratified flows, which are
φ which takes the value 1 in the regions filled with fluid, and
rather common in environmental problems. This approach
the value 0 elsewhere. At each time step the transport equa-
is called 3D shallow water computations, and it has been
tion for φ is solved in a fixed mesh, the free surface is recon-
used by Bijvelds et al. [10] to compute the turbulent flow
structed according to the values of φ at the mesh nodes, and
in square harbours, by Lloyd and Stansby [55, 56] to model
suitable boundary conditions are applied at the free bound-
the shallow flow around conical islands, and by Stansby and
ary. In order to avoid loss of accuracy in the definition of
Lloyd [84, 85] to compute the wake around islands in oscil-
the interface due to numerical diffusion when advecting the
latory laminar and turbulent shallow flow. All the previous
function φ, a sharp reconstruction of the interface must be
simulations were done over simple or simplified geometries.
done at each time step. Several methods of reconstruction
Nonetheless, some works have been done using the 3D shal-
exist in the literature [3, 82], being that step an important low water equations to model the free surface flow in envi-
feature of the method. The fully Eulerian methods capture ronmental problems with complex geometries. In particular,
the interface using the function φ, that is why they are of- we should cite the works by Miglio et al. [60] and by Casulli
ten called interface capturing methods. On the other hand, and Stelling [18]. The latter model was extended to non-
combined Eulerian-Lagrangian methods use a front-tracking hydrostatic free surface flows by Casulli and Zanolli in [19].
technique (like Lagrangian methods) in order to follow the All of them use a height-function formulation in order to
free surface, although all the computations are done over a track the free surface.
fixed grid (like Eulerian methods). In this category fall the
height-function methods [63, 64] (which use a single-valued 2.3 Depth Averaged Models for Free Surface Flow
height function zs = f (x, y) in order to define the free sur-
face), the surface marker methods [63, 64], and the volume Further simplifications can be done in order to derive the
marker methods. depth averaged shallow water equations, also known as
From the previous description, which does not include all St. Venant equations or 2D shallow water equations (2D-
the existing methods for tracking the free interface, but only SWE), which are obtained after vertical integration of the
308 L. Cea et al.

3D shallow water equations (3D-SWE). The depth averaged the 2D-SWE. Uijttewaal and Tukker [92] studied the devel-
formulation has been successfully applied to different prob- opment of quasi-2D structures in a shallow free-surface mix-
lems, obtaining quite accurate results with a relatively low ing layer, concluding that the turbulence model used should
computational cost when compared with the 3D approach. account for both, the 3D turbulent structures created by bed
It has also the advantage of being very robust for comput- friction and the quasi-2D large structures originated by hor-
ing accurately the water depth, even in unsteady problems izontal strain. Hence, the effects of bed friction in the turbu-
with free surface shocks, as it happens to be the case in hy- lent kinetic energy must be included in the model. A classi-
draulic jumps and in dam break simulations. The treatment cal depth averaged turbulence model proposed by Rastogi
of unsteady wet-dry fronts, which appear typically in coastal and Rodi [76] is a depth averaged version of the famous
regions and in flooding problems, is also much simpler and k–ε model of Jones and Launder [52], with additional source
stable than in the 3D approach. The 2D depth averaged for- terms which account for the bed generated turbulence. Booij
mulation has been extensively used to model the dam break [12] proposed the modification of some constants in the bed
problem [7, 15], the propagation and runup of shallow water friction production terms of the Rastogi and Rodi model.
waves [21, 34, 46], flooding and drying problems [13, 16, Babarutsi and Chu [4] proposed a two-length-scale depth av-
17, 70], free surface flow in hydraulic structures [11, 23, 62], eraged version of the k–ε model which accounts in a differ-
ent way for the 3D bed generated turbulence. The same two-
flow in rivers and estuaries [22, 33, 83, 97, 98, 103], flow in
length-scale model, but with a slight modification on the dis-
coastal regions [81], and sediment transport in open chan-
sipation equation, was used by Babarutsi and Chu in [5]. All
nels and reservoirs [65, 100]. Another simplification step
these models will be presented and discussed in Sect. 3.
can still be given in order to obtain the 1D St. Venant equa-
tions, which can be applied to channels, hydraulic structures
2.4 Mathematical Approximations in Depth Averaged
or rivers, when the transversal effects are of little impor-
Models
tance.
Due to the physical assumptions done in the shallow wa- The mathematical derivation of the 2D-SWE can be found
ter models, the accuracy of the results is problem dependent. in many hydrodynamic books, and has already been pre-
Lloyd and Stansby [55] used both the 2D-SWE and the 3D- sented by many authors. Minor differences appear from one
SWE to compute the shallow water flow around conical is- derivation to another, but basically the process consists in as-
lands, and they found that in some cases the 2D model gives suming an hydrostatic pressure distribution, integrating the
more accurate results than the 3D model. The authors at- horizontal RANS equations over the water depth, applying
tributed those results to the fact that in a 2D model the verti- Leibnitz’s rule, and using the kinematic free surface and bed
cal mixing is instantaneous, while in a 3D model it depends surface conditions. Several approximations are done through
on the turbulence model used. the mathematical derivation. It is important to have in mind
Regarding turbulence, simple models are often used in these approximations in order to understand the limitations
hydraulic and coastal engineering. The extension of the spa- of the equations and to interpret the results obtained from
tial domain and the different scales involved preclude the them correctly. The approximations made are summarised
use of fully 3D-LES and DNS techniques in practical prob- in the following:
lems at the present time. The 2D features which are present
• Incompressible flow:
in most shallow water flows have induced Uittenbogaard
The incompressible RANS equations are taken as the base
and van Vossen to simulate the turbulent scales larger than equations for the derivation. By assuming incompress-
the water depth, rather than modelling them [93]. This ap- ible flow the density variations with the pressure gradients
proach leads to 2D-LES computations. This may seem con- are neglected, which is a reasonable hypothesis in water
tradictory with the idea of LES being always three dimen- flows.
sional, but it can be justified by using a filter size approx- • Hydrostatic pressure:
imately equal to the water depth, and by the fact that the The hydrostatic pressure distribution is mainly a conse-
horizontal large scales have 2D characteristics. Of course, quence of assuming a separation of the vertical and hor-
this approach is only appropriate when the water depth is izontal characteristic scales. This occurs when both the
extremely small, otherwise the filter size would not be ade- horizontal length scale Ln and the horizontal velocity
quate for LES purposes. scale Un are larger than the vertical ones Hn , Wn , which
The 2D-SWE are often used without any turbulence mod- is a typical characteristic of quasi-2D flows. The defini-
elling at all, in many cases directly neglecting the turbulent tion of the horizontal and vertical scales is not trivial, and
diffusive forces [13, 34, 46] or using a constant eddy viscos- depends on the flow conditions and geometry. In a long
ity [62]. Two-equation eddy viscosity models, specially k–ε shallow wave propagation problem the horizontal scale is
models, are usually the more sophisticated models used with given by the wave length, while the vertical scale is given
Depth Averaged Modelling of Turbulent Shallow Water Flow with Wet-Dry Fronts 309

by the water depth, since it is over those distances that in the vertical profile of the horizontal velocities u and v.
the velocity and pressure changes occur. In some cases Rastogi and Rodi [76] computed the shallow water flow in
the vertical length scale is given by the variations in the a longitudinal channel with secondary flows due to buoy-
bed and free surface elevation rather than by the water ancy effects, and found that the bed friction tends to damp
depth, and therefore, the separation of scales condition the secondary flows, diminishing in this way the relative
(Hn /Ln  1) is actually a restriction on the free surface importance of the dispersion terms. Hence, a better agree-
and bed slopes. Therefore, sometimes the 2D-SWE may ment between a 2D and a 3D model is obtained for rough
be applied to flows with a large water depth, assuming beds rather than for smooth beds.
that the bed and free surface slopes are small. Apart from The dispersion terms Duu and Duv are generally ne-
the separation of scales, there are another two conditions glected in the depth averaged equations, assuming the fol-
which must be fulfilled in order to assure an hydrostatic lowing approximation:
pressure distribution: the Reynolds number must be much
larger than 1, and the turbulence intensity must be smaller ∂Duu ∂Duv
+ = 0. (10)
than 1. Both conditions are usually fulfilled in shallow ∂x ∂y
water flows. Finally, it should be noticed that the fact of
In recent works, Duan [36] and Lien et al. [54] proposed
assuming a separation of length scales is equivalent to ne-
approximated expressions for the dispersion terms, based
glect the vertical accelerations, but it does not mean that
on experimental velocity profiles obtained in curved chan-
the vertical velocity is neglected.
nels. Both formulations are specially suitable for bend
• Homogeneous behaviour over the water depth:
channels for two reasons: first, due to the specific velocity
This approximation implies that the values of both the ve-
profile assumed over the water depth; and second, because
locity and the Reynolds stresses are almost independent
they need the definition of a streamwise direction, as well
of the vertical coordinate. This cannot be assumed as fre-
as an inner and outer bank which cannot be defined in a
quently as the hydrostatic pressure assumption. Vertical
general geometry. Both formulations are quite recent and,
integration over the water depth of the convective terms
to our knowledge, they have only been tested and used
in the horizontal RANS equations gives:
in bend channels, but not in general flow conditions, due
∂hU ∂hU 2 ∂hU V ∂Duu ∂Duv to the two reasons mentioned above. An alternative way
+ + + + (8) to account for the dispersion terms is to increase the ef-
∂t ∂x ∂y ∂x ∂y
fective eddy viscosity. However, in that case the value of
with the effective viscosity must be fixed on empirical basis,
 zs  zs and it depends on the problem geometry, as well as on the
2
Duu = u dz, Duv = u v  dz, flow conditions. Molls and Chaudry [62] use and effective
zb zb
  stress including the laminar stress, the turbulent stress and
zs zs
U= udz, V= vdz, the dispersion stress. The same idea is used by Minh-Duc
(9) et al. [61], who introduced a coefficient in the k–ε model
zb zb

u(x, y, z) = U (x, y) + u (x, y, z), in order to increase or decrease the eddy viscosity value.
• Depth averaged viscous and turbulent stresses:
v(x, y, z) = V (x, y) + v  (x, y, z), Slightly different formulations of the viscous and turbu-
lent diffusion terms can be found in the literature depend-
where h is the water depth, U and V are the depth aver-
ing on the author’s preference. The most common formu-
aged horizontal velocities, u and v  are the deviation of
lation is given by:
the horizontal velocities with respect to the depth aver-
aged velocities, zb is the bed elevation, and zs is the free  
∂ ∂Ui
surface elevation. (ν + νt )h , i = 1, 2, (11)
∂xj ∂xj
The fluxes Duu and Duv are often called longitudinal and
lateral dispersion stresses respectively. Their relative im- but it is also common to see the diffusion term expressed
portance respect to the convective and turbulent stress as:
terms depends on the magnitude of the velocities u and  
∂ ∂hUi
v  . In the limit case of a complete uniform velocity pro- (ν + νt ) , i = 1, 2, (12)
file over the water depth the dispersion terms vanish. In a ∂xj ∂xj
general case, their value is strongly dependent on the ex- or even as:
istence of secondary vertical currents, which appear typ-  
ically when the curvature effects in the velocity field are ∂ ∂Ui
h (ν + νt ) , i = 1, 2. (13)
important. These secondary flows create non-uniformities ∂xj ∂xj
310 L. Cea et al.

Even though there are some authors who justify one or the of a depth averaged turbulence model, which are usually de-
another formulation based on mathematical assumptions rived from RANS turbulence models by introducing in some
and approximations, in practise (11) is the most extended, way the effects of bed friction and shallowness in the turbu-
and to our criterion the most adequate, and therefore, it lence field. The most simple turbulence models used with
has been the one adopted in this paper. the 2D-SWE are the depth averaged parabolic model and
It has been assumed in the previous expressions that the the depth averaged mixing length model. Both of them con-
turbulent stresses are computed with an eddy viscosity sider an equilibrium state of turbulence, and compute the
turbulence model. In the case that an algebraic stress eddy viscosity with simple algebraic expressions assuming
model or a Reynolds stress model was used, the shear a turbulent length scale dependent on the water depth. Both
stress source term would be given by: are local models, in the sense that the eddy viscosity de-
  pends only on the local flow field, and it is not transported
∂ ∂Ui ∂    with the flow. The parabolic model accounts only for the tur-
νh − hui uj , i = 1, 2, (14)
∂xj ∂xj ∂xj bulence generated by bed friction, while the mixing length
model accounts also for the horizontal strain-rate produc-
where ui uj are the depth averaged Reynolds stresses. tion of turbulence. The parabolic model was used by Ding
and Wang [33], Duan [36], and Hsieh and Yang [45] in the
2.5 The Depth Averaged Equations computation of several open channel flows. In order to ac-
count for the production, transport, and dissipation of turbu-
With the former approximations, the turbulent depth aver- lence, Rastogi and Rodi [76] proposed a depth averaged k–ε
aged shallow water equations are obtained as: model for the 2D-SWE, which was used, among many oth-
ers, by Minh-Duc et al. [61] and Wu [100] to compute the
∂h ∂hUj sediment transport in open channels, and by Rameshwaran
+ = 0, (15)
∂t ∂xj and Shiono [75] to model the flow in meandering channels.

3.2 The Depth Averaged Parabolic Eddy Viscosity Model

∂hUi ∂hUi Uj
+
∂t ∂xj This is the simplest turbulence model for shallow flows. It
  assumes an equilibrium state of the flow, in which the verti-
∂h ∂zb τb,i ∂ ∂Ui
= −gh − gh − + νh cal shear stress and the hydrostatic pressure gradient are in
∂xi ∂xi ρ ∂xj ∂xj
perfect balance. By assuming uniform flow the convection
∂hui uj terms in the RANS equations can be neglected, and the x-
− , i = 1, 2. momentum equation reduces to a balance between vertical
∂xj
shear and pressure gradient:
Additional terms can be added to (15) in order to account
1 ∂P ∂τxz ∂u
for differences in atmospheric pressure and wind shear stress = = νt , (16)
in the free surface. Those terms have not been included in ρ ∂x ∂z ∂z
(15) for clarity in the notation. where ρ is the density of water, P = ρg(zs − z) is the hydro-
Leaving the turbulent terms apart, which must be com- static pressure, zs is the free surface elevation, and νt (z) is
puted with a suitable depth averaged turbulence model, the the eddy viscosity (not its depth averaged value). Assuming
2D-SWE are a system of 3 partial differential equations a logarithmic velocity profile in depth, the model predicts a
with 3 unknowns (U, V , h), which are defined over a two- parabolic profile in the vertical coordinate for the eddy vis-
dimensional spatial domain. This is an important reduction cosity [20], which is given by:
in the computational cost with respect to the original 3D-  
z
RANS equations, which are 4 equations defined over a 3D νt (z) = uf κz 1 − , (17)
spatial domain, with the additional inconvenience of the free h
surface moving boundary. where uf is the bed friction velocity κ = 0.41 is the von
Karman’s constant, and h is the water depth. The depth av-
eraged eddy viscosity is obtained after vertical integration
3 Turbulence Modelling in Shallow Water Flows of (17) over the water depth, which yields:
  
3.1 Depth Averaged Turbulence Models for Shallow Flows 1 zs z
νt = uf κz 1 − dz
h zb h
The depth averaged horizontal Reynolds stresses u2 , u v  , 1
= κuf h = 0.068uf h. (18)
v 2 appearing in the 2D-SWE need to be computed by means 6
Depth Averaged Modelling of Turbulent Shallow Water Flow with Wet-Dry Fronts 311

The parabolic eddy viscosity model does not account for ls dependent on the water depth can often lead to an under-
the effect of horizontal velocity gradients, but only for the estimation of νt , specially in flows in which turbulence is
turbulence generated by bed friction. It leads to very low mostly 2D.
eddy viscosity values when the turbulence production due to Equation (23) is not valid near the walls, as it would pre-
horizontal shear is important. It does not account for trans- dict very large length scales. Instead, it is more correct to use
port and dissipation processes. Despite its simplicity, it has the wall distance as the turbulent length scale in the near wall
been used by Hsieh and Yang [45], Lien et al. [54] and regions. The final expression for the eddy viscosity given by
Duan [36] in the computation of simple channel flows. It the model is:
is also used by more complex models in order to account for
       
the turbulence produced by bed friction. ∂U 2 ∂V 2 ∂U ∂V 2 uf 2
νt = ls2 2 +2 + + + 2.34 ,
∂x ∂y ∂y ∂x κh
3.3 The Depth Averaged Mixing Length Model ls = min(0.267κh, κdwall ), (24)

This is a depth averaged version of the original mixing where dwall is the distance to the nearest wall. In the ab-
length model. In order to account for both the horizontal and sence of horizontal velocity gradients, (24) and (18) are
the vertical production of turbulence, the total eddy viscos- equal. Hence, the depth averaged mixing length model tends
ity is split into an horizontal νth and a vertical νtv component. to the parabolic eddy viscosity model when turbulence is
The horizontal eddy viscosity accounts for the turbulence mainly produced by bed friction. This is the case of flows
produced by horizontal shear. It is computed as: with a very low water depth, a relatively rough bed surface,
  and small horizontal velocity gradients. The depth averaged
νth = ls2 2Sij Sij = ls2 2(S11 2 + S 2 + 2S 2 ),
22 12 (19) mixing length model was used in channel flow computations
by Jia and Wang [50], among others.
where ls is the characteristic horizontal turbulent length
scale, and Sij is the horizontal mean strain-rate tensor com-
3.4 k–ε Models
puted from the depth averaged velocity as:
 
1 ∂Ui ∂Uj 3.4.1 The k–ε Model of Rastogi and Rodi
Sij = + . (20)
2 ∂xj ∂xi
The first k–ε model for shallow water flows was proposed by
On the other hand, the vertical eddy viscosity is generated
Rastogi and Rodi [76] as a depth averaged version for quasi-
by the vertical velocity gradient produced by bed friction.
2D flows of the original k–ε model of Jones and Launder
Therefore, it is computed from the parabolic eddy viscosity
[52]. Instead of solving for the three-dimensional turbulence
model as:
and dissipation, it solves for their depth averaged value.
1 The equations of the model are very similar to those of
νtv = κuf h. (21)
6 the standard 2D k–ε model. The main difference is that it
is necessary to introduce a production term in order to ac-
The total eddy viscosity is evaluated from the horizontal
count for the production of turbulence due to bed friction.
and vertical values as:
The equations of the model are given by:
 
   
1 κuf h 2 ∂hk ∂Uj hk ∂ νt ∂k
νt = (νth )2 + (νtv )2 = ls2 2Sij Sij + . (22) + = ν+ h
6 ls2 ∂t ∂xj ∂xj σk ∂xj
In the inner part of the domain the turbulent length scale + hPk + hPkv − hε,
ls is assumed to be dependent on the water depth, since it   
∂hε ∂Uj hε ∂ νt ∂ε
is a restriction on the size of the turbulent eddies. In order + = ν+ h
to estimate ls , a fictitious turbulent length derived from the ∂t ∂xj ∂xj σε ∂xj
parabolic eddy viscosity model is averaged over the vertical ε ε2
coordinate, which yields [20]: + hcε1 Pk + hPεv − hcε2 ,
k k
ls ≈ 0.267κh. (23) k2
νt = cμ , Pk = 2νt (Suu
2
+ Svv
2
+ 2Suv
2
), (25)
ε
It should be remarked that, as it was pointed out in u3f 1
Sect. 2, there might be 2D turbulent structures with a ls Pkv = ck , ck = 1/2
,
h cf
larger than the water depth. Therefore, the assumption of
312 L. Cea et al.

u4f 1/2
cε2 cμ value of e∗ = 0.11 in (29). In the results presented in this pa-
Pεv = cε , cε = 3.6 , per the original definition of cε was used, since (29) would
h2 3/4
cf
need a previous calibration for each practical application.
cμ = 0.09, cε1 = 1.44, cε2 = 1.92,
3.4.2 The Modified Constants of Booij
σk = 1.0, σε = 1.31,

where cf is the bed friction coefficient (u2f = cf |U|2 ). The The production terms Pkv and Pεv in (25) depend on two
coefficients, ck and cε . The value of these coefficients deter-
five constants of the model (cμ , σk , σε , c1ε , c2ε ) are assumed
mine the turbulence level in uniform channel flow conditions
to have the same values as in the original k–ε model. The
(see (27)). Booij [12] proposed the following modification in
term Pk accounts for the production of turbulent energy due
the value of ck and cε :
to horizontal velocity gradients, and it has the same expres-
sion as in the 2D model. The source terms Pkv and Pεv are 1 1
responsible for modelling the 3D turbulence generated by ckB = ck , cεB = cε , (30)
10 44
bed friction.
where the superindex B refers to Booij. The effect of dimin-
In uniform channel flow conditions turbulence is in an
ishing the coefficients ck and cε in such a way is that the
equilibrium state. In that case the production of turbulent
production terms of both k and ε due to bed friction (Pkv
kinetic energy is due to bed friction, and the k–ε equations
and Pεv ) are strongly reduced. Using the modified coeffi-
reduce to:
cients of Booij to evaluate the turbulent kinetic energy and
ε2 dissipation in uniform channel flow (see (27)) gives:
Pkv = ε, Pεv = cε2 . (26)
k
kuB = 0.44ku , εuB = 0.10εu , B
νt,u = 1.94νt,u , (31)
Using the expressions for Pkv and Pεv given by (25), the
values for the turbulent energy and dissipation predicted by where the values ku , εu , νt,u are those given by the original
the model in uniform channel flow are given by: constants used by Rastogi and Rodi [76]. Hence, in uniform
channel flow conditions, the modification proposed by Booij
3/2
cε2 2 2 |U|1/2 uf reduces the turbulent kinetic energy by a factor of 0.44, but
ku = ck uf = , it almost duplicates the eddy viscosity. On the other hand,
cε 1.08
(27) in flows dominated by horizontal shear the differences with
u3f |U|u2f the Rastogi and Rodi model smear, because the term Pkv
εu = ck =
h h becomes much smaller than Pk in (25).
and the eddy viscosity is equal to: 3.4.3 The k–ε Model of Babarutsi and Chu
2
cε2
νt,u = cμ ck3 uf h = 0.08uf h (28) Babarutsi and Chu [4] proposed a two-length-scale depth
cε2 averaged k–ε model. The main difference with Rastogi’s
model is the way in which the effects of the 3D turbulence
which is a very similar value to that one given by the depth
produced by bed friction are introduced. Babarutsi and Chu
averaged parabolic and mixing length models (see (18)). On
propose to split the total eddy viscosity into a 3D part and a
the other hand, when the horizontal shear is much larger than
2D part as:
the vertical shear, the vertical production terms are negligi-
ble compared to the horizontal shear production (Pkv ≈ 0, νt = νt2D + νt3D . (32)
Pεv ≈ 0), and the model reduces to the standard 2D k–ε
model. The 3D eddy viscosity accounts for the small scale
A slight modification of the coefficient cε has been used bed generated turbulence. The equation for νt3D used by
in recent works by Minh-Duc et al. [61]. The expression Babarutsi and Chu is given by:
used is given by:
νt3D = 0.08uf h (33)
1/2
1 cε2 cμ
cε = , (29) which is the same value given by the model of Rastogi and
(e∗ σt )1/2 c3/4 Rodi for uniform channel flow (see (28)). The 2D part is
f
computed from k  and ε  , which account for the large-scale
where σt is a Prandtl/Schmidt number relating the eddy vis- 2D turbulence generated by horizontal shear, as:
cosity and diffusivity for the transport of scalars (σt = 0.7 in
[61]), and e∗ is an adjustable empirical parameter. The con- k 2
νt2D = cμ . (34)
stant 3.6 used in the definition of cε in (25), corresponds to a ε
Depth Averaged Modelling of Turbulent Shallow Water Flow with Wet-Dry Fronts 313

It should be remarked that the turbulent kinetic energy k  the turbulence is 2D, and generated by horizontal shear. The
does not include the 3D bed generated turbulence. Hence, equations for k  and ε  reduce to the standard 2D k–ε model.
the transport equations used to compute k  and ε  are similar This model has been used by Babarutsi et al. in [6], with
to those used in Rastogi’s model (see (25)), but the vertical the constant cε3 = 1, to compute shallow recirculating flows
production terms are zero (Pkv = Pεv = 0) and a new term dominated by friction. It has also been used by Babarutsi
F  is introduced in order to account for the transfer of en- and Chu to model shallow mixing layers [5], where they ob-
ergy between the large 2D turbulent scales and the small 3D tained the same results with cε3 = 1 and with cε3 = 0.8.
scales. The equations for k  and ε  are given by:
   3.4.4 Differences between the 3 Versions of the k–ε Model
∂hk  ∂Uj hk  ∂ νt ∂k 
+ = ν+ h
∂t ∂xj ∂xj σk ∂xj The main differences between the previous depth averaged
  k–ε models appear in friction dominated flows. The modifi-
+ h(Pk  − F ) − hε ,
cation of Booij reduces the bed generated turbulence, and as
∂hε  ∂Uj hε  a consequence, the large scale turbulence as well as the total
+ (35)
∂t ∂xj turbulent diffusion of the model are increased. We have not
  
∂ νt ∂ε  found concluding results showing that the modification pro-
= ν+ h posed by Booij improves the results of the original model of
∂xj σε ∂xj
Rastogi in general flow conditions. Furthermore, the work
ε  ε 2 by Minh-Duc et al. [61] seems to show that the optimum
+ hcε1 (Pk  − (1 − cε3 )F ) − hcε2 .
k k value of the constant cε might be problem dependent.
The version of Babarutsi and Chu is more interesting
A value of cε3 = 0.8 is proposed in [4]. Note that in (35)
from a conceptual point of view, since it distinguishes be-
the diffusion term is modelled with the total eddy viscosity
tween two different turbulent length scales. The main differ-
νt . On the other hand, the production term Pk  is computed
ences with Rastogi’s model appear in shear flows dominated
with (25), but using only the 2D eddy viscosity νt2D . All the
by friction. Both models converge to the standard 2D k–ε
constants of the model are the same as in Rastogi’s model.
model for zero bed friction. On the other hand, in uniform
The sink term F  is associated with the negative work done
channel flow, where all the turbulence is generated by bed
by the large scale turbulent fluctuations against the bed fric-
friction, both models give the same eddy viscosity value (see
tion, i.e. the bed friction takes out energy from the large 2D
(28) and (37)). It is in horizontal shear layers dominated by
turbulent structures. It is computed as:
bed friction where both models show some differences. Al-
cf [u2 (2U 2 + V 2 ) + 2u v  U V + v 2 (U 2 + 2V 2 )] tai et al. [2] estimate that, in order to produce the modelled
F = √ , effects introduced by the term F  , the water depth should be
h U2 + V 2 (36) approximately 100 times smaller than the horizontal turbu-
lent length scale, i.e. the effects of F  are only significant
where the Reynolds stresses u2 , u v  , v 2 are computed from in flows which are strongly dominated by friction and at the
the Boussinesq approximation using the large scale eddy same time have large 2D turbulent structures. Babarutsi et
viscosity νt2D . al. [6] compared the results given by the models of Babarutsi
In uniform channel flow, the production of large scale and Rastogi in a recirculating shallow water flow dominated
turbulence due to horizontal shear is zero (Pk  = 0). The by friction. Despite the different eddy viscosity predicted by
large scale turbulent kinetic energy equation reduces to 0 = the models, no significant differences in the mean flow field,
−F  − ε  . Therefore, all the large scale turbulent variables neither in the recirculation length, were found, and there-
are zero, i.e.: fore, no conclusion was obtained regarding which model
performs better. Both models were compared again by
ku = 0, εu = 0, 2D
νt,u = 0, Babarutsi and Chu [5] in order to model transverse mixing
(37)
νt,u = νt3D = 0.08uf h, layers in shallow flows dominated by friction. In that case
they found that the results given by the model of Babarutsi
where the subindex u refers to uniform channel flow con- were in better agreement with the experimental data.
ditions. In this case all the turbulence is 3D, and generated All the results presented in this work have been obtained
by bed friction. The total eddy viscosity is similar to that with the original depth averaged model of Rastogi and Rodi.
one given by the model of Rastogi. On the other hand, when The reason of using the original model of Rastogi and Rodi
the bed friction velocity uf tends to zero, the flow is domi- is that it has been used in many river flow and channel sim-
nated by horizontal shear. In such a case the term F  and the ulations, showing a good behaviour [10, 61, 68, 75, 97, 100,
3D eddy viscosity may be neglected (F  ≈ 0, νt3D ≈ 0), all 102]. On the other hand, the versions of Booij and Babarutsi
314 L. Cea et al.

have not been so extensively used. In any case, in the appli- 3.4.6 Realizability in the Eddy Viscosity Models
cations presented in this work no differences are expected to
be found in the mean flow field obtained with the different Another important drawback to Boussinesq assumption,
versions of the model. This is because the flow in the verti- which is somehow related to the excessive turbulence pre-
cal slot fishway, as well as in the open channel 90° bend, is diction in non-isotropic turbulence, is that it can produce
not dominated by bed friction, while in the Crouch estuary, negative values of the normal Reynolds stresses, specially
although the turbulence is mainly generated by bed friction, in stagnation regions. This was pointed out by Durbin in
there are not strong horizontal shear stresses. [37]. In order to avoid negative normal stresses, a constraint
should be used in the eddy viscosity value. This procedure
3.4.5 Limiter to the Production of Turbulent Kinetic Energy also helps to avoid the excessive turbulent production near
stagnation regions.
It is well known that Boussinesq assumption does not work In 2D flow, the constraint over the eddy viscosity which
well when the turbulence is strongly non-isotropic, as it hap- assures that all the normal turbulent stresses remain positive
pens in swirling flows, in regions with strong adverse pres- is given by:
sure gradients, or near stagnation points. The main cause of  1/2
these problems is that the eddy viscosity is assumed to be k 2
νt < . (40)
isotropic. In depth averaged models, the non-isotropy be- 3 Sij Sij
tween the vertical and horizontal directions is somehow in-
The limiter given by (40) should always be used in all
troduced in the equations by the different ways of account-
turbulence models based on Boussinesq assumption.
ing for the production of 2D horizontal turbulence and 3D
turbulence generated by bed friction. However, it still re-
mains the problem with non-isotropy in the horizontal plane.
4 A Depth Averaged Algebraic Stress Turbulence
It is well known that in stagnation points the produc-
Model
tion of turbulent kinetic energy is overpredicted by the k–ε
model, which gives much larger values of turbulent kinetic
One of the main weakness of eddy viscosity models is
energy than the measured ones [59]. For this reason, a lim-
Boussinesq approximation, which assumes an isotropic
iter in the production term is introduced in almost all the
eddy viscosity when computing the Reynolds stresses. This
eddy viscosity models in order to avoid too large values of
assumption causes a reduction in the accuracy of the model
the turbulent energy in stagnation regions. Menter [59] lim-
when the turbulence field is strongly non-isotropic, as for
its the ratio between the production and dissipation of turbu-
example in swirling and recirculating flows, in stagnation
lence as:
regions, or in flows with strong adverse pressure gradients.
P = min(Pk , cl ε), (38) In order to improve the turbulence modelling in those flow
conditions, the Reynolds Stress Models (RSM) solve one
where cl is a constant with a value between 10 and 20. This transport equation for each Reynolds stress, which avoids
limiter is a rough estimation based on empirical observations the need of using Boussinesq assumption in the computa-
which does not imply a good modelling of the stagnation tion of the turbulent stresses. Reynolds Stress Models ac-
region, but it can improve the numerical results, avoiding an count much better for non-isotropic turbulence, but they are
extremely large turbulence level which may propagate and more expensive and unstable from a computational point
affect the solution in the whole numerical domain. of view. An alternative approach are the Algebraic Stress
In this work the turbulent production limiter has been ap- Models (ASM) introduced by Rodi in [79], in which the
plied independently to both, the horizontal shear production Reynolds stresses are computed by means of algebraic non-
and the bed friction production as: linear equations. This approach reduces the computational
cost considerably. Algebraic Stress Models can be regarded
Pk = min(Pk , cl ε), either as a simplification of RSM or as a non-linear exten-
(39)
Pkv = min(Pkv , cl ε) sion of eddy viscosity models.
Algebraic Stress Models are quite unstable, specially in
with a value of cl = 10. The limit imposed in the bed friction 3D flow computations, due to the highly non-linear implicit
production Pkv is specially suited for problems with wet-dry expressions used to evaluate the Reynolds stresses. The sta-
fronts. The term Pkv can take quite large values when the bility of the model is increased if explicit rather than implicit
water depth is very small. In that case, the second limiter in expressions are used to compute the turbulent stresses [30].
(39) avoids instabilities in the solution, and allows the use In this section an explicit depth averaged Algebraic Stress
of smaller values of the wet-dry tolerance parameter (see Model (DASM) which accounts for the turbulence produc-
Sect. 5). tion due to bed friction is proposed. The model proposed
Depth Averaged Modelling of Turbulent Shallow Water Flow with Wet-Dry Fronts 315

shows a numerical stability and computational cost similar energy is computed directly from the Reynolds stresses as:
to that of the k–ε model.
Puu + Pvv + Pww
Pk =
2
4.1 Algebraic Stresses Model for 3D Flow  
∂U ∂U ∂V
= − u2 − u v  +
∂x ∂y ∂x
The expressions of the Algebraic Stress Model proposed in ∂V Puu,V + Pvv,V
[79] for 3D flow are given by: − v 2 + . (43)
∂y 2
Introducing (42) and (43) into the algebraic expressions for
k (1 − c2 )(Pij − 3 Pk δij ) + ij,R + ij,S
2 w w
2
ui uj = kδij + . the Reynolds stresses (see (41)) and neglecting the wall-
3 ε c1 + Pεk − 1 correction pressure-strain terms yields:

(41)
 
k (1 − c2 ) ∂V ∂U ∂U
uv = −u2 − v 2 − u v 
ε c11 ∂x ∂y ∂x
The production terms Pk and Pij , as well as the wall- 
∂V
correction pressure strain terms, w w
ij,R and ij,R , depend on − u v  + Puv,V ,
∂y
ui uj . Hence, there is a strong non-linear relation between all 
the Reynolds stresses, which promotes the instability of the 2 k (1 − c2 ) 4 ∂U 4 ∂U
u2 = k + − u2 − u v 
model, being the numerical convergence usually more unsta- 3 ε c11 3 ∂x 3 ∂y
 (44)
ble than in RSM. Explicit expressions have been derived for 2 2 ∂V 2   ∂V 2 1
2D [71] and 3D [40] flow. The expressions for the Reynolds + v + uv + Puu,V − Pvv,V
3 ∂y 3 ∂x 3 3
stresses in the explicit ASM are more complex, specially in 
2 k (1 − c2 ) 4 ∂V 4 ∂V
3D flow, but the numerical stability is improved [30]. v 2 = k + − v 2 − u v 
3 ε c11 3 ∂y 3 ∂x

2 2 ∂U 2   ∂U 2 1
4.2 Algebraic Stresses in 2D Shallow Flows + u + uv + Pvv,V − Puu,V
3 ∂x 3 ∂y 3 3

Only three Reynolds stresses (u2 , u v  and v 2 ) appear in with the constants:
the 2D-SWE. In order to account for the bed friction effect, Pk
the production term for each Reynolds stress will be split c1 = 1.8, c2 = 0.6, c11 = c1 + − 1. (45)
ε
into an horizontal 2D production (Puu,H , Pvv,H , Puv,H ) and
The system of equations 44 can be written in matrix form
a vertical production due to bed shear (Puu,V , Pvv,V , Puv,V ).
as:
Assuming the separation between horizontal and vertical ⎛ ⎞ ⎛ ⎞ ⎛ 2 ⎞ ⎛ ⎞
scales in shallow flows, i.e. that the vertical velocity is much u2 a11 a12 a13 u b1
⎜ 2 ⎟ 1 − c k
2 ⎝ ⎜ ⎟
smaller than the horizontal one, the expression for each pro- ⎝ v ⎠= a21 a22 a23 ⎠ ⎝ v 2 ⎠ + ⎝ b2 ⎠
 
c11 ε
duction term is given by: uv a31 a32 a33 u v  b3
(46)
∂U ∂U
Puu = Puu,H + Puu,V = −2u2 − 2u v  + Puu,V , with the coefficients aij and bi equal to:
∂x ∂y
∂V ∂V 4 ∂U 2 ∂V
Pvv = Pvv,H + Pvv,V = −2u v  − 2v 2 + Pvv,V , a11 = − , a12 = ,
∂x ∂y 3 ∂x 3 ∂y
(42) 4 ∂U 2 ∂V
Puv = Puv,H + Puv,V a13 = − + ,
3 ∂y 3 ∂x
 
∂V ∂U ∂V ∂U
= −u2 
−u v  + − v 2 + Puv,V . 2 ∂U 4 ∂V
∂x ∂x ∂y ∂y a21 = , a22 = − ,
3 ∂x 3 ∂y
4 ∂V 2 ∂U (47)
Since the vertical velocity is assumed to be negligible, a23 = − + ,
3 ∂x 3 ∂y
the production of w 2 is zero (Pww = 0). The evaluation of
the terms Puu,V , Pvv,V , Puv,V , which account for the pro- ∂V ∂U
a31 = − , a32 = − ,
duction due to vertical shear, will be treated later on in this ∂x  ∂y
section. Instead of using the Boussinesq eddy viscosity ap- ∂U ∂V
a33 = − + ,
proximation, in the ASM the production of turbulent kinetic ∂x ∂y
316 L. Cea et al.
 
2 k (1 − c2 ) 2 1 4.3 Estimation of the Production of Turbulence Due to Bed
b1 = k + Puu,V − Pvv,V ,
3 ε c11 3 3 Friction
 
2 k (1 − c2 ) 2 1
b2 = k + Pvv,V − Puu,V , (48) The only thing which remains to specify in the present
3 ε c11 3 3
model is the evaluation of the vertical production terms due
k (1 − c2 ) to bed friction (Puu,V , Pvv,V , Puv,V ). The vertical produc-
b3 = Puv,V .
ε c11 tion of Reynolds stresses in (42) is given by:
After simple algebraic manipulation, the system of equa- ∂u ∂v
tions (46) can be expressed as: Puu,V = −2u w  , Pvv,V = −2v  w  ,
∂z ∂z
⎛ ⎞
  u2 Pww,V = 0,
k ⎜ ⎟ (53)
c11 I − (1 − c2 ) A ⎝ v 2 ⎠ = c11 b, (49) ∂v ∂u ∂u
ε Puv,V = −u w  − v  w , Puw,V = −w 2 ,
u v  ∂z ∂z ∂z
∂v
where A is a matrix given by the expressions (47), and b is Pvw,V = −w 2 .
∂z
the vector given by expressions (48). The coefficient c11 is
given by (45). At this point two possibilities arise when com- Since the production of w 2 is zero, all the generation
puting the production of turbulent kinetic energy Pk in (45). of turbulent kinetic energy is distributed on the horizontal
The first option is to compute Pk using (43). If this is done
Reynolds stresses u2 and v 2 . Hence, the vertical production
the system of equations (49) is non-linear, and an iterative
of k is given by:
procedure is needed to solve it. This leads to a more unsta-
ble and expensive numerical scheme, because the coupling Puu,V + Pvv,V
Pkv = . (54)
between Reynolds stresses in very strong. Alternatively, Pk 2
can be computed using the eddy viscosity assumption as:
It should be remarked that even if the production of w 2
      
∂U 2 1 ∂U ∂V 2 ∂V 2 is zero, this does not mean that its value is zero, since the
Pkνt = 2νt + + +
∂x 2 ∂y ∂x ∂y algebraic expressions (41) account for pressure-strain distri-
bution between the normal stresses.
+ Pkv (50)
In order to approximate the three vertical production
with the eddy viscosity νt computed from the k–ε model. terms (Puu,V , Pvv,V , Puv,V ), uniform channel flow condi-
With this choice, the system of equations (49) can be solved tions will be assumed. Under these flow conditions the verti-
explicitly, improving in this way the convergence of the nu- cal velocity, as well as the horizontal velocity gradients, are
merical scheme. The approximation given by (50) is not ex- zero (w = 0, ∂x∂
= ∂y
∂
= 0), and all the production of turbu-
pected to influence importantly the results, since it is just lence is due to vertical shear created by bed friction. A hor-
used in order to compute the constant c11 and therefore, its izontal rotation of the reference coordinate system will be
influence is smoothed. Therefore, (50) will be used for com- considered, in such a way that in the new coordinate system
puting the turbulent kinetic energy production term in the the transverse velocity is zero (v  = 0). The relation between
coefficient c11 (see (45)). The system of equations (49) can the velocities and Reynolds stresses in both coordinate sys-
be written in tensorial form as: tems is given by a horizontal rotation of α degrees as:
⎛ ⎞ ⎛ ⎞
mij rj = c11 bi , (51) u u
⎝ v ⎠ = NT ⎝ v  ⎠ , (55)
where rj (j = 1, 3) are the horizontal Reynolds stresses w w
(u2 , v 2 , u v  ), and the tensor mij is equal to:
where NT is the transposed of the rotation matrix N:
k ⎛ ⎞
mij = c11 δij − (1 − c2 ) aij . (52) cos α −sinα 0
ε
N = ⎝ sin α cos α 0 ⎠ . (56)
The Reynolds stresses are obtained after solving the sys-
0 0 1
tem of equations (51). It should be remarked that even if the
Reynolds stresses u w  , v  w  and w 2 do not appear in the In a similar way, the relation between the turbulent stresses
2D-SWE, that does not mean that their value is assumed to is obtained by a horizontal rotation of the Reynolds stress
be zero. These turbulent stresses can be evaluated from the tensor as:
general equation (41) in a similar way as it has been done
for the horizontal stresses. T = NT T N, 
tij = nki tkl nlj , (57)
Depth Averaged Modelling of Turbulent Shallow Water Flow with Wet-Dry Fronts 317

where T and T are the Reynolds stress tensors in both co- Finally the vertical production terms are evaluated as:
ordinate systems:
u2f u2x u2f u2y
⎛ ⎞ ⎛ 2 ⎞ Puu,V = 2 , Pvv,V = 2 ,
t11 t12 t13 u u v  u w  h|U| h|U|
⎜ ⎟
T = ⎝ t21 t22 t23 ⎠ = ⎝ u v   v 2 v  w ⎠ . (58) (65)
u2f ux uy
t31 t32 t33 u w  v  w  w 2 Puv,V = −2 .
h|U|
The main property of the new reference system is that the
4.4 Realizability Condition in the DASM
transverse velocity vanishes (v  = 0), and as a consequence
 also vanishes:
the turbulent stress t23 It is necessary to use a realizability condition with the
  DASM, since there is nothing in the equations of the model
 k 1 − c2 ∂v 
t23 = −w 2 = 0. (59) which assures that the normal Reynolds stresses will re-
ε c11 ∂z main positive. Hence, the following restrictions should be
imposed in the model:
Considering this property, the vertical production of u2
can be computed as: u2 ≥ 0, v 2 ≥ 0, u2 + v 2 ≤ 2k. (66)
⎛ ⎞
0 The third condition in (66) assures that the vertical nor-
∂u
Puu,V = −2t13 = −2(cos α, sin α, 0)T ⎝ 0 ⎠ mal stress w 2 remains positive, and avoids excessively large
∂z
1 values of the horizontal normal stresses. If any of the normal
   stresses is negative its value is set to zero. At the same time,
∂u ∂v 
× cos α + sin α in order to keep constant the total turbulent kinetic energy,
∂z ∂z
⎛  ⎞ (60) it is necessary to subtract from the other normal stresses the
t13  same amount of energy. This is done, for the stress u2 , in
 ⎠ ∂u cos α
= −2 (cos α, sin α, 0) ⎝ t23 the following way:
 ∂z
t33
⎧
 ⎪
⎪ u2 = 0,
 ∂u ⎪
= −2t13 2
cos α. ⎪
⎪
∂z ⎪
⎨ 2 v 2
v = v 2 − a 2 ,
If u2 = −a 2 < 0 then v 2 + w 2
In a similar way, the vertical productions of v 2 and u v  ⎪
⎪
⎪
⎪
are obtained as: ⎪
⎪ w 2
⎩w 2 = w 2 − a 2 .
∂v  v 2 + w 2
 ∂u
Pvv,V = −2t23 = −2t13 sin2 α, (67)
∂z ∂z
(61)

∂v ∂u  ∂u It is straightforward to show that (67) conserves the turbu-
Puv,V = −t13 − t23 = 2t13 sin α cos α.
∂z ∂z ∂z lent kinetic energy. A simple way to implement the realiz-
ability condition in the numerical solver is to define the fol-
From (54), (60) and (61), the vertical production of turbulent
lowing ratios:
kinetic energy is equal to:

Puu,V + Pvv,V  ∂u Ru =
2k
Rv =
2k
Rw =
2k
Pkv = = −t13 . (62) , ,
2 ∂z 2k − u2 2k − v 2 2k − w 2
Using (62) in (60) and (61) yields: (68)

Puu,V = 2Pkv cos2 α, and to recompute the Reynolds stresses as:

Pvv,V = 2Pkv sin2 α, (63) u2 = max{0, min(u2 , Rv u2 , Rw u2 )},
(69)
Puv,V = −2Pkv sin α cos α. v 2 = max{0, min(v 2 , Ru v 2 , Rw v 2 )}.

The production of turbulent kinetic energy due to bed

friction will be approximated for uniform channel flow in 5 A Finite Volume Solver for the Turbulent Shallow
the same way as it is done in the k–ε model: Water Equations

|U|u2f This section presents an unstructured finite volume model

Pkv = . (64)
h for solving the 2D-SWE and turbulence models presented
318 L. Cea et al.
⎛ ⎞ ⎛ ⎞
0 0
⎜ −gh ∂zb ⎟ ⎜ − τb,x ⎟
S=⎝ ∂x ⎠ , G = ⎝ ρ ⎠,
τb,y
−gh ∂z ∂y
b
− ρ
⎛ ⎞
0
⎜ ∂ (νe h ∂Ux ) ⎟
D = ⎝ ∂xj ∂xj ⎠ , (70)
∂ ∂U
∂xj (νe h ∂xjy )
(a) Original triangular mesh.
where νe is the effective viscosity (νe = ν + νt ). The vectors
Fx and Fy account for the convective flux. The vector S ac-
counts for the bed slope, G for the bed friction, and D for the
turbulent and viscous diffusion. An explicit discretisation of
(70) is given by:

wn+1 − wn ∂Fx n ∂Fy n

+ (w ) + (w ) = Sn + Gn + Dn , (71)
t ∂x ∂y
(b) Control volume Ci (defined by dashed region).
where the super-index n refers to time t n . The extension of
Fig. 1 Control volume generation from a triangular mesh the scheme to second order in time can be easily imple-
mented for example by a two-step method [14, 20]. Since
in Sects. 3 and 4. A comprehensive description of finite vol- the solver is explicit in time, the time step is limited by the
ume methods for fluid dynamics can be found in [53, 87]. CFL condition, which has been implemented as:
The solver is up to second order accurate in space and time,  
di
and uses upwind schemes for the convective flux and bed ti = CFL min √ , (72)
|U|i + ghi
slope source term. All the terms are treated explicitly in time
except the viscous diagonal, which is treated implicitly. where di = Ai /Pi is the ratio between the area Ai and the
perimeter Pi on each cell Ci .
5.1 Discretisation of the Spatial Domain Integration of (71) over a control volume Ci and yields:

In order to discretise the spatial domain, the solver uses wn+1 − wni 
edge-type finite volumes, which were introduced by Bermú-
i
Ai + (Fx ñx + Fy ñy ) dL
t
j ∈Ki Lij
dez et al. in [9]. These kind of control volumes are generated

from a previous triangulation of the numerical domain, as
= (Sn + Gn + Dn ) dA, (73)
sketched in Fig. 1. The nodes (Ni and Nj ) are placed at the Ci
midpoint of the edges of the triangulation (Fig. 1(a)), and
the cells Ci are built from the vertex (V) and the barycen- where wni is the mean value of w in the cell Ci at time t n , Lij
tre (B) of the triangles (Fig. 1(b)). The control volumes built is the common face to the volumes Ci and Cj , ñ = (ñx , ñy )
in such a way have 4 faces except in the boundaries, where is the unit normal vector to the cell face, and Ki accounts
they have 3 faces. This approach permits an easy definition for all the cells Cj which share any face with the cell Ci .
of the normal vector at the boundary faces, avoiding the in-
determination which may appear in certain boundary nodes 5.3 Convective Flux
when using vertex-type control volumes in complex geome-
tries [32]. The convective flux in (73) is discretised with a Roe up-
wind scheme [80]. Harten’s regularisation [43] is applied
5.2 Discretisation of the Mean Flow Equations when necessary. In order to obtain second order accuracy
in space, the left and right states of the Riemann problem
The 2D-SWE can be written in vectorial form as: are obtained after a linear reconstruction of the variables w
∂w ∂Fx ∂Fy from the nodes to the cell faces, which can be computed as:
+ + = S + G + D,
∂t ∂x ∂y 1 1
⎛ ⎞ ⎛ ⎞ ⎛ ⎞ wIj = wi + ∗i , wiJ = wj + ∗j , (74)
h qx qy 2 2
⎜ 2 2 ⎟ ⎜ qx qy ⎟
w = ⎝ qx ⎠ , Fx = ⎝ qx + gh ⎠ , Fy = ⎝ h ⎠, where ∗i , ∗j are the limited slopes [88] at the nodes Ni
h q q 2 qy2 2
gh
qy x y
h + h 2 and Nj , wIj is the extrapolated value of wi to the cell face
Depth Averaged Modelling of Turbulent Shallow Water Flow with Wet-Dry Fronts 319
⎛ ⎞
|λ̃1 | 0 0
|D| = ⎝ 0 |λ̃2 | 0 ⎠,
0 0 |λ̃3 |

λ̃1 = nx Ũx + ny Ũy , λ̃2 = λ̃1 + c̃ n2x + n2y ,

λ̃3 = λ̃1 − c̃ n2x + n2y , (77)

where wL = wIj and wR = wiJ . In (76–77) nij is a normal

vector to the cell face Lij with the same length as the cell
(a) Triangles used to compute the gradients ∇wi and ∇wj . face. The matrices X, X−1 and |D| are evaluated at the mean
state of Roe [80], which is given by:

 hL + hR
h̃ = hL hR , c̃ = g ,
2
√ √
hL Ux,L + hR Ux,R
Ũx = √ √ , (78)
hL + hR
(b) Linear reconstruction from nodes to faces. 1D view. √ √
hL Uy,L + hR Uy,R
Fig. 2 Second order scheme. Reconstruction of the conservative vari- Ũy = √ √ .
hL + hR
ables from the cell nodes to the cell faces
5.4 Bed Slope Source Term

Lij , and wiJ is the extrapolated value of wj to the cell face When the upwind scheme given by (74–78) is used with a
Lij (Fig. 2(b)). The limited slopes can be computed as: centred discretisation of the bed slope source term S, spuri-
 ous oscillations are generated under hydrostatic flow con-
max[0, min(∇wi rij , ij )], if ij > 0, ditions. In order to avoid this unphysical oscillations, the
∗i =
min[0, max(∇wi rij , ij )], if ij < 0 bed slope source term must be discretised with an upwind
scheme [9, 38, 39, 96]. A suitable upwind discretisation of
with ij = wj − wi (75) S which is free of spurious oscillations when used with the
first order scheme of Roe is given by [9, 38]:
with an analogous expression for ∗j . In (75) rij is the dis-
tance vector between the two nodes Ni and Nj . The gradi- 1  d⊥,ij
Si = SC
i − |nij ||Q|ij Q−1
ij S̃ij ,
ents ∇wi and ∇wj are computed from the values of w at the Ai 2
j ∈Ki
nodes (Ni , Ni1 , Ni2 ) and (Nj , Nj1 , Nj2 ) respectively [95]
1  d⊥,ij
i =
SC |nij |S̃ij ,
(Fig. 2(a)). The limited slopes computed from (75) repro- (79)
duce the Minmod limiter [88]. Ai 2
j ∈Ki
The boundary integral of the convective flux in (73) is ⎛ ⎞
0
approximated by the numerical flux φ ij as: hi + hj zb,j − zb,i
S̃ij = −g ⎝ ñx,ij ⎠ ,
 2 d⊥,ij
ñy,ij
(Fx ñx + Fy ñy ) dL ≈ φ ij (wL , wR , nij ) (76)
Lij where Si is the upwind discretisation of the bed slope source
term, SC i is a centred discretisation of the source term in the
with cell Ci , S̃ij is a centred approximation of the source term at
the cell face Lij , Ai is the area of the cell Ci , and d⊥,ij =
Z(wL , nij ) + Z(wR , nij )
φ ij (wL , wR , nij ) = rij ñij is the projection of the distance between the nodes Ni
2
and Nj (dij ) on to the unit normal vector ñij .
1 In order to prove the convenience of upwinding the bed
− |Q(wL , wR , nij )|(wR − wL ),
2 slope source term, the first order scheme of Roe will be ap-
Z = Fx nx + Fy ny , |Q| = X|D|X−1 , plied to hydrostatic flow conditions. In the hydrostatic as-
⎛ ⎞ sumption, the discrete stationary equations reduce to:
0 1 1

X = ⎝ −ñy Ũx + c̃ñx Ũx − c̃ñx ⎠ , φ ij = Si Ai . (80)
ñx Ũy + c̃ñy Ũy − c̃ñy j ∈Ki
320 L. Cea et al.

The numerical flux given by the scheme for the cell Ci is The discretisation of the bed slope term is obtained intro-
given by: ducing expressions (83) into (79), which yields:
⎛ ⎞
Zi + Zj 1  0
φ ij = − Xij |D|ij X−1
ij (wj − wi ), (81) |nij | hi + hj
2 2 Si Ai = − g (zb,j − zb,i ) ⎝ ñx,ij ⎠
2 2
j ∈Ki ñy,ij
where the subindices ij indicate Roe’s mean state, i.e. Xij = ⎛ ⎞
X(w̃ij ) and |D|ij = |D|(w̃ij ). Considering that in hydrosta-  |nij | 1
tic conditions qx = qy = 0 in all cells, the mean state of Roe + c̃ij (zb,j − zb,i ) ⎝ 0 ⎠ , (85)
2
is reduced to: j ∈Ki 0

 hi + hj where the first addend accounts for the centred contribution
h̃ = hi hj , c̃ = g , and the second addend for the upwind contribution.
2
(82) Considering that in hydrostatic flow the water surface el-
q̃x = 0, q̃y = 0. evation is constant (zb,i + hi = zb,j + hj ), it is straightfor-
ward to show that the upwind contribution in the convective
The values of the different vectors and matrices in (81) flux (second addend in (84)) and in the source term (second
are given by: addend in (85)) are equal. The balance of the centred contri-
⎛ ⎞ butions follows directly considering that for any closed vol-
0 1 1
ume the following property (see (86)) applies for both the x
Xij = ⎝ −ñy,ij c̃ij ñx,ij −c̃ij ñx,ij ⎠ ,
and the y component of the normal vector:
ñx,ij c̃ij ñy,ij −c̃ij ñy,ij
⎛ ⎞  
0 −2c̃ij ñy,ij 2c̃ij ñx,ij h2i ñx,ij |nij | = h2i nx,ij = 0. (86)
1 ⎝
−1
Xij = c̃ij ñx,ij ñy,ij ⎠ , j ∈Ki j ∈Ki
2c̃ij
c̃ij −ñx,ij −ñy,ij
⎛ ⎞ The previous argument is valid for the first order upwind
0 0 0 scheme. If the second order extension of Roe scheme is
Dij = c̃ij |nij | ⎝ 0 1 0 ⎠ , used, the exact balance between the numerical flux and the
0 0 −1 bed slope source term is broken, even if a second order ex-
⎛ ⎞ (83) trapolation is used to evaluate the bed elevation at the cell
0 0 0
|D|ij = c̃ij |nij | ⎝ 0 1 0 ⎠ , faces. This is because the extrapolated value of hi is differ-
0 0 1 ent at each cell face, and therefore,
⎛ ⎞ 
0 0 0 h2Ij ñx,ij |nij | = 0, (87)
1
D−1
ij =
⎝0 1 0 ⎠, j ∈Ki
c̃ij |nij |
0 0 −1
⎛ ⎞ where hIj is the extrapolated value of the flux from the node
0 ⎛ ⎞
⎜ h2i ⎟ hi Ni to the cell face Lij . For this reason, the centred con-
Zi = |nij | ⎜⎝
g 2 ñx,ij ⎟ ,
⎠ w i = ⎝ 0 ⎠. tributions in the flux term discretisation (see (84)) and in
h2i 0 the source term discretisation (see (85)) do not balance any
g 2 ñy,ij
more. There is actually a deficit in the source term given by:
Introducing expressions (83) into (81) gives, after some ⎛ ⎞
 0
mathematical manipulation, the total flux for the cell Ci as: |nij | 2
Ri = −g h ⎝ ñx,ij ⎠ (88)
⎛ ⎞ 2 Ij
j ∈Ki ñy,ij
  |nij | h2i + h2j 0
φ ij = g ⎝ ñx,ij ⎠
2 2 and the following relation applies:
j ∈Ki j ∈Ki ñy,ij
⎛ ⎞ 
 |nij | 1 φ ij = Si Ai + Ri . (89)
− c̃ij (hj − hi ) ⎝ 0 ⎠ , (84) j ∈Ki
2
j ∈Ki 0
We have found that a simple and efficient solution to this
where the first addend accounts for the centred contribution problem is to use a first order approximation for the water
and the second addend for the upwind contribution. depth (hIj = hi ), and a second order approximation for the
Depth Averaged Modelling of Turbulent Shallow Water Flow with Wet-Dry Fronts 321

unit discharges. In this case the deficit Ri is zero:

⎛ ⎞ ⎛ ⎞
 0 0
|nij | 2
Ri = −g hi ⎝ ñx,ij ⎠ = ⎝ 0 ⎠ (90)
2
j ∈Ki ñy,ij 0

and the exact balance between flux and source is recovered.

This approach leads to an hybrid scheme which is able to
compute the hydrostatic flow solution exactly, and at the
same time reduces in a very significant way the numerical Fig. 3 Discretisation of the diffusion term
diffusion.
Hubbard and García-Navarro [47] proposed a high order
correction of the source term in order to obtain a fully sec- where S∗i is the new discretisation of the source term. In
(93), the upper-case subindices refer to extrapolated vari-
ond order scheme which gives an exact balance between
ables at the cell faces, and the lower-case subindices refer
convective flux and bed slope in the hydrostatic case. Fol-
to node values. It is straightforward to show that the source
lowing the ideas presented in [47], when using a second or-
term S∗i balances exactly the convective flux in the hydrosta-
der scheme (88) can be replaced by the more general rela-
tic case.
tion:
All of the applications presented in this paper have been
 computed with both the hybrid scheme (first order for the
Ri = − [Fx,Ij nx,ij + Fy,Ij ny,ij ] water depth and second order for the unit discharges) as
j ∈Ki well as with the fully second order scheme. No significant
 differences have been found in the results given by both
=− [(Fx,Ij − Fx,i )nx,ij + (Fy,Ij − Fy,i )ny,ij ]. schemes, probably because the water depth gradients are
j ∈Ki
not large enough. On the other hand, the hybrid scheme has
(91) been found to be much more stable than the fully second or-
der scheme, especially when the bathymetry is irregular and
Equation (91) is actually a discretisation of the flux gra- there are unsteady wet-dry fronts in the solution.
dient within the node Ni and the face Lij . Hence, it can be
rewritten as: 5.5 Turbulent and Viscous Diffusion
  ∂Fx ∂Fy

d⊥,ij The diffusion term is discretised with a semi-implicit cen-
Ri = − + |nij |
∂x ∂y iI 2 tred scheme. In the following we will refer only to the
j ∈Ki
x-momentum equation, with analogous expressions for the
 d⊥,ij y-momentum equation. Integration of the diffusion term in
=− SiI |nij | , (92)
2 the x-momentum equation over the cell Ci gives:
j ∈Ki
    
∂ ∂Ux ∂ ∂Ux
where it has been considered that we want to maintain a bal- Dtot = νe h + νe h dA
Ci ∂x ∂x ∂y ∂y
ance between the convective flux and the bed slope source  
term in the hydrostatic case. The vector SiI in (92) accounts  ∂Ux ∂Ux
≈ νe,ij hij nx + ny . (94)
for a discretisation of the bed slope term between the node ∂x ∂y ij
j ∈Ki
Ni and the face Lij . Having this in mind, and considering
(89), the discretisation of the source term applied to the fi- When the eddy viscosity is large, it is convenient to im-
nite volume defined in Fig. 1 is given by: plicit the main diagonal of the diffusion source term in order
to relax the stability condition over the time step. In order
1 1  d⊥,ij to do so, following the ideas of Davidson [29], the total dif-
S∗i = Si + R i = Si − SiI |nij |
Ai Ai 2 fusive flux (Dtot ) is split in two parts: an orthogonal dif-
j ∈Ki
fusion (D⊥ ) and a non-orthogonal diffusion (D ), so that
1  d⊥,ij hi + hIj zb,Ij − zb,i

Dtot = D⊥ + D .
= Si − −|nij | g 2
Ai 2 2 d⊥,ij The gradient of the velocity at the cell face Lij is com-
j ∈Ki
⎛ ⎞ puted applying the Gauss theorem on the volume Aij , which
0 is defined by the shaded area in Fig. 3. The two edges of the
× ⎝ ñx,ij ⎠ , (93) volume Aij which pass through the nodes Ni and Nj are de-
ñy,ij fined by the same normal vector as the cell face Lij , i.e. nij .
322 L. Cea et al.
 νt

The other two edges are parallel to the line which joins the ∂
∂xj ((ν + ∂k
σk )h ∂xj )
H1 = νt ,
nodes Ni and Nj , and are defined by the normal vector α ij , ∂
∂xj ((ν + ∂ε
σε )h ∂xj )
with modulus |α ij | = |rij |, so that α ij rij = 0. The area of  
the volume is given by Aij = |nij |d⊥,ij = rij nij . With these min(2νt Sij Sij h, 10εh)
H2 =
definitions, the gradient of Ux at the cell face Lij is evalu- c1ε kε 2νt Sij Sij h
ated as:    
   min(ck u3f , 10εh) −εh
∂Ux  1 ∂Ux 1 H3 = u4 , H4 = . (98)
≈ dA = Ux ñx dL 2
−c2ε εk h
∂x ij Aij Aij ∂x Aij L cε hf

1 The source terms Hm (m = 1, 4) account respectively for

≈ (Ux,j nx,ij + Ux,B αx,ij − Ux,i nx,ij the viscous and turbulent diffusion (H1 ), the production due
Aij
to horizontal velocity gradients (H2 ), the production due to
− Ux,V αx,ij ) (95)
bed friction (H3 ), and the dissipation rate (H4 ). The same
∂Ux
with an analogous expression for ∂y |ij .
The following ex- second order discretisation scheme described for the mean
pression is obtained for the discrete diffusive flux at the flow equations is used for the k–ε equations. However, in
cell Ci : this case the normal convective flux to the face Lij (Z,ij )
 |nij |
Dtot ≈ νe,ij hij (Ux,j − Ux,i ) depends linearly on the depth averaged velocity as:
d⊥,ij
j ∈Ki
   Z,ij = (Ux ñx + Uy ñy )ij ij = Un,ij ij . (99)
D⊥ ≡orthogonal
At each cell face the numerical flux Z∗,ij is computed
 dij as:
+ νe,ij hij (Ux,B − Ux,V )(α̃x,ij ñx,ij + α̃y,ij ñy,ij ) .
d⊥,ij Ij + iJ 1
j ∈Ki
   Z∗,ij = Un,ij − |Un,ij |(iJ − Ij ), (100)
2 2
D ≡non-orthogonal
(96) where Un,ij is a centred discretisation of the normal veloc-
ity to the cell face, and iJ , Ij are linearly extrapolated
The non-orthogonal part (D|| ) is treated explicitly with
values from the cell nodes to the cell faces using (74).
the rest of source terms, while the orthogonal part (D⊥ ) is
All the source terms are discretised at the cell nodes using
split as:
a centred scheme. In order to reinforce the stability of the
|nij |
D⊥,ij = νe,ij hij (Ux,j − Ux,i ) scheme and to help the turbulent quantities k and ε to remain
d⊥,ij
positive during the computation, the following semi-implicit
D⊥,ij linearization of the source terms is used in the solver [29]:
= D⊥,ij Ux,j − qx,i , (97)
hi
H = HnN n+1 + HnP . (101)
ν h |n |
where D⊥,ij = e,ijd⊥,ijij ij
is the orthogonal diffusion coef-
All the negative source terms are discretised in the form
ficient. In (97) all the variables are evaluated at time t n ex-
HnN n+1 , while the positive source terms are included in
cept the unit discharge qx,i , which is evaluated at time t n+1 .
In this way no additional computational cost is introduced, HnP . The production source terms H2 and H3 are always
since there is no need to solve any system of equations. In positive and therefore, they are included in HnP . The dissi-
orthogonal meshes the vectors nij and α ij are perpendicular pation source term H4 is always negative, so it is discretised
and therefore, the non-orthogonal diffusion in (96) vanishes as:
   
(D = 0). −εh −( kε )n (kh)n+1
H4 = 2 = . (102)
−c2ε εk h −c2ε ( kε )n (εh)n+1
5.6 Discretisation of the k–ε Equations

The modelled k–ε equations are written in vectorial form as: The diffusion term H1 can be positive or negative and
thus, it is included in HnP or HnN n+1 depending on its sign.
∂ ∂F,x ∂F,y 
4
+ + = Hm , With these considerations the terms HnN and HnP in (101) are
∂t ∂x ∂y given by:
m=1
   
hk hkUx  n  n
= , F,x = = Ux ,
hε hεUx H1 − kε
HnN = min ,0 + ,
   −c2ε kε
hkUy (103)
F,y = = Uy ,
hεUy HnP = max(H1 , 0)n + Hn2 + Hn3 .
Depth Averaged Modelling of Turbulent Shallow Water Flow with Wet-Dry Fronts 323

5.7 Wet-Dry Fronts

A fixed finite volume mesh is used to discretise the whole

spatial domain, and the control volumes are allowed to wet
and dry during the simulation in order to model unsteady
wet-dry fronts. A wet-dry tolerance parameter (εwd ) is de- (a) Bed redefinition zb,ij = hi . Reflection condition qn,ij = 0.
fined, such that if the water depth in a cell is lower than εwd
the cell is considered to be dry. In the same way, if the water
depth at the face Lij is lower than εwd , the face is consid-
ered to be dry and it does not participate in the calculation.
The water depth is never forced to be zero, in order to keep
the mass conservation property of the scheme. The lowest
(b) No redefinition. No reflection condition.
value of εwd is desired in order to obtain accurate solutions.
However, an excessively low value of εwd promotes numeri- Fig. 4 Wet-dry front
cal instabilities and necessitates the use of a very small CFL,
especially when dealing with a very irregular bathymetry.
The reflection condition sets to zero the normal unit dis-
The wet-dry tolerance parameter is also used in the tur-
charge at the cell face where the wet-dry front occurs (qn,ij ),
bulence models. In the dry cells (hi ≤ εwd ) the eddy viscos-
which is the kinetic condition at a wall boundary:
ity is set to zero, and the turbulent energy and dissipation
are fixed to a residual value (εk and εeps ), which is several qn,ij = qx,ij ñx,ij + qy,ij ñy,ij = 0. (105)
orders of magnitude smaller than the minimum values ex-
pected in the flow. At the same time the production terms Condition (105) only applies when hi < zb,j − zb,i
are set to zero (Pkv = Pk = Pεv = Pε = 0). When a cell is (Fig. 4(a)). Setting the normal unit discharge to zero at the
wet (hi > εwd ) but the water depth is very small, the verti- interface is justified by the assumption of a piecewise con-
cal production term Pkv might become very large. In those stant bed elevation, which is similar to a representation of
situations the production limiter given by (39) helps to con- the bed by a set of small vertical walls. In this way the
trol the production of turbulence, giving more stable results. wet-dry front is only allowed to advance when the water
In such a way, smaller values of εwd can be used without depth in the wet cell is larger than the bed step between cells
producing instabilities in the k–ε model. (hi > zb,j − zb,i , Fig. 4(b)). It should be noticed that the
The fixed mesh approach needs a suitable wet-dry con- unit discharge is not set to zero in the left cell Ci , but only
dition at the fluid interface which ensures the conservation at the face Lij when computing the convective flux. Condi-
of mass and momentum, and is at the same time not diffu- tion (105) also assures that the convective transport of k and
sive and free of spurious oscillations. Assuming a piecewise ε at the wet-dry front is zero when hi < zb,j − zb,i . In this
constant distribution of the bed elevation, two conditions are case the diffusive flux of k and ε over the face Lij is also set
imposed at the wet-dry front: (1) redefinition of the bed ele- to zero.
vation; (2) reflection condition.
The aim of redefining the bed elevation is to obtain an
exact balance at the wet-dry front between the bed slope and 6 Some Applications and Results
the hydrostatic pressure term for hydrostatic conditions. As
detailed by Brufau in [14], if the bed slope is not redefined 6.1 Open Channel Flow in a 90° Bend
spurious waves are generated at the front. If the wet-dry
front occurs between the cells Ci and Cj (assuming with- Shallow flow in open channels is one potential application
out loss of generality that the cell Ci is the wet one), the of the 2D-SWE. In the absence of abrupt obstacles which
modified bed slope at the front (zb,ij ) is defined as [14]: might induce 3D flow features, the flow is usually well
 described by the shallow water approximation. However,
hi − hj if,
zb,ij = strong changes in the direction of the channel can produce
zb,j − zb,i otherwise, 2D recirculation eddies, as well as 3D flow patterns which
might deteriorate the results given by the 2D-SWE. In this
hj ≤ εwd and hi < zb,j − zb,i . (104)
section the free surface flow in an open channel with a 90°
The treatment of the wet-dry fronts given by (104) gives bend is computed and compared with experimental velocity
the exact hydrostatic flow solution for any bed elevation data obtained by Bonillo [11].
without diffusing the front, and without generating any spu- The open channel has two rectilinear sections joined by
rious oscillations of the free surface. a 90° bend (Fig. 5). The first section is 0.86 m wide with a
324 L. Cea et al.

(a) First order scheme.

Fig. 5 Spatial domain in the 90° bend

flat bed. At the end of the first section, just before the bend,
there is a small step with a change in the bed elevation of
zb = −0.013 m. The second section is 0.72 m wide with a
flat bed.
The water discharge in the channel is 29.5 l/s, the mean
water depth is around 0.175 m, and the horizontal velocities
are of the order of 0.20 m/s. The ratio between water depth
and horizontal length scale is approximately H /L ≈ 0.2,
and the Froude number FR ≈ 0.15. The wake stability para- (b) Second order scheme.
meter for this flow conditions can be evaluated using Man-
ning formula as: Fig. 6 Influence of the numerical scheme on the velocity field. Longi-
tudinal velocity Vx (m/s) field
L gn2 L
S = 2cf = 2 0.33 ≈ 0.035, (106)
H h h length model, the k–ε model of Rastogi and Rodi, and the
where g is the gravity acceleration, n is the Manning coef- depth averaged algebraic stress model presented in Sect. 4.
ficient, h is the water depth, and L is the horizontal length All the models have been used with the original constants,
scale, taken to be equal to the width of the channel. Although without any modification or previous calibration.
it is out of the scope of this work, it is worth to mention It is very important to use an accurate numerical scheme
that the wake stability parameter could be used in order to with a low numerical diffusion. Otherwise the solution may
characterise the properties of the recirculation eddy that is be too diffusive, and the velocity profiles too flattened. In
formed just downstream the bend, in the inner part of it. Af- order to check the influence of the numerical scheme on the
ter the results of Chen and Jirka [24, 25] and Lloyd et al. [55, solution, two different numerical schemes have been used:
57], it is expected that the length of the recirculation region, the first order Roe scheme and its second order extension.
its stability and steadiness will depend on the value of S. In The velocity fields shown in Fig. 6 have been obtained with-
any case, in order to verify the former assertion a compre- out using any turbulence model. Further, the effective vis-
hensive experimental campaign should be undertaken. cosity was set to zero, and therefore, all the diffusion comes
The most challenging point of the numerical simulation is from the numerical scheme. The first order scheme is too
the correct prediction of the separated region which appears diffusive, and thus, it is not able to generate a recirculation
in the inner corner of the bend. It is well known that RANS region downstream the bend (Fig. 6(a)). On the other hand,
models usually fail to give accurate predictions of strong the second order scheme, with a much lower numerical dif-
separated flows over 3D geometries. However, in this case fusion, produces a large recirculation bubble (Fig. 6(b)). The
the 2D forcing produced by the shallowness of the flow is less diffusive the numerical scheme is, the larger the recir-
expected to reduce the 3D features of the flow field, helping culation bubble is. Since the water surface elevation and the
to improve the accuracy of the results. In order to correctly bed are almost flat, in this case the hybrid first/second or-
predict the recirculation eddy, three depth averaged turbu- der and the fully second order schemes produce exactly the
lence models have been used in the simulations: the mixing same results.
Depth Averaged Modelling of Turbulent Shallow Water Flow with Wet-Dry Fronts 325

(a) Mesh 1. (b) Mesh 2.

(c) Mesh 3. (d) Mesh 4.

Fig. 7 Numerical meshes in the near bend region

A mesh convergence analysis has been done in four dif- Table 1 Characteristics of the computational grids for the 90° bend.
ywall : distance from the first inner node to the wall in the recirculation
ferent grids, with different mesh size near the walls and region
in the separated region. The spatial domain covered by the
meshes is shown in Fig. 5. The finite volume mesh is gen- Triangular mesh Finite volume mesh
erated from a triangular mesh by the procedure presented in Vertex Elements Volumes Faces ywall
Sect. 5. Figure 7 shows the near bend region of the four tri-
angular meshes which have been used to build the respective M1 738 1296 2033 3888 0.031 m
finite volume meshes. The characteristics of both the trian- M2 1697 3124 4820 9372 0.018 m
gular and the finite volume grids are shown in Table 1. The M3 2495 4712 7206 14136 0.008 m
main difference between them is the mesh size near the wall M4 3003 5700 8702 17100 0.004 m
and in the separated region. A refined grid is needed in order
to avoid an excessive numerical diffusion in this area and to
be able to resolve the recirculation bubble. The grid size in culation length computed with several meshes and turbu-
the recirculation region differs approximately by a factor 2 lence models. The bubble width can be approximately de-
between each mesh. termined from the experimental data, showing a very good
The mesh convergence has been analysed according to agreement with the numerical results given by the three tur-
the separation length and the velocity profiles in the recircu- bulence models (Table 3).
lation region, since it is there where the largest differences The velocity field does not reach convergence until grid
between the four meshes appear. Table 2 shows the recir- M3. The results are quite insensitive to further mesh refine-
326 L. Cea et al.

ment. No significant differences appear between grids M3 is more sensitive to mesh refinement, appearing still some
and M4 (Fig. 8(b,c)). On the other hand, the ML model differences between the velocity profiles obtained with the
grids M3 and M4 (Fig. 8(a)).
The comparison between the numerical results and the
Table 2 Recirculation length computed with several meshes and tur- experimental data will be made in the recirculation region,
bulence models
since the first section of the channel does not present any in-
ML k–ε DASM

M1 0.38 m 0.43 m 0.35 m Table 3 Recirculation width computed with mesh M4

M2 1.12 m 1.05 m 1.12 m Exp. ML k–ε DASM
M3 1.46 m 1.10 m 1.09 m
M4 1.71 m 1.15 m 1.15 m M4 0.13 m 0.12 m 0.11 m 0.12 m

(a) ML model. (b) k–ε model.

(c) DASM.

Fig. 8 Mesh convergence in the 90° bend. Velocity profile Vx (m/s) at x = −1 m

Depth Averaged Modelling of Turbulent Shallow Water Flow with Wet-Dry Fronts 327

(a) Experimental. (b) ML model.

(c) k–ε model. (d) DASM.

Fig. 9 Velocity field Vx (m/s) in the 90° bend. Several turbulence models and experimental results. White colour accounts for positive Vx

teresting flow features. The following facts should be taken all agreement between the experimental and the numerical
into consideration when comparing the numerical and ex- fields is very satisfactory, as it is shown in Fig. 9. The ML
perimental results: (1) the numerical velocity is the depth model overpredicts the recirculation region (white colour
averaged velocity, while the experimental velocity is the hor- in Fig. 9). This is because the model does not account for
izontal velocity at approximately the mid-point between the transport processes, and therefore, all the turbulence gener-
bed of the channel and the free surface (zexp ≈ h/2); (2) ac- ated in the inner corner of the bend is not convected, nei-
cording to the turbulence level and the number of instanta- ther diffused, downstream the bend, as it is in the DASM
neous samples which have been used to compute the sta- and k–ε models. The size of the recirculation bubble com-
tistical estimators of the mean variables, the variability of puted with the DASM and k–ε models shows a better agree-
the experimental mean velocity in the recirculation region is ment with the experimental one, although it is still slightly
about 10%. larger.
It should be noticed that in uniform channel flow the ve- A more detailed comparison of the experimental and nu-
locity profile is quite homogeneous in the vertical direc- merical fields is done in cross sections downstream the bend.
tion, and the velocity at zexp is approximately equal to the Considering the variability of the experimental data, the
depth averaged velocity. However, near the bend a verti- agreement in the velocity is very satisfactory (Fig. 10). From
cal eddy generates in the flow, which tends to destroy the the cross section at x = −1.5 m (Fig. 10(a)), it seems that the
homogeneity in the vertical direction. This has two main length of the recirculation bubble is overestimated by all the
effects. First, the differences between the depth averaged models, specially by the ML model, which is the one which
velocity and the experimental velocity measured at zexp predicts the largest recirculation region. The cross sections
may increase, depending on the actual vertical profile. Sec- at x = −0.8 m and x = −1.3 m show a small bump in the ex-
ond, a modelling error is introduced in the depth averaged perimental velocity profile near the outer wall, which seems
equations, which neglect the dispersion terms due to the to propagate toward the wall. This bump occurs in the re-
non-homogeneity of the vertical profiles. Despite these 3D gion of maximum vertical velocity, which means that it is
flow features appearing in the near bend region, the over- probably generated by the secondary flow originated in the
328 L. Cea et al.

(a) Section x = −1.5 m. (a) Section x = −1.5 m.

(b) Section x = −1.3 m.

(b) Section x = −1.3 m.

(c) Section x = −0.8 m.

(c) Section x = −0.8 m.
Fig. 11 Turbulent kinetic energy k (m2 /s2 ) at several cross sections in
Fig. 10 Longitudinal velocity Vx (m/s) at several cross sections in the the 90° bend
90° bend

Regarding the performance of the turbulence models, the

bend region. It should be remarked again that, since the ex- differences between the computed velocity fields, which ap-
perimental velocity has been measured at mid-water depth, pear mainly in the recirculation bubble, are really small
it may not be representative of the depth averaged velocity compared with the scatter in the experimental data (Fig. 11).
at points with non-uniform vertical profiles. The ML model gives a slighter poorer prediction of the
Depth Averaged Modelling of Turbulent Shallow Water Flow with Wet-Dry Fronts 329

bubble length, but its degree of accuracy is rather satisfac- shear stress, and therefore its prediction is very dependent
tory, specially considering the simplicity of the model. The on the turbulence model used. The flow separates behind
DASM and k–ε models give very similar results. Larger dif- the inlet slot baffle and the so-called lower eddy appears.
ferences appear in the turbulence intensity near the bend Going through the outlet slot the flow is forced to reattach.
(cross section x = −0.8 m), where the maximum turbulent The correct identification of these eddies is of great impor-
kinetic energy is better predicted by the DASM than by the tance from a practical point of view, since they have a big
k–ε model (Fig. 11(c)). However, the differences between influence guiding the fish along the pool. These are also ar-
models smear out really fast downstream, and no differ- eas where the fish may take a break on its way through the
ence et all can be observed at the cross section x = −1.3 m fishway.
(Fig. 11(b)). Although both models slightly over-predict the The maximum velocity, which occurs just downstream
turbulence intensity inside the bubble, the overall agreement the inlet slot, is another important flow feature which should
and downstream evolution is very satisfactory. be correctly computed. It determines whether or not the
fishes are able to pass from one pool to another. If the slot ve-
6.2 Turbulent Flow in a Vertical Slot Fishway locity is too high the fish may not be able to swim against it.
A correct turbulence modelling is of great importance to
predict all the previously described features. An excessive
Fishways are hydraulic structures which allow the upstream
turbulence level will diffuse too much the velocity profiles,
migration of fishes through engineering constructions and
tending to eliminate the recirculation eddies. On the other
natural obstructions in rivers. A vertical slot fishway is
hand, a low level of turbulent energy will sharpen the pro-
a channel divided into several pools separated by vertical
files too much, predicting excessively high velocities.
slots through which the water flows downstream. The per-
The only empirical parameter which shows up in the
formance of the fishway depends directly on the velocity
mean flow equations is the bed friction coefficient. It turns
field and turbulence intensity in the pools. Therefore, a cor- out that its influence on the mean velocity field is mini-
rect evaluation of these variables has a great importance in mal, mainly because the experimental model walls are very
the design process. In this section the free surface flow in smooth. For the same reason, the turbulence production due
two different designs of vertical slot fishways (Designs T1 to bottom friction is also low compared to the horizontal
and T2) is computed with the solver presented in Sect. 5. shear production. A sensibility analysis of the results to the
The numerical velocity and turbulent fields are compared bed friction coefficient was done using the Manning’s for-
with the experimental results of Puertas et al. [72]. mula. No differences were found for values of the Man-
It is not straightforward that the 2D-SWE are able to re- ning’s coefficient in the range n = 0–0.02 s/m1/3 , which
produce the flow pattern in a vertical slot fishway. In fact, cover widely the possible values for the experimental model.
most of the hypotheses which are made in the derivation
of the 2D-SWE (Sect. 2) are broken in some regions of the 6.2.1 Numerical Meshes
flow, specially in the near slot region. However, as it will be
shown from the results, the numerical and experimental data Each numerical mesh contains an inlet pool, three active
agree quite well. pools with either a T2 or a T1 design, and an outlet pool
From the experimental results of several researchers [69, (Fig. 12). The first pool has a transition role, allowing the
turbulent kinetic energy and velocity fields to develop. The
72–74] it follows that for uniform flow conditions the ve-
second pool has been used in order to compare the numeri-
locity field in vertical slot fishways is almost independent of
cal and experimental flow fields. In all the cases it has been
the total flow discharge. On the other hand, the water depth
checked that there are no significant differences between the
is proportional to the flow discharge with an almost linear re-
flow fields in the second and third pools (Fig. 12). Therefore,
lation. The first feature we should expect from the numerical
a uniform state of the flow in pools 2 and 3 can be assumed.
model is to reproduce this behaviour. In order to prove so,
In order to obtain a mesh independent solution three
three different discharges covering almost the full range of
meshes with different spatial resolution were tested for each
experimental flow conditions were used in the simulations: pool design, which are addressed as T1m0, T1m1, T1m2 for
35 l/s, 65 l/s and 105 l/s. design T1, and T2m0, T2m1, T2m2 for design T2 (Table 4).
The flow pattern in both designs can be described as a Figure 13 shows the resolution of meshes T1m1 and T2m1,
main jet which crosses the pool from the inlet to the outlet which are the ones chosen to perform the definitive simula-
slot (Fig. 12). At each side of the jet a recirculation eddy ap- tions. The mesh convergence analysis was done for the dis-
pears. The size and strength of these eddies depends on the charge of 65 l/s with both the mixing length and k–ε models.
pool design. The numerical model should predict the size of No significant differences were found in the results obtained
the two recirculation regions which occur on both sides of with the intermediate (T1m1, T2m1) and the finest meshes
the main jet. The upper eddy is generated by shear friction, (T1m2, T2m2). All the results presented in this paper have
resembling a cavity flow. The eddy is driven by turbulent been obtained with the meshes T1m1 and T2m1.
330 L. Cea et al.

Depth averaged speed |V | (m/s). Design T1

Turbulent kinetic energy k (m2 /s2 ). Design T2

Fig. 12 Depth averaged speed and turbulent kinetic energy fields. Q = 65 l/s. k–ε model

Table 4 Characteristics of the numerical meshes for the pool de-

signs T1 and T2

Triangular mesh Finite volume mesh

Vertex Elements Volumes Faces

T1m0 1462 2545 4006 7635

T1m1 5610 10334 15943 31002
T1m2 8932 16926 25857 50778
T2m0 3588 6692 10279 20076
(a) Mesh T2m1
T2m1 6182 11716 17897 35148
T2m2 8692 16610 25301 49830

6.2.2 Boundary Conditions

In all the simulations the flow is subcritical at the open

boundaries. At the inlet boundary the total water discharge
is distributed over a uniform discharge profile normal to the
open boundary, while the water depth is computed by the
numerical scheme. Three discharges were used in the sim- (b) Mesh T1m1
ulations: 35 l/s, 65 l/s and 105 l/s. At the outlet boundary
Fig. 13 Numerical meshes for the fishway. Designs T1 and T2
the free surface elevation is fixed for each discharge in order
to obtain approximately the same water depth at the mid-
dle section of all the pools. This is the same procedure by
which the downstream boundary condition was fixed in the as it passes through the first pool. In all the simulations
experimental model. the turbulent variables at the inlet boundary were fixed to
The actual value of the turbulent kinetic energy and dis- kin = 10−2 m2 /s2 and εin = 10−2 m2 /s3 , which gives an
sipation rate at the inlet boundary is difficult to estimate. eddy viscosity of approximately νt ≈ 10−3 m2 /s.
A too high turbulence level is not desirable since it would Wall functions are used at the solid boundaries. The wall
take too much time to dissipate and therefore, it would af- mesh size varies widely with the location due to the different
fect strongly the results in the downstream pools. It is more local flow conditions, but was always kept below an upper
adequate to let the flow produce the real state of turbulence limit of y + ≈ 100.
Depth Averaged Modelling of Turbulent Shallow Water Flow with Wet-Dry Fronts 331
Fig. 14 Numerical and
experimental depth averaged
longitudinal velocity Vx (m/s) at
several cross sections.
Design T2

(a) Q = 65 l/s. x = 0.26 m (b) Q = 65 l/s. x = 0.46 m

(c) Q = 105 l/s. x = 0.26 m (d) Q = 105 l/s. x = 0.46 m

6.2.3 Velocity Field experimental data. The mixing length (ML) model is more
sensitive to the total discharge. This is because the length
A comparison of the depth averaged longitudinal veloc- scale used to compute the eddy viscosity with the ML model
ity at several cross sections reveals a quite satisfactory depends directly on the water depth (see (24)). Therefore,
agreement between the experimental and numerical results with larger depths the turbulent stresses (diffusive forces) in-
(Figs. 14 and 15). The experimental velocity has been com- crease and the velocity profiles become more flattened. The
puted from the results of Puertas et al. [72] as a weight dependence of the turbulent length scale on the water depth
average of the velocity at several elevations (2 elevations is a weakness of the ML model.
for Q = 35 l/s, 3 elevations for Q = 65 l/s and 5 eleva-
tions for Q = 105 l/s). The ASM and k–ε models pre-
6.2.4 Turbulence Field
dict fairly well the width of the recirculation eddies as
well as the evolution of the maximum velocity as the flow
crosses the pool. Although both models give slightly differ- The ASM gives better predictions of the turbulence field
ent results, none of them can be considered in better agree- than the k–ε model, specially regarding the Reynolds
ment with the experimental data. The predictions of the ML stresses. The k–ε model fails to predict the anisotropy be-
model are somewhat worse, specially as the discharge in- tween the two horizontal Reynolds stresses, and gives exces-
creases. sively large values of the transversal component v 2 (Figs. 17
The DASM and k–ε models predict maximum velocity and 18). On the other hand, the agreement between the ASM
and recirculation regions which are quite insensitive to the predictions of v 2 and the experimental data is remarkable,
total discharge. This result is in direct agreement with the specially for the pool design T1.
332 L. Cea et al.

(a) Q = 35 l/s. x = 0.16 m (b) Q = 35 l/s. x = 0.46 m

(c) Q = 105 l/s. x = 0.16 m (d) Q = 105 l/s. x = 0.46 m

Fig. 15 Numerical and experimental depth averaged longitudinal velocity Vx (m/s) at several cross sections. Design T1

Nonetheless, the global agreement between the exper- value of the turbulent kinetic energy given by the k–ε model
imental and numerical normal Reynolds stresses is quite is inside the eddy region, with values around 0.24 m2 /s2 .
satisfactory. The largest differences appear in the main jet With the ASM the peak values are lowered to 0.18 m2 /s2
stream, where the ASM tends to underpredict the turbulent and moved downstream.
energy, while the k–ε model overpredicts it (Fig. 17). The The strong anisotropy between the longitudinal and
excessively large turbulence level given by the k–ε model transversal Reynolds stresses downstream the inlet slot is
just downstream the inlet slot might be explained by the well predicted by the ASM model. While the longitudinal
strong shear strain, with velocity gradients of the order of component (u2 ) is underpredicted, probably because the
20 s−1 . The flow conditions in this region, with a very strong turbulence produced by the air bubbles in the slot region
swirl, as well as separation and reattachment in a very short is not considered in the model, the transversal component
distance, are beyond the capabilities of the eddy viscosity (v 2 ) is quite accurately computed. On the other hand, the
turbulence models. The largest source of turbulence in the k–ε model overpredicts both components in this region, spe-
pool is due to the strong shear strain in the inlet slot and cially the transversal one.
lower eddy. Although the turbulence production in the k–ε 6.3 Tidal Flow in the Crouch Estuary
model is controlled by the production limiter (see (38)) and
the realizability condition (see (40)), the ASM still gives The Crouch estuary, Essex, UK (Fig. 19), is a challenging
lower values for the turbulent stresses in the slot stagnation test case if we consider its complex bathymetry, with ex-
regions and in the lower recirculation eddy. The maximum tensive intertidal flats. Tidal length is approximately 25 Km
Depth Averaged Modelling of Turbulent Shallow Water Flow with Wet-Dry Fronts 333

(a) Design T1. ML model. (b) Design T2. ML model.

(c) Design T1. k–ε model. (d) Design T2. k–ε model.

(e) Design T1. ASM. (f) Design T2. ASM.

Fig. 16 Depth averaged velocity field. Q = 105 l/s. The black line separates the regions with positive and negative longitudinal velocity

along its largest arm. The tidal range is approximately 4.8 m The flow in the model is completely driven by the wa-
at the mouth. The maximum depth is about 15 m at high ter surface elevation at the mouth, which is the only open
tides. Freshwater influx is negligible in the whole estuary, boundary condition to be imposed. The rest of the bound-
having a mean flow of only 0.4 m3 s−1 , and therefore, all aries are treated as walls. Due to the size and complexity of
the water in the estuary is salt water. The large separation the computational domain (27.65 Km2 ) the mesh size near
between horizontal and vertical length scales justify the use the walls is relatively coarse and therefore, a slip condition
of a 2D-SWE model. has been used at the wall boundaries. Accordingly, the nor-
334 L. Cea et al.

(a) Q = 65 l/s. x = 0.26 m (a) Q = 65 l/s. x = 0.16 m

(b) Q = 65 l/s. x = 0.46 m

(b) Q = 65 l/s. x = 0.46 m

(c) Q = 65 l/s. x = 0.86 m

(c) Q = 65 l/s. x = 0.86 m
Fig. 18 Numerical and experimental depth averaged horizontal
Fig. 17 Numerical and experimental depth averaged horizontal Reynolds stresses u2 and v 2 (m2 /s2 ) at several cross sections. Design
Reynolds stresses u2 and v 2 (m2 /s2 ) at several cross sections. Design T1. Q = 65 l/s
T2. Q = 65 l/s
they do not participate in the solution. On the other hand
the wet-dry front, which defines the fluid extension, plays a
mal velocity as well as the diffusion of k and ε are set to more important role in the flow field.
zero at the boundary faces. Nonetheless, most of the time The water surface elevation at the mouth is assumed to
during the simulation the wall boundaries are dry and thus, be constant all along the boundary. Its level is directly ob-
Depth Averaged Modelling of Turbulent Shallow Water Flow with Wet-Dry Fronts 335

Fig. 19 Bathymetry of the Crouch-Roach estuary for the numerical model, relative to mean sea level at the mouth. Locations of tide gauges and
current meter deployments are shown

tained from a tidal gauge. The flow at the mouth is always Table 5 Characteristics of the finite volume mesh and its associated
triangular mesh
subcritical. Hence, only the water depth needs to be imposed
during the ebb tide (subcritical outlet). During the flood tide Triangular mesh Finite volume mesh
(subcritical inlet) the flow is forced to be perpendicular to Vertices Elements Volumes Faces
the boundary, as an additional condition to the imposed wa-
ter depth. Some regions of the mouth boundary may become Mesh 1 17258 31733 48995 95199
dry or wet depending on the water surface level. For this rea-
son special care must be taken at the boundary line, because
it may be the case that the imposed water depth is too small, Approximately 17 tidal cycles (9 days) were simulated,
and the flow becomes supercritical. This is likely to occur and the results were compared with experimental water level
when the flow enters the domain and the bed slope at the and current speed time series. In order to check the sensi-
boundary is large. This situation can be avoided by limiting tivity of the numerical simulation to the wet-dry tolerance
the maximum Froude number at the inlet boundary, ensuring parameter εwd , three different values were tried in the com-
that the flow is always subcritical. Accordingly, the normal putations, with particular reference to the performance near
unit discharge at the inlet boundary is limited as: Bridgemarsh island. The smallest value of the wet-dry tol-
 erance parameter εwd is desirable from an accuracy point of
qn,in ≤ 0.95 gh3in , (107) view. But too small values of εwd promote numerical insta-
bilities, especially in domains with a very irregular bathym-
where qn,in is the normal unit discharge at the inlet, and hin etry, and require the use of very small time steps. A sensitiv-
is the imposed water depth. Condition (107) only needs to be ity analysis using several values of εwd (10−2 m, 10−3 m,
applied to a very few boundary cells, and at very few time 10−4 m) revealed no differences in the results for values
steps during the simulation. of εwd < 10−3 m (Fig. 21). This value was therefore used
A Manning coefficient of n = 0.02 ms−1/3 was used in in all the other simulations.
the inter-tidal and sub-tidal regions, where the bottom sed- The time step used in the present simulations is rather
iments are predominantly muddy. There are also several small (around 1 s) and thus, the results obtained with the
marsh regions, which are assigned a Manning coefficient first order and the second order discretisation in time are
n = 0.05 ms−1/3 . Nevertheless, the sensitivity of the numer- identical to all practical effects. The hybrid scheme was used
ical solution to bed friction is low with regard to both water in the computations, which proved to be much more stable
depth and tidal current velocities. than the fully second order scheme, while producing almost
After using several meshes with different resolution, a identical results.
mesh with approximately 50000 volumes was chosen. The Both the ML and k–ε models give very similar eddy vis-
numerical mesh used (Table 5) covers the whole estuary us- cosity fields (Fig. 22). This is because turbulence is mainly
ing 48995 control volumes. The area of the control volumes produced by bed friction, and both models account in a
in the main channel is about 300 m2 . similar way for this process. For uniform channel flow the
336 L. Cea et al.
Fig. 20 Detail of Mesh 1 in the
vicinity of Bridgemarsh Island

(a) Wallasea. (b) Fambridge.

Fig. 21 Dependence of the horizontal velocity on the wet-dry tolerance parameter (εwd ). νt = 0. The differences between the results obtained
with εwd = 0.001 m and with εwd = 0.0001 m cannot be appreciated in the plots

eddy viscosity given by the k–ε model reduces to νtk –ε = bed friction, by Davies et al. [31] when computing the tidal
0.08uf h, while the ML model gives νtML = 0.068uf h (see flow in the Irish Sea, and by Lloyd and Stansby [55] when
Sect. 3). For this reason the velocity and water depth fields modelling the flow around conical islands.
are very insensitive to the turbulence model used. The velocity field in the whole estuary is very dependent
Even in the mouth of the estuary, where the turbulence on the bathymetry, with velocities being highest in the deep-
intensities are highest (Fig. 22), the turbulent horizontal est regions of the estuary. The sensitivity of the numerical
Reynolds number is rather large (Rt = U L/νt ≈ 21000). results to the bed friction coefficient and turbulence model
This means that the turbulent diffusion forces are small is much smaller than the differences between numerical and
when compared to the convective forces, which diminishes experimental data (Figs. 23, 24 and 25). Some of the dis-
the influence of the turbulence model on the mean veloc- agreements between numerical and experimental results are
ity field. Velocity and water depth fields independent of the probably due to errors in the local bathymetry of the model.
turbulence model were also obtained by Babarutsi et al. [6] The fact of comparing the experimental velocity measured
when modelling shallow recirculating flows dominated by at a fixed height of 2 m above the bed, with the numeri-
Depth Averaged Modelling of Turbulent Shallow Water Flow with Wet-Dry Fronts 337

(a) ML model. (b) k–ε model.

Fig. 22 Eddy viscosity field at the mouth of the estuary. Ebb tide. ML and k–ε models. Dry regions are shown in black

(a) Water depth. (b) Velocity.

Fig. 23 Water depth and horizontal velocity time series at Wallasea. ML and k–ε turbulence models, and zero eddy viscosity (νt = 0)

cal depth-averaged velocity, might also account for some of schemes used to solve them, practical applications and ex-
these discrepancies. perimental validation. Special attention has been placed in
Figures 23, 24 and 25 show several numerical and ex- the modelling of turbulent flows and in the treatment of
perimental time-series for water depth and current speed. In wet-dry fronts. A depth averaged algebraic stress model
general water depth is rather well predicted by the model (DASM)has been proposed as an extension of the 2D alge-
(Figs. 23(a), 24(a) and 25(a)) except at Fambridge, where braic stress model, with additional source terms which ac-
the tidal amplitude is slightly underpredicted. Regarding the count for the production of Reynolds stresses due to verti-
current speed, both ebb and flood flows at Wallasea are very cal shear. In addition to the DASM, depth averaged mixing
well predicted. The agreement at Holliwell and Fambridge length and k–ε models have been discussed and used sistem-
is also satisfactory, although the maximum ebb velocity is atically in the practical applications. The computational cost
slightly underpredicted (Figs. 24(b) and 25(b)).
and numerical stability of the DASM are similar to those
of the k–ε model. Additional limiters to the production of
7 Conclusions turbulence proposed by Menter [59] and Durbin [37] have
been introduced in the original models in order to control
In this paper we have presented a comprehensive study on the excessively large turbulence production given by eddy
the numerical modelling of turbulent shallow flows, includ- viscosity models in stagnation regions. Menter’s limiter has
ing a critical description of the equations and numerical been adapted in order to control also the turbulence produc-
338 L. Cea et al.

(a) Water depth. (a) Water depth.

(b) Velocity. (b) Velocity.

Fig. 24 Water depth and horizontal velocity time series at Holliwell. Fig. 25 Water depth and horizontal velocity time series at Fambridge.
ML and k–ε turbulence models, and zero eddy viscosity (νt = 0) ML and k–ε turbulence models, and zero eddy viscosity (νt = 0)

tion due to bed friction near wet-dry fronts. This turns up irregular bathimetry, not diffussive even in the presence of
in an improved stability of the k–ε solver in the presence discontinuities in the bed elevation, and it does not generate
of wet-dry fronts over irregular bathimetry. With this limiter spurious oscillations which could contaminate the numerical
and a semi-implicit linearisation of the source terms in the results.
k–ε equations, a very stable scheme has been obtained. It turned up that, even if the flow is fully turbulent in all
In order to avoid spurious oscillations when the bathime- the applications, the importance of the turbulence modelling
try is irregular, an upwind discretisation of the bed slope is much dependent on the problem, to the point that in the
term is used in the solver [9], with second order corrections Crouch estuary turbulence has not importance at all in the
proposed by Hubbard and García-Navarro [47]. In this way global flow field, while in the fishway its role is critical.
a fully second order upwind scheme free of spurious oscil- Still, it should be remarked that the turbulence field might
lations is obtained. Alternatively to the fully second order be important in order to compute sediment transport or dis-
scheme, a rather simple, free of spurious oscillations and sta- persion of pollutants, which are problems usually dominated
ble scheme, which reduces the numerical diffusion in a sig- by turbulent diffusion.
nificant way, has been proposed by just using a second order The comparison between the numerical and experimental
discretisation for the two unit discharge components, whilst results in the applications presented in Sect. 6 showed that
keeping a first order discretisation for the water depth and the depth averaged shallow water equations, coupled with a
bed elevation. The resulting hybrid scheme has been used in suitable turbulence model, may be used in order to compute
the practical applications, giving accurate and stable results. the free surface flow in vertical slot fishways as well as the
Nevertheless, it is important to remark that the water depth tidal flow in complex estuaries with extensive tidal flats. In
gradient in all the applications was quite smooth. the case of the fishway, even if there are some regions in the
The treatment of wet-dry fronts proposed in Sect. 5 is flow where the modelling assumptions are broken, specially
based on the ideas of Brufau et al. [16], with a slight mod- near the slot region, the numerical results show a satisfactory
ification in the reflection condition. The wet-dry condition global agreement with the experiments. The most important
used in the solver is numerically stable, even with a very features of the flow field (maximum velocity in the pool, size
Depth Averaged Modelling of Turbulent Shallow Water Flow with Wet-Dry Fronts 339

of the recirculation regions, turbulence level, dependence of 10. Bijvelds MDJP, Kranenbutg C, Stelling GS (1999) 3D numerical
the solution on the total discharge . . .) are well predicted by simulation of turbulent shallow-water flow in square harbor. J
Hydraul Eng 125(1):26–31
the solver. This is an interesting finding because from the
11. Bonillo J (2000) Un modelo de transporte de sustancias solubles
nature of the fishway flow (3D turbulence, water depth sim- para flujos turbulentos en lámina libre. Tesis doctoral, Área de
ilar to the horizontal length scale, 3D mean flow features Ingeniería Hidráulica, Universidad de A Coruña
in the slot region, . . . ), it seems that only a 3D model ap- 12. Booij R (1989) Depth-averaged k–ε modelling. In: Proc 23rd
IAHR congr, Delft, The Netherlands, vol A. IAHR, pp 198–206
proach would give successful and accurate results. However,
13. Bradford F, Sanders BF (2002) Finite-volume model for
a depth average model can give very satisfactory predictions shallow-water flooding of arbitrary topography. J Hydraul Eng
with a relatively low computational cost. The importance of 128(3):289–298
a correct turbulence modelling in the fishway was confirmed 14. Brufau P (2000) Simulación bidimensional de flujos hidrod-
inámicos transitorios en geometrías irregulares. Tesis doctoral,
by the numerical results. The ML model failed to predict Área de Mécanica de Fluidos, Universidad de Zaragoza
the velocity field accurately, specially for large water dis- 15. Brufau P, García-Navarro P (2000) Two-dimensional dam break
charges, because the model assumes a turbulent length scale flow simulation. Int J Numer Methods Fluids 33:35–57
proportional to the water depth, which is not correct in the 16. Brufau P, Vázquez-Cendón ME, García-Navarro P (2002) A nu-
merical model for the flooding and drying of irregular domains.
fishway flow. The DASM and k–ε models give fairly ac- Int J Numer Methods Fluids 39(3):247–275
curate velocity fields. The original purpose of using the al- 17. Brufau P, García-Navarro P, Vázquez-Cendón ME (2004) Zero
gebraic stress model was to account for the effects of the mass error using unsteady wetting-drying conditions shallow
anisotropic Reynolds stresses in the flow. When compared to flows over dry of irregular topography. Int J Numer Methods Flu-
ids 45:1047–1082
the k–ε model, which assumes an isotropic eddy viscosity, 18. Casulli V, Stelling GS (1998) Numerical simulation of 3D quasi-
the DASM improves the prediction of the turbulent stresses. hydrostatic free-surface flows. J Hydraul Eng 124(7):678–686
On the other hand, no improvement has been found in the 19. Casulli V, Zanolli P (2002) Semi-implicit numerical modeling
water depth and velocity fields, which are predicted with a of nonhydrostatic free-surface flows for environmental problems.
Math Comput Model 36:1131–1149
similar degree of accuracy with either the DASM or the k–ε 20. Cea L (2005) An unstructured finite volume model for un-
model. Probably the main disagreements between the exper- steady turbulent shallow water flow with wet-dry fronts: numer-
imental and numerical results are due to the approximations ical solver and experimental validation. Doctoral Thesis, Depar-
assumed in the depth averaged equations, and may only be tamento de Métodos Matemáticos y de Representación, Univer-
sidad de A Coruña
improved by using a 3D model rather than implementing 21. Cea L, Ferreiro A, Vázquez-Cendón ME, Puertas J (2004) Ex-
more complex depth averaged turbulence models. perimental and numerical analysis of solitary waves generated
by bed and boundary movements. Int J Numer Methods Fluids
46(8):793–813
22. Cea L, French J, Vázquez-Cendón ME (2006) Numerical mod-
References elling of tidal flows in complex estuaries including turbulence:
an unstructured finite volume solver and experimental validation.
1. Akanbi AA, Katopodes ND (1988) Model for flood propagation Int J Numer Meth Eng 67(13):1909–1932
on initially dry land. J Hydraul Eng 114(7):689–706 23. Cea L, Pena L, Puertas J, Vázquez-Cendón ME, Peña E (2007)
2. Altai S, Zhang J, Chu VH (1999) Shallow turbulent flow simula- Application of several depth-averaged turbulence models to sim-
tion using two-length-scale model. J Eng Mech 125(7):780–788 ulate flow in vertical slot fishways. J Hydraul Eng 133(2):160–
3. Aulisa E, Manservisi S, Scardovelli R (2003) A mixed markers 172
and volume-of-fluid method for the reconstruction and advection 24. Chen D, Jirka GH (1995) Experimental study of plane turbulent
of interfaces in two-phase and free-boundary flows. J Comput wakes in a shallow water layer. Fluid Dyn Res 16:11–41
Phys 188:611–639 25. Chen D, Jirka GH (1997) Absolute and convective instabilities
4. Babarutsi S, Chu VH (1991) A two-length-scale model for quasi- of plane turbulent wakes in a shallow water layer. J Fluid Mech
two-dimensional turbulent shear flows. In: Proc 24th congr of 338:157–172
IAHR, vol C, Madrid, Spain, pp 51–60 26. Chen D, Jirka GH (1998) Linear instability analysis of turbulent
5. Babarutsi S, Chu VH (1998) Modelling transverse mixing layer mixing layers and jets in shallow water layers. J Hydraul Res
in shallow open-channel flows. J Hydraul Eng 124(7):718–727 36(5)
6. Babarutsi S, Nassiri M, Chu VH (1996) Computation of shal- 27. Chu VH, Babarutsi S (1988) Confinement and bed-friction
low recirculating flow dominated by friction. J Hydraul Eng effects in shallow turbulent mixing layers. J Hydraul Eng
122(7):367–372 114(10):1257–1274
7. Bellos CV, Soulis JV, Sakkas JG (1991) Computation of two- 28. Chu VH, Wu JH, Khayat RE (1991) Stability of transverse shear
dimensional dam break induced flows. Adv Water Res 14(1):31– flows in shallow open channels. J Hydraul Eng 117(10):1370–
41 1388
8. Bermúdez A, Vázquez-Cendón ME (1994) Upwind methods for 29. Davidson L (1993) Implementation of a k–ε model and a
hyperbolic conservation laws with source terms. Comput Fluids Reynolds stress model into a multiblock code. Tech Rep CRS4-
23(8):1049–1071 APPMATH-93-21, Applied Mathematics and Simulation Group
9. Bermúdez A, Dervieux A, Desideri JA, Vázquez-Cendón ME CRS4, Cagliary, Italy
(1998) Upwind schemes for the two-dimensional shallow water 30. Davidson L (1997) An introduction to turbulence models. Tech
equations with variable depth using unstructured meshes. Com- Rep 97/2, Dept of Thermo and Fluid Dynamics, Chalmers Uni-
put Methods Appl Mech Eng 155:49–72 versity of Technology
340 L. Cea et al.

31. Davies AM, Jones JE, Xing J (1997) Review of recent devel- 56. Lloyd PM, Stansby PK (1997) Shallow water flow around model
opments in tidal hydrodynamic modeling. II: Turbulence energy conical islands of small side slope. Part 2: Submerged. J Hydraul
models. J Hydraul Eng 123(4):278–292 Eng 123:1068–1077
32. Dervieux A, Desideri JA (1992) Compressible flow solvers using 57. Lloyd PM, Stansby PK, Chen D (2001) Wake formation around
unstructured grids. Rapports de Recherche 1732, INRIA islands in oscillatory laminar shallow water flows. Part 1. Exper-
33. Ding Y, Jia Y, Wang SSY (2004) Identification of Manning’s imental investigation. J Fluid Mech 429:217–238
roughness coefficients in shallow water flows. J Hydraul Eng 58. Maronnier V, Picasso M, Rappaz J (1999) Numerical simulation
130(6):501–510 of free surface flows. J Comput Phys 155:439–455
34. Dodd N (1998) Numerical model of wave run-up, overtopping, 59. Menter FR (1993) Zonal two-equation k–ω turbulence models
and regeneration. J Waterw, Port, Coastal, Ocean Eng 124(2):73– for aerodynamic flows. AIAA paper 93-2906. In 24th fluid dy-
81 namics conference, 6–9 July 1993, Orlando, FL
35. Dracos T, Giger M, Jirka GH (1992) Plane turbulent jets in a 60. Miglio E, Quarteroni A, Saleri F (1999) Finite element approx-
bounded fluid layer. J Fluid Mech 241:587–614 imation of quasi-3D shallow water equations. Comput Methods
36. Duan JG (2004) Simulation of flow and mass dispersion in me- Appl Mech Eng 174:355–369
andering channels. J Hydraul Eng 130(10):964–976 61. Minh-Duc B, Wenka T, Rodi W (2004) Numerical model-
37. Durbin P (1996) On the k–ε stagnation point anomaly. Int J Heat ing of bed deformation in laboratory channels. J Hydraul Eng
Fluid Flow 17:89–90 130(9):894–904
38. Fernández-Nieto ED (2003) Aproximación numérica de leyes 62. Molls T, Chaudhry MH (1995) Depth-averaged open-channel
de conservación hiperbólicas no homogéneas. Aplicación a las flow model. J Hydraul Eng 121(6):453–465
ecuaciones de aguas someras. Tesis doctoral, Universidad de 63. Nichols BD, Hirt CW (1973) Calculating three-dimensional free
Sevilla surface flows in the vicinity of submerged and exposed struc-
39. García-Navarro P, Vázquez-Cendón ME (2000) On numerical tures. J Comput Phys 12:234–246
treatment of the source terms in the shallow water equations. 64. Nichols BD, Hirt CW (1975) Methods for calculating multidi-
Comput Fluids 29:951–979 mensional, transient free surface flows past bodies. In: Proc 1st
40. Gatski T, Speziale C (1993) On explicit algebraic stress models int conf on numer ship hydrodyn, Gaithersburg, MD
for complex turbulent flows. J Fluid Mech 154:59–78 65. Olsen NRB (1999) Two-dimensional numerical modelling of
41. Gerrits J (2001) Dynamics of liquid-filled spacecraft. PhD thesis, flushing processes in water reservoirs. J Hydraul Res 37(1):3–16
University of Groningen, The Netherlands 66. Olsen NRB, Kjellesvig HM (1998) Three-dimensional numerical
42. Giger M, Dracos T, Jirka GH (1991) Entrainment and mixing in flow modelling for estimation of spillway capacity. J Hydraul Res
36(5):775–784
plane turbulent jets in shallow water. J Hydraul Res 29(4):615–
67. Olsen NRB, Melaaen MC (1998) Three-dimensional calculation
643
of scour around cylinders. J Hydraul Eng 119(9):1049–1054
43. Harten A, Lax P, van Leer A (1983) On upstream differencing
68. Peltier E, Duplex J, Latteux B, Pechon P, Chausson P (1991)
and Godunov-type schemes for hyperbolic conservation laws.
Finite element model for bed-load transport and morphological
SIAM Rev 25:35–61
evolution. Comput Model Ocean Eng 91
44. Hirt CW, Nichols BD (1981) Volume of fluid (VOF) method for
69. Pena L, Cea L, Puertas J (2004) Turbulent flow: an experimental
the dynamics of free boundaries. J Comput Phys 39:201–225
analysis in vertical slot fishways. In: Fifth international sympo-
45. Hsieh TY, Yang JC (2003) Investigation on the suitability of two-
sium on ecohydraulics, Madrid, Spain, IAHR, pp 881–888
dimensional depth-averaged models for bend-flow simulation. J
70. Playán E, Walker WR, Merkley GP (1994) Two-dimensional
Hydraul Eng 129(8):597–612
simulation of basin irrigation I: theory. J Irrig Drain Eng
46. Hubbard ME, Dodd N (2002) A 2D numerical model of wave 120(5):837–856
runup and overtopping. Coast Eng 47:1–26 71. Pope (1975) A more general effective-viscosity hypothesis. J
47. Hubbard ME, García-Navarro P (2000) Flux difference splitting Fluid Mech 472:331–340
and the balancing of source terms and flux gradients. J Comput 72. Puertas J, Pena L, Teijeiro T (2004) An experimental ap-
Phys 165:89–125 proach to the hydraulics of vertical slot fishways. J Hydraul Eng
48. Idelsohn SR, Storti MA, Oñate E (2001) Lagrangian formula- 130(1):10–23
tions to solve free surface incompressible inviscid fluid flows. 73. Rajaratnam N, Katopodis C, Solanski S (1992) New designs for
Comput Methods Appl Mech Eng 191:583–593 vertical slot fishways. Can J Civ Eng 19(3):402–414
49. Ingram RG, Chu VH (1987) Flow around islands in Rupert Bay: 74. Rajaratnam N, van der Vinne G, Katopodis C (1986) Hydraulics
an investigation of the bottom friction effect. J Geophys Res of vertical slot fishways. J Hydraul Eng 112(10):909–927
92(C13):14521–14533 75. Rameshwaran P, Shiono K (2003) Computer modelling of two-
50. Jia Y, Wang SSY (1999) Numerical model for channel flow and stage meandering channel flows. Water Marit Eng 156(4):325–
morphological change studies. J Hydraul Eng 125(9):924–933 339
51. Jirka GH (2001) Large scale flow structures and mixing 76. Rastogi AK, Rodi W (1978) Predictions of heat and mass transfer
processes in shallow flows. J Hydraul Res 39(6):567–573 in open channels. J Hydraul Div HY 3:397–420
52. Jones WP, Launder B (1972) The prediction of laminarization 77. Richardson JE, Panchang VG (1998) Three-dimensional sim-
with a two-equation model of turbulence. Int J Heat Mass Transf ulation of scour-inducing flow at bridge piers. J Hydraul Eng
15:301–314 124(5):530–540
53. LeVeque RJ (2002) Finite volume methods for hyperbolic prob- 78. Rider WJ, Kothe DB (1998) Reconstructing volume tracking. J
lems. Cambridge texts in applied mathematics, vol 31. Cam- Comput Phys 141:112–152
bridge University Press, Cambridge 79. Rodi W (1976) A new algebraic relation for calculating the
54. Lien HC, Hsieh TY, Yang JC, Yeh KC (1999) Bend-flow Reynolds stresses. Z Angew Math Mech 56:219–221
simulation using 2D depth-averaged model. J Hydraul Eng 80. Roe PL (1986) Discrete models for the numerical analysis of
125(10):1097–1108 time-dependent multidimensional gas dynamics. J Comput Phys
55. Lloyd PM, Stansby PK (1997) Shallow water flow around model 63:458–476
conical islands of small side slope. Part 1: Surface piercing. J 81. Sauvaget P, David E, Soares CG (2000) Modelling tidal currents
Hydraul Eng 123:1057–1067 on the coast of Portugal. Coast Eng 40:393–409
Depth Averaged Modelling of Turbulent Shallow Water Flow with Wet-Dry Fronts 341

82. Scardovelli R, Zaleski S (1999) Direct numerical simulation of 94. van Prooijen BC, Booij R, Uijttewaal WSJ (2000) Measurement
free surface and interfacial flows. Annu Rev Fluid Mech 31:567– and analysis methods of large scale horizontal coherent structures
603 in a wide shallow channel. In: 10th international symposium on
83. Sleigh PA, Gaskell PH, Berzins M, Wright NG (1998) An un- applications of laser techniques to fluid mechanics. Calouste Gul-
structured finite-volume algorithm for predicting flow in rivers benkian Foundation, Lisbon
and estuaries. Comput Fluids 27(4):479–508 95. Vázquez-Cendón ME (1996) An efficient upwind scheme with
84. Stansby PK (1997) Semi-implicit finite volume shallow-water finite volumes of the edge-type for the bidimensional shallow
flow and solute transport solver with k–ε turbulence model. Int J water equations. In: Benkhanldoun F, Vilsmeier R (eds) Finite
Numer Methods Fluids 25:285–313 volumes for complex applications. Problems and perspectives.
85. Stansby PK, Lloyd PM (2001) Wake formation around is- Hermes, Paris, pp 605–612
lands in oscillatory laminar shallow-water flows. Part 2. Three- 96. Vázquez-Cendón ME (1999) Improved treatment of source terms
dimensional boundary-layer modelling. J Fluid Mech 429:239– in upwind schemes for the shallow water equations in channels
254 with irregular geometry. J Comput Phys 148:497–526
86. Thomas FO, Goldschmidt VW (1986) Structural characteristics 97. Wilson CAME, Bates PD, Hervouet JM (2002) Comparison of
of developing turbulent planar jet. J Fluid Mech 63:227–256 turbulence models for stage-discharge rating curve prediction in
87. Toro EF (1999) Riemann solvers and numerical methods for fluid
reach-scale compound channel flows using two-dimensional fi-
dynamics, 2nd edn. Springer, New York
nite element methods. J Hydrol 257:42–58
88. Toro EF (2001) Shock-capturing methods for free-surface shal-
98. Winterwerp JC, Wang ZB, van Kester ATM, Verweij JF (2002)
low flows. Wiley, Chichester
89. Tryggvason G, Bunner B, Esmaeeli A, Juric D, Al-Rawahi N, Far-field impact of water injection dredging in the Crouch River.
Tauber W, Nas JHS, Yan YJ (2001) A front tracking method for Water Marit Eng 154(4):285–296
the computations of multiphase flow. J Comput Phys 169:708– 99. Wolanski E, Imberger J, Heron ML (1984) Island wakes in shal-
759 low coastal waters. J Geophys Res 89:10553–10569
90. Uijttewaal WSJ, Booij R (2000) Effects of shallowness on the de- 100. Wu W (2004) Depth-averaged two-dimensional numerical mod-
velopment of free-surface mixing layers. Phys Fluids 12(2):392– eling of unsteady flow and nonuniform sediment transport in
402 open channels. J Hydraul Eng 130(10):1013–1024
91. Uijttewaal WSJ, Girka GH (2003) Grid turbulence in shallow 101. Wu W, Rodi W, Wenka T (2000) 3D numerical modeling of
flows. J Fluid Mech 489:325–344 flow and sediment transport in open channels. J Hydraul Eng
92. Uijttewaal WSJ, Tukker J (1998) Development of quasi two- 126(1):4–15
dimensional structures in a shallow free-surface mixing layer. 102. Ye J, McCorquodale JA (1997) Depth-averaged hydrodynamic
Exp Fluids 24:192–200 model in curvilinear collocated grid. J Hydraul Eng 123(5):380–
93. Uittenbogaard R, van Vossen B (2001) 2D-DNS of quasi-2D tur- 388
bulence in shallow water. In: 3rd AFOSR international confer- 103. Yoon TH (2004) Finite volume model for two-dimensional
ence on direct numerical simulation and large eddy simulation shallow water flows on unstructured grids. J Hydraul Eng
(TAICDL), University of Texas at Arlington 130(7):678–688