Computational fluid dynamics and the physical modelling of an upland urban river

Geomorphology 44 (2002) 375 – 391 www.elsevier.com/locate/geomorph

Computational f luid dynamics and the physical modelling of an upland urban river Lin Ma a, Philip J. Ashworth b,*, James L. Best c, Lionel Elliott d, Derek B. Ingham d, Leslie J. Whitcombe c a

School of the Environment, University of Leeds, Leeds, LS2 9JT, West Yorkshire, UK b School of Geography, University of Leeds, Leeds, LS2 9JT, West Yorkshire, UK c School of Earth Sciences, University of Leeds, Leeds, LS2 9JT, West Yorkshire, UK d Department of Applied Mathematics, University of Leeds, Leeds, LS2 9JT, West Yorkshire, UK Received 18 April 2000; received in revised form 15 April 2001; accepted 16 April 2001

Abstract This paper describes the application of a commercially available, three-dimensional computational fluid dynamic (CFD) model to simulate the flow structure in an upland river that is prone to flooding. Simulations use a rectangular channel geometry, smooth sidewalls and a bed topography obtained from the field site that contains a subdued pool – riffle sequence. The CFD model uses the RNG j – e turbulence closure scheme of Yakhot and Orszag (J. Sci. Comput. 1 (1986) 1), as implemented in FLUENT 4.4.4, with a free surface. Results are shown for numerical runs simulating a 1:100 year return interval flood. Output from the numerical model is compared to a physical model experiment that uses a 1:35 scale fibreglass mould of the field study reach and measures velocity using ultrasonic Doppler velocity profiling (UDVP). Results are presented from the numerical and flume models for the water surface and streamwise velocity pattern and for the secondary flows simulated in the numerical model. A good agreement is achieved between the CFD model output and the physical model results for the downstream velocities. Results suggest that the streamwise velocity is the main influence on the flow structure at the discharge and channel configuration studied. Secondary flows are, in general, very weak being below the resolution of measurement in the physical model and less than 10% of the streamwise velocity in the numerical model. Consequently, there is no evidence for a ‘velocity dip’. It is suggested that the subdued topography or inlet morphology may inhibit the development of secondary flows that have been recorded in previous flat-bed, rectangular open channel flows. A significant corollary of these results is that the morphological evolution of the pool – riffle sequence at high discharges may be controlled primarily by the downstream distribution of velocity and sediment transport with little role for lateral sorting and sediment routing by secondary flows. This paper also raises a number of issues that may be of use in future CFD modelling of three-dimensional flow in open channels within the geomorphological community. D 2002 Elsevier Science B.V. All rights reserved. Keywords: Computational fluid dynamics; Pool – riffle; Secondary flows; Upland river; Physical model

*

Corresponding author. Division of Geography, School of the Environment, University of Brighton, Brighton, Sussex, BN2 4GJ, UK. E-mail addresses: [email protected] (L. Ma), [email protected] (P.J. Ashworth), [email protected] (J.L. Best), [email protected] (L. Elliott), [email protected] (D.B. Ingham), [email protected] (L.J. Whitcombe). 0169-555X/02/$ - see front matter D 2002 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 9 - 5 5 5 X ( 0 1 ) 0 0 1 8 4 - 2

376

L. Ma et al. / Geomorphology 44 (2002) 375–391

1. Introduction Upland urban channels in the UK are typically characterised by steep slopes, coarse bed material, straight or artificially lined banks and often experience high magnitude, short-lived floods. With increasing pressure from construction on land adjoining the floodplains of urban rivers, and the tendency for more of the catchment to be covered by an impervious surface, the risk from flooding in upland urbanised areas is becoming particularly acute. One zone of the catchment that is particularly vulnerable to flooding is

the point at which the upland river exits from the foothills into the low-lying floodplain. Traditionally, this location has been the optimal site for settlements because it represents a point for river crossing and a source of water and power for industry. Todmorden in northern England (Fig. 1a – b) is typical of many upland areas in the UK that have become established around a major river (in this case, the River Calder) but have suffered frequent and extreme flooding over the past two centuries. This paper presents the first results from a combined numerical and physical-scale modelling project that is investigating the flow and sediment

Fig. 1. Locations of (a) the study site in the UK, (b) the Calder catchment and town of Todmorden, and (c) study reach, 200 m downstream of the confluence of the Walsden Water and Calder rivers.

L. Ma et al. / Geomorphology 44 (2002) 375–391

dynamics of the River Calder as it passes through a straight reach with non-erodible banks and minimal spatial variation in bed topography. The study of two- and three-dimensional flow in open channels has recently experienced a surge of interest in the application of computational fluid dynamics (CFD) to hydrological, and in particular, geomorphological problems (e.g. Bates et al., 1997; Bradbrook et al., 1998; Hodskinson and Ferguson, 1998; Lane, 1998; Thompson et al., 1998; Lane et al., 1999; Nicholas and Sambrook Smith, 1999; Lane and Bates, 2000; Nicholas and McLelland, 2000; Wu et al., 2000). These modelling studies have shown that there is great potential in using CFD in complex channel geometries (e.g. Hodskinson and Ferguson, 1998; Wu et al., 2000), highly three-dimensional and non-uniform flow (e.g. Bradbrook et al., 1998, 2000) and at the interface between within-channel and overbank flows (e.g. Thomas and Williams, 1995a,b; Sofialidis and Prinos, 1999). However, the application of CFD models to natural river flows is not without problems, because (i) it involves an averaging of the Reynolds stresses and selection of an appropriate turbulence closure model (Lane, 1998), (ii) it is sometimes difficult to illustrate full mesh or grid independence (e.g. Hardy et al., 1999; Habashi et al., 2000), (iii) it is difficult to parameterise bed roughness (e.g. Nicholas, 2001), and (iv) there are few reliable field or laboratory databases for channels with non-uniform beds with which to compare and validate CFD output. Nevertheless, CFD modelling has tremendous potential for application to geomorphological and hydrological problems where the spatio-temporal boundary conditions may vary over a wide range, thus often making field study either impracticable or impossible. This paper illustrates the application of CFD to a semi-natural upland river channel that is also the focus of a series of physical scale-modelling experiments. Although flume models also involve assumptions and simplifications concerning flow and the treatment of wall roughness (cf. Peakall et al., 1996), such experiments do provide measurements of the boundary conditions for the CFD model and also quantify water surface elevation and flow structure for comparison to the numerical output. Previous studies have shown that a combined numerical – physical modelling approach is a useful method to validate numerical output,

377

although it can be difficult to isolate whether any differences between the two are caused by numerical simplifications or experimental errors (e.g. Moulin and Ben Slama, 1998; Olsen and Kjellesvig, 1999; Rajendran et al., 1999). The main objectives of this paper are therefore to (i) describe the new components and application of a CFD model for use in open channel flows, (ii) discuss some numerical issues that require further consideration by the CFD modelling community, (iii) compare CFD output to flume data derived for the same channel configuration and flow boundary conditions, and (iv) use the CFD model to describe the flow structure in a rectangular channel with bed undulations at high stage, a channel configuration typical of many upland UK urban channels.

2. Field site Fig. 1a –c shows the location and geometry of the field site on the River Calder, in the town of Todmorden, at the base of the Pennine Hills. The study reach is 175 m long, 10 m wide with 3.5 m high, near-vertical, stone sidewalls (Fig. 2). The high artificial banks easily contain the mean annual flood and prevent regular flooding of the nearby properties (see Fig. 2). Although the fixed sidewalls prevent the flow spreading laterally, the channel is free to rework the coarse alluvial bed material (mean surface D50 = 0.067 m) and creates a subtle, but obvious pool –riffle sequence (Fig. 3), with average bed slope of 0.0035. Vegetation may grow on the sidewalls (see Fig. 2), but this is regularly trimmed by the Environment Agency, and for the purposes of this study, the walls are assumed to be both vertical and smooth. Todmorden can experience high magnitude, short-lived floods (Fig. 4) that may cause extensive flooding of the nearby domestic and commercial properties and town high street. Peak discharges for the 2.33-, 10-, 50- and 100-year return interval flood at the Walsden Water/Calder confluence (Fig. 1c) are 40.7, 54.0, 80.0 and 87.6 m3 s
1, respectively. The 1:100 year return interval flood should just be contained within the channel crosssection with the 3.5-m-high stone walls, but because the top 1 m of the walls are permeable, extensive flooding may occur during floods of less than 1:50 return interval. An EDM survey of the channel bed

378

L. Ma et al. / Geomorphology 44 (2002) 375–391

Fig. 2. View of the study reach. Channel width is approximately 10 m. Picture taken from the footbridge in the study reach (location in Fig. 1c) looking upstream to the Walsden Water/Calder river confluence seen in the background of the picture. Note the near-vertical walls and vegetation that border the channel and the adjacent properties that are vulnerable to flooding (see text). Picture taken during winter with mean flow depth of approximately 0.5 m.

was undertaken at the start of the project with 717 points ( F 0.005 m) defining the three-dimensional bed topography for the flume model (Fig. 3, see below).

3. Physical scale modelling The three-dimensional EDM survey data were used to produce a precision engineered 1:35 scale model of the Todmorden study reach. The model was constructed of smooth fibreglass so that the main form

roughness was preserved but not the grain roughness. The model was fixed onto the base of a 0.3 m wide, 10 m long, 0.5 m deep, smooth-sided flume. The 5-m long model of the study reach was inserted in the central section of the flume and joined to two 2.5-m-long fibreglass sections at the up and downstream ends. Each 2.5 m joining section had the same cross-section throughout its length that was based on the model entrance (for the upstream block) and exit (downstream block) topography, with the join between all sections being smooth.

Fig. 3. Flume bed topography produced from interpolation of 717 EDM survey points of the Todmorden study reach and scaled by 1:35. Note the limited range in bed topography but also the clear demarcation of shallow pool (labelled P) and riffle (labelled R) sites.

L. Ma et al. / Geomorphology 44 (2002) 375–391

379

Fig. 4. Rainfall and discharge record for a characteristic flood at Todmorden on 4 – 8 March 1998. Only the water stage is recorded at Todmorden. Discharge is predicted using a runoff and river flow model by the North – East Environment Agency. The 1:100 year return period flood at Todmorden is approximately 88 m3 s
1 (see text).

Flume bed slope was set at 0.0025, which is lower than the field to ensure uniform flow within the 5 m model section. Discharge was held constant at 12.4 l s
1 (scale discharge of 89.9 m3 s
1), which is approximately the scaled equivalent of the 1:100 year flood. The three-dimensional flow structure was measured using ultrasonic Doppler velocity profilers (UDVP; see Takeda, 1991, 1993; Best et al., 2001) recording at either 4 MHz (downstream and vertical components) or 2 MHz (cross-stream component). This paper only reports on the downstream (u) velocity component for the flume data (see later). The downstream velocity was measured by fixing five UDVP’s in a horizontal plane facing upstream. The frame holding the five UDVP’s was moved laterally and vertically to measure 10 points across the flume width at 0.03-m intervals and 10 points in the vertical at 0.0075-m intervals. The UDVP rig enabled velocity measurements to be taken to within 0.015 m from the sidewalls and 0.025 m from the bed. Each UDVP measures the one-dimensional velocity (the component of flow towards the UDVP transducer) for 128 bins of 0.74 mm length away from the probehead at a resolution of F 4.8 mm s
1. Velocities from the first seven bins, which started 30 mm away from the probehead to avoid any possible stagnation of flow

at the probehead, were averaged to give a mean bin length and width of approximately 5 and 7 mm, respectively. Although velocity measurements were collected at 0.5-m sections in the 5-m flume model, this paper only uses the data for 0 (inlet), 2.0 (poolhead), 3.0 (riffle) and 5.0 m (outlet) sections. The water surface was measured using a point gauge accurate to F 0.0001 m and the mean Froude number (Fr) and Flow Reynolds Number (Re) for the flume run were 0.65 and 25,865, respectively.

4. Numerical modelling The three-dimensional CFD modelling of the fluid flow over the smooth but undulating bed topography of the Todmorden reach was performed using the commercial software package FLUENT 4.4.4. While a steady water discharge of 12.4 l s
1 has been used for the results presented in this paper, the solution technique that has been implemented allows simulation of the fluid flow when the discharge changes through a flood (e.g. Fig. 4) with only minor modifications to the code. In all cases, the simulated fluid flow with a steady water discharge converges to a steady water level from any reasonable given initial position.

380

L. Ma et al. / Geomorphology 44 (2002) 375–391

4.1. General considerations Flow in alluvial channels is inherently three-dimensional, with a complex relationship between the water surface, turbulence generation and bed topography. Most numerical simulations of river flows use a twodimensional depth-averaged model which diminishes the secondary flow occurring near the bed due to the 2D simulation of fluid flow over a three-dimensional bed topography (see review in Hervouet and Van Haren, 1996). In order to take into consideration the full three-dimensional aspect of the flow and to simulate more accurately the effect of bed topography, the approach of solving the full three-dimensional form of the Navier – Stokes equations is considered. Except for a few, relatively simple, fluid flow situations, such as plane mixing layers (e.g. Rogers and Moser, 1992), direct solutions of the Navier – Stokes equations using the Direct Simulation Method are very difficult given the present level of computer techniques and power, because the number of cells required to obtain meaningful and accurate numerical results is too large. With the advancement of the technique known as Large Eddy Simulation, some investigators have now simulated open channel flows in the absence of bed forms (Shi et al., 1999; see review in Meneveau and Katz, 2000). However, in this paper, turbulence closure models together with the Control Volume Method (see Rhie and Chow, 1983) have been employed in order to solve the Reynolds-averaged form of the Navier – Stokes equations. When solving for river flows, the water surfaces themselves will usually be part of the numerical solutions. Although the water surface and slope may be obtained by field survey or laboratory measurements for a steady river flow, and thus a ‘fixed lid’ may be superimposed on the free surface of the river in the numerical simulation, this is not always easy to undertake. When the river flow is in an unsteady state, such as during a flood and the water level is constantly changing, the water surface has to be predicted numerically. Even for a steady fluid flow, when the water surface elevation is not precisely available, an approximate ‘fixed lid’ may be specified on the water surface. Unfortunately, any incorrect ‘fixed lid’ will affect the distribution of the mass and momentum of the fluid flow in the numerical simulation, such that the fluid velocity may be over/

underestimated, hence leading to errors in the prediction of bed shear stress, and consequently, bedload transport. Therefore, in this paper, the Volume of Fluid Technique (VOF, Hirt and Nichols, 1981), has been incorporated into the solution of the Navier – Stokes equations in order to predict numerically the water surface elevation. In the commercial package FLUENT (see FLUENT User’s Guide Volume 4, 1996), the water flow and some of the air flow above the water are solved simultaneously. Governing fluid flow equations for air and water are expressed in a single form but with different physical properties. For a particular cell, the physical properties that are included in the equations are determined according to the mass fractions of the air and water in that cell. In addition, the complexity of implementing boundary conditions on the water surface has been avoided since the water surface is the interface between the air and the water. For an unsteady fluid flow, those cells currently filled with air provide the space for the water when the water level rises. 4.2. Governing fluid flow equations The governing Navier – Stokes equations for the unsteady, incompressible and turbulent fluid flow in a Cartesian co-ordinate system may be written as 4.2.1. Navier –Stokes equation 
i  i @ui @p @ @u @u j j @u q þ qu ¼
iþ j l þ @x @x @t @x j @x j @xi @sij ð1Þ þ j þ qgi @x 4.2.2. Continuity equation @ui ¼0 @xi

ð2Þ

where ui, p, q and l are the mean fluid velocity components, pressure, density and molecular viscosity of the fluid flow, respectively. sij is known as the Reynolds stress tensor, which represents the effects of the turbulent fluctuations in the fluid flow. As far as predicting the free surface by implementing the VOF, the volume fraction of the water, F, has

L. Ma et al. / Geomorphology 44 (2002) 375–391

been introduced in a computational cell and is defined as follows: F¼

dXwater dXcell

ð3Þ

where dVcell is the volume of the computational cell and dVwater is the fraction of the volume of the cell filled with water. Therefore, we have 8 F ¼ 1; cell is full of water > > > > < ð4Þ F ¼ 0; cell is full of air : > > > > : 0 < F < 1; cell contains a free surface The governing fluid flow equations for the water flow and the air flow above the water are expressed in a single form as given by Eqs. (1) and (2), but the physical properties that appear in Eqs. (1) and (2) are different and they are defined by the volume-fractionweighted average of physical properties of the air and water as follows: q ¼ Fqwater þ ð1
FÞqair

ð5Þ

l ¼ Flwater þ ð1
FÞlair :

ð6Þ

4.2.3. The Volume Fraction equation According to the law of mass conservation of air and water, the volume fraction of the water satisfies @F @F þ ui i ¼ 0 @t @x

ð7Þ

By numerically solving the volume fraction Eq. (7), the volume fractions of water, F, and of air, 1
F, in a control volume cell may be obtained. For more details of the implementation of the VOF, readers are referred to the FLUENT User’s Guide Volume 4 (1996). 4.3. Turbulence closure For closure of the governing equations, a turbulence model has to be introduced in order to model the Reynolds stress tensor sij. A number of different turbulence models are available within FLUENT, each involving different levels of complexity. The most commonly used is the standard j –e turbulence model,

381

which is based on the assumption of isotropic turbulence that makes the standard j – e model produce a large turbulent viscosity. However, the relatively newly developed RNG j – e turbulence model (see Yakhot and Orszag, 1986) has some improvements over the standard j – e model in simulating threedimensional fluid flows involving a large degree of rate of strain and boundary curvature. Since this situation is similar to that expected in the threedimensional, undulating bed of the River Calder at Todmorden, the RNG j –e model is used in this paper.

5. Boundary conditions and solution strategy for the CFD model Since the governing fluid flow equation is elliptic, boundary conditions around the solution domain are required in order to determine the solution. The upstream boundary conditions may be given according to the velocity distribution at the inlet of the computational domain to justify the rate of the discharge of the fluid flow. If the downstream fluid flow is fully developed, i.e. the water level does not change after the downstream boundary, then fully developed fluid flow conditions may be specified at the end of the computational domain. In the situation where the fluid flow cannot be regarded as being fully developed at the downstream boundary location, such as when simulating the fluid flow in a flume with a limited channel length, then a fluid velocity distribution may have to be specified at the downstream boundary. Correct boundary conditions at the downstream end are important as these may have a strong effect on the water level in the whole computational domain. For the results that are presented in this paper, the fluid velocity and water levels have been specified at the upstream and downstream boundaries using the data measured in the physical model. The turbulent kinetic energy and dissipation rate of the fluid flow at the upstream and downstream boundaries, which are usually not known, may be determined by specifying the turbulent intensity of the fluid flow. Nallasamy (1987) recommends using an empirical value of 10
20% of the downstream fluid velocity for the turbulent intensity to represent a wide variety of fluid flows that have not experienced intense interruptions to the three-dimensional flow

382

L. Ma et al. / Geomorphology 44 (2002) 375–391

or do not require any special treatment of the turbulence. Various values between 5% and 20% for the turbulence intensity at the upstream and downstream boundary locations have been investigated in this study. It is noted that the differences in the fluid velocities predicted using 5%, 10% and 20% turbulence intensity are graphically indistinguishable. However, employing different values for the turbulent intensity at the upstream and downstream boundaries do affect the predicted magnitude of the turbulent kinetic energy and dissipation rate of the fluid flow, but this appears to be limited to the region near the upstream and downstream boundaries. In this paper, a value of 10% is assumed for the turbulent intensity of the fluid flow at the upstream end of the domain and a value of 15% for the downstream end, to take into account the effect of the non-uniform bed topography on the fluid flow. Fluid flow near the bed and banks of a river is often very complex in terms of both its mean and turbulent structure. Although research is still required to improve understanding of these regions, the general principles of fluid mechanics show that fluid flow near these boundaries experience a rapid decrease in velocity in order to comply with the no-slip condition. Therefore, an extremely fine grid near these boundaries, together with excessive computational effort, is usually required in order to simulate their effect on the fluid velocity. In the case where the major interest has been on the characteristics of the main fluid flow rather than on the flow near the solid boundaries, wall functions have usually been employed in order to simulate these effects, and hence, the use of fine grids has been unnecessary. In this study, the standard wall function has been employed, which may be expressed as follows: uþ ¼

1 lnðEyþ Þ
DBðKsþ Þ j

uþ ¼ u=u*; yþ ¼ y=y*; Ksþ ¼ Ks =y*; pffiffiffiffiffiffiffiffiffiffi y* ¼ m=u* and u* ¼ sw =q

ð8Þ

where u is the fluid velocity parallel to the solid wall, y is the distance from the wall, u* is the wall friction velocity, Ks is the wall roughness height, E is a constant taking a value of 9.8, DB is an expression dealing with the hydraulically rough bed (see below),

j is Von Karman’s constant and usually takes a value of 0.4187, and m and sw are the dynamic viscosity of the fluid and the shear stress at the solid wall boundaries, respectively. Most river flows are in the fully rough regime, i.e. Ks+ > 60 and Eq. (8) may be employed. However, there is considerable uncertainty in the selection of an appropriate value for Ks in numerical simulations of river flows with strongly non-uniform bed material (cf. Nicholas, 2001). For example Hey (1979) and Bray (1980) found Ks = 3.5 D84, while Ferguson et al. (1989) reported Ks ranging from a minimum of 1.45 D50 to a median value of Ks = 4 D50. Clearly, further research is required on this issue, especially since the parameterisation of bed roughness and its representation in different mesh sizes is fundamental to the solution of near-boundary flow, shear stress and hence prediction of particle entrainment. In the study reported in this paper, a smooth bed and vertical banks have been employed in the flume model so that the contribution of grain roughness to Ks is zero and form roughness is negligible. This leads to DB = 0. In addition to the fluid velocities, the water level at the inlet boundary needs to be given by specifying the water volume fraction at the boundary, and this should be consistent with the water discharge rate through the channel. The water level at the downstream end may be predicted numerically according to the downstream boundary conditions (e.g. velocity boundary conditions) and in order to satisfy mass conservation. The initial conditions for the fluid flow field may be specified in an arbitrary way. However, for a steady discharge the experimental water level has been used in order to save on computational time. If the initial water surface is not given, then it is assumed that no water is in the computational domain at the beginning of the computation and the steady water level may be produced as time progresses. When the VOF technique is activated, the air – water two phase flow equations are solved. Therefore, boundary conditions for both the air and the water flows are required, although the air flow above the water is not of much interest. However, it should be noted that because the viscosity of air is approximately 50 times smaller than that of water, the air flow above the water has little effect on the water flow, but it permits space for the water level to rise and fall. On the other hand, allowing some airflow through the

L. Ma et al. / Geomorphology 44 (2002) 375–391

domain may help in maintaining computational stability, although the airflow conditions may be arbitrarily specified. In the numerical simulations performed in this study, a uniform air flow velocity, of magnitude equivalent to the water speed, has been specified at the inlet boundary of the channel and this allows air flow through the channel above the water. Boundary conditions for fully developed fluid flows have been imposed for the air flow at the downstream end. It is noted that if no air is allowed to enter through the upstream boundary, or if the air flow rate is too small, then large recirculating regions of air may occur above the water. This reversal of air flow back up the channel then causes computational instabilities. A Body Fitted Co-ordinate (BFC) grid system has been employed so that the undulating bed topography may be more accurately represented than when using a Cartesian grid. Fig. 5 illustrates the computational domain along the centre of the river channel that has been employed in the simulation. In all the results presented in this paper, the three-dimensional computational domain consists of 202  22  16 (x, y, z) grid cells, where the first, second and third integers denote the number of cells used in the downstream, lateral and vertical directions of the computational domain, respectively. Relatively fine grids have been distributed in the regions where large fluid velocity gradients are to be expected, such as in the regions near the bed and walls. Having set-up the grid and specified all the necessary boundary conditions, then the control volume method with the SIMPLE (Patankar, 1980) pressure – velocity iteration procedure, together with the commonly used Power Law scheme (Patankar, 1980; Olsen and Kjellesvig, 1999) for the discretisation of the convection terms, has been employed in order to solve the governing fluid flow equations numerically.

383

The water fraction equation is explicitly solved and the water surface updated at each time step before/ after the Navier – Stokes equations are solved. Although the second-order QUICK (Leonard, 1979) interpolation scheme can be activated, it is sometimes more difficult to obtain a stable water surface in the case studied than when using the Power Law scheme. A value of 0.001 has been employed as the convergence criteria at each time step for the sum of each of the normalised residuals over the whole fluid domain for all the governing fluid flow equations, including the mass residuals. For a steady discharge rate, such as that investigated in this paper, the steady state solution is assumed to have been reached when the changes in the water level are graphically indistinguishable.

6. Results and discussion of the CFD model A number of runs have been completed for both steady and unsteady fluid flows in the numerical and physical model of the river channel. The results presented in this paper are for a steady discharge rate of 12.4 l s
1 in the flume, which approximates to a 1:100 year discharge in the River Calder. A bed slope of 0.0025 was used in the flume and numerical model. The right-hand system co-ordinate system has been chosen such that the x-axis follows the average slope of the river, y-axis is across the river and z-axis is vertical. The u-component of the fluid velocities at the up and downstream boundaries has been specified using the velocity data collected by UDVP in the flume model at distances of 0 and 5 m downstream. The fluid velocities and water surfaces in the study reach are then obtained numerically and compared with the flume data at sections taken at 2 and 3 m.

Fig. 5. Computational grid employed by CFD model (202  22  16 including boundary cells). Note that model solves the position of the interface between water and air (see text). The grid shown is for the 5-m-long physical and numerical model.

384

L. Ma et al. / Geomorphology 44 (2002) 375–391

Three alternative grids with different resolutions have been tested and the results compared to the output for the 202  22  16 grid reported in this paper. Using a total of four different mesh sizes is consistent with the recommendations of Hardy et al. (1999), who suggested that any geomorphological or environmental modelling projects should construct at least four meshes of different spatial resolution. Grids of 101  11  8, 101  22  16, 202  22  16 and 202  22  32 have been employed and Fig. 6a shows that the numerically predicted water surface levels for the two grid sizes 101  11  8 and 202  22  16 are graphically indistinguishable. In addition, for these grid sizes, the streamwise velocities at the mean depth differ by a maximum of 5%. In strict terms, the solution is therefore, non-grid-independent, although the magnitude of the numerical errors is low. Clearly, further model runs are necessary in order to demonstrate a fully grid-independent solution, although it is sometimes difficult to illustrate full grid independence (e.g. Hardy et al., 1999; Habashi et al., 2000). In addition to the limitation imposed by the available computer capacity, different grid resolutions may pick up the influence of different scales of grain and form resistance and hence make it difficult to select an appropriate Ks value in the model. The discussion below refers only to the numerical results for the 202  22  16 grid. Fig. 6a – b show a reasonably good agreement (differences of less than 5%) between the numerically predicted water surface and the experimental data. The water surface across the channel is, in general, quite flat. This is a consequence of the large relative

roughness being modelled (in this case, ‘form’ roughness only). Fig. 7a – d show comparisons between the numerically predicted contours of the streamwise velocity and those of the flume data at two different downstream locations. The two cross-sections are located at 2 and 3 m from the location of the upstream boundary conditions and are positions that correspond to the poolhead and the downstream riffle, respectively (see Fig. 3). A reasonable agreement has been obtained between the numerical and experimental data for the streamwise fluid velocity that shows a simple distribution. Both the numerical and experimental results show a core of high fluid velocity near the centre of the channel, together with low velocity areas near to the solid boundaries of the channel. The core of the high velocity is slightly to the left-hand side of the channel centre at both sections and corresponds to flow going into and out of the pool located at f 2.5 m (see Fig. 3). Although the skew of the isovels to the left-side of the channel is slight, this may indicate some topographic steering of the flow. This spatial pattern of velocity is also consistent with the moderate channel asymmetry at the upstream section of the scale model that produces an uneven fluid input. It should be recalled that the 5-m model of the study reach had an upstream adjoining 2-m block with a uniform cross-section along its length that was identical to the upstream cross-section of the study reach model. Hence, the velocity distribution at the entrance to the model study reach already has an inherent asymmetry to the distribution of velocity across the channel. It may therefore be possible that the upstream morphology has more of an influence on the primary

Fig. 6. Position of the water surface measured in the 1:35 flume model and produced by the numerical model. (a) View downstream along the centreline of the channel and using two different grid sizes in the numerical model (see text), and (b) view across the channel (downstream is into the paper) at a location halfway down the 5-m-long model with the 202  22  16 grid only (see text).

L. Ma et al. / Geomorphology 44 (2002) 375–391

385

Fig. 7. Streamwise (u) isovels for two locations in the 5-m-long flume and numerical model. Isovels are shown for 0.025 m s
1 intervals. The view shown is looking downstream at the sections. The downstream positions of 2 m (a – b) and 3 m (c – d) correspond to a poolhead and riffle, respectively (see Fig. 3).

isovels than the subdued pool – riffle morphology at this flow discharge and depth. Furthermore, the high Reynolds numbers of the fluid flow in the numerical model make the governing equations almost parabolic in nature, which therefore makes the downstream boundary condition relatively unimportant. It should be stressed that the numerical results presented are non-grid-independent due to the limited number of computational cells employed in the simulation. Hence, errors in the prediction of the water level, as well as the detailed velocity distributions in the cross-sections of the channel are inevitable. A comparison between the numerically predicted streamwise fluid velocities and those of the flume data is shown in Fig. 8. Each data point in Fig. 8 represents the location of a UDVP measurement point at sections 2.0 and 3.0 m downstream in the physical model. The u-velocity at the same point in the numerical model is produced by interpolation of the CFD output onto the

10  10 UDVP grid. A very close clustering of points around the perfect agreement line indicates a reasonably good correlation (r 2 = 0.79) between the numerical and experimental data. One of the features of fully developed turbulent flow in open channels with no bed topography is that the maximum downstream velocity at the channel centre is not at the water surface, but rather just below it, i.e. at about 0.6 – 0.7 of the water depth. This phenomenon, known as the ‘velocity-dip’, has been found both experimentally and theoretically (Rajaratnam and Muralidhar, 1969; Sarma et al., 1983; Tominaga et al., 1989; Nezu and Nakagawa, 1993). The main cause of the velocity-dip is due to the momentum transport caused by secondary flow that occurs in the fully developed turbulent open channel flow (Tominaga et al., 1989; see review in Younis, 1996). Fig. 9a – b show the secondary flow vectors for sections 2 and 3 m for the numerical simulations only.

386

L. Ma et al. / Geomorphology 44 (2002) 375–391

Fig. 8. Correlation between the u velocity measured in the flume model and the CFD model. Each data point represents the location of a UDVP measurement point at sections 2.0 and 3.0 m downstream in the physical model and the interpolation of the CFD model output onto the 10  10 UDVP grid.

It was not possible to quantify the secondary flow vectors in the flume because the cross-stream (v) velocity component was so weak that it was below the resolution of the UDVP. Fig. 9a– b show several key features. Firstly, the magnitude of secondary flows is very small, reaching a maximum of 10% of the downstream flow, with many points being characterised by values of less than 5% of the primary velocity. The secondary flow is strongest near the bed and regions of the bed where there is the greatest change in topography. The weak secondary flow, except near the bed, is similar to the values of less than 2% and 4% of the maximum downstream velocity found by Nezu and Nakagawa (1984) and Nezu et al. (1993), respectively. Plots of the planform uv vectors (Fig. 10) near the bed and surface illustrate that these very weak secondary flows have very little effect on the resolved flow direction, with flow essentially parallel to the banks at both the surface and near the bed. A slight convergence of flow into the pool at 2.5 m is apparent, but any topographic effects on the secondary flow are minimal. A significant corollary of this result is that sediment transport would be expected to occur parallel to the banks, the

pattern of transport probably being dominated by any spatial variation in bed shear stress rather than routing by secondary flows. Fig. 9a –b show no evidence for a velocity dip or associated surface convergent cells near the surface at the centre of the channel. As discussed above, this finding is in contrast to previous work on straight channels with no bed topography (e.g. Nezu and Rodi, 1986). Furthermore, the downstream isovels and secondary flows vectors (Figs. 7 and 9) suggest that there are no ‘corner cells’ generated by this bed topography, this also being dissimilar to past work on straight open channels (e.g. Nezu and Nakagawa, 1984; Nezu et al., 1988, 1993). Although of low magnitude, the secondary flows in the numerical results do show a general circulation upwards from the riffle and down into the poolhead on the left-side of the channel (Fig. 9a) and then circulation in an opposite direction from the downstream riffle and into the next poolhead (Fig. 9b) that is located on the right-hand side of the channel (see Fig. 3). Such a reversal in the circulation direction of the secondary flow has been observed by many workers for flow around meander bends and from pool-topool (e.g. Hey and Thorne, 1975; Dietrich and Smith, 1983). The precise location of the point at which the secondary flow changes circulation direction is usually associated with the riffle (cf. Thompson, 1986: fig. 2, p. 634) or the position where the bed begins to slope into the next pool. The secondary flow vectors for the study reach examined here show that the flow circulation does change direction from pool to riffle/ poolhead but that the magnitude of the secondary flows is very low and strongest near the bed. In anisotropic turbulence-driven secondary circulation, turbulence is damped near the free surface and a damping factor of 0.4 – 0.8, or a modified value for the constant cl is employed in the j –e model (Nezu and Nakagawa, 1993). In contrast, a surface-correction model for the pressure-strain term is usually used in the Reynolds-stress model (Gibson and Rodi, 1989). In this study, the inclusion of air above the water has made the free surface an interface between the air and water. This automatically introduces some kind of damping effect into the numerical scheme due to the differences in the physical properties of the two phases. Therefore, a vertical parabolic type of behaviour for the eddy viscosity (Nezu and Rodi, 1986)

L. Ma et al. / Geomorphology 44 (2002) 375–391

387

Fig. 9. Secondary flow vectors for sections (a) 2 m, and (b) 3 m from the numerical model. Primary flow direction is into the page. Note that the secondary flows are generally less than 5% of the primary velocities shown in Fig. 7a and c.

appears in the results of the numerical simulations reported here. The numerical simulations and physical model both show that secondary flows are extremely weak in the subdued pool – riffle topography of the study reach. One reason why this may be so for the physical model is because the flume may have an insufficient length of inlet before the study reach is encountered, thereby inhibiting the full development of turbulence-driven secondary flows. However, this is one of the reasons why the RNG j –e model was employed in the numerical scheme that also showed weak secondary flows. On the other hand, the study reach model is straight and has a relatively subtle bed topography that prevents strong ‘pressure-driven’ secondary flows from occurring. Hence, there is no clear evidence of a turbulencedriven ‘velocity dip’ at approximately 70% of the water depth (see Fig. 7a– d), because there are no significant

turbulence-driven secondary flow vortices that are supposed to be responsible for the ‘velocity dip’. The dominant control on the velocity distribution in the study reach is the primary u-velocity component. However, the shift of this core of high-speed fluid towards the left-hand side of the channel and below the water surface is suggested to be a consequence of the asymmetric profile of the input fluid velocity. The asymmetry of the bed topography in the study reach model and prior to the upstream end of the scale model appears to dominate the fluid flow structure throughout the flume. Undoubtedly, this is partly because the pool – riffle topography is so subdued (although it has a maximum amplitude of 0.35 m in the field, see Fig. 3), but also because the 1:100 recurrence interval flood drowns out the bed form roughness. Numerical and experimental runs are currently underway with an enhanced pool – riffle top-

388

L. Ma et al. / Geomorphology 44 (2002) 375–391

Fig. 10. Planform (uv) velocity vectors from the numerical model for (a) near-surface (second flow cell at 10 mm from water surface), and (b) near-bed (bottom flow cell, 2 mm from bed) in the 5 m study reach. The bed topography displayed in Fig. 3 is superimposed to show the locations of pools (P) and riffles (R). Note the difference in velocity vector scales for (a) and (b). Results show minimal flow convergence or divergence over the undulating bed topography.

ography from a different field prototype to investigate the impact of relative bed emergence and more pronounced bed morphology on the fluid flow structure. Finally, the results from this study of a laterally confined channel may have wider application to other rivers with stable or fixed banks (e.g. bedrock channels or engineered structures) or channels with low width/ depth ratios at high flow (e.g. incised rivers). Such rivers may also experience negligible secondary flows in rare, high floods (e.g. 1:100 recurrence interval) and the bed evolution and sediment transport may largely reflect the spatial and temporal pattern of bed shear stress imposed solely by the downstream flow.

7. Conclusions An experimental (flume) and numerical (computational fluid dynamics, CFD) model has been produced of an upland urban channel that experiences frequent

and severe flooding. For the purposes of the model, the channel geometry is assumed to be straight with smooth vertical walls. The channel bed has a subdued pool –riffle sequence and both the flume and CFD models use a fixed bed that has no bed roughness. Results are presented for model runs using a 1:100 year return period flood. The main conclusions are: (1) Realistic simulations of the streamwise velocity are produced both by the flume and numerical models. The magnitude and spatial distribution of streamwise velocity is in close agreement between flume and numerical model. (2) The use of four different mesh sizes in the numerical model illustrates that errors in the water surface and mean velocity less than 5%. In strict terms, the numerical solution is non-grid-independent, although the magnitude of the numerical errors is low. (3) Both numerical and physical model results show a core of high streamwise velocity that is offcentre. This asymmetry is probably caused by the topography of the channel cross-section at the up-

L. Ma et al. / Geomorphology 44 (2002) 375–391

stream input point. Highest streamwise velocities occur near the water surface. (4) Secondary flows in the numerical model are weak (less than 10% of the streamwise velocity) and confined to regions near the bed and where there is the greatest change in bed topography. There is no evidence for either a ‘velocity dip’ or significant ‘corner cells’ as has been documented in past work on rectangular open channels with no bed topography. (5) Results from both the physical and numerical models suggest that the flow structure at high discharges, in laterally confined channels with a subdued bed topography, is dominated by the primary velocity. Hence, bed morphological changes and sediment transport will be predominantly in the downstream direction. A number of issues are currently under investigation in an ongoing physical and numerical modelling programme of the Todmorden reach. Priorities are to investigate the effects of grain roughness and changing topographic emergence (relative relief) on the flow structure at different discharges and model bedload transport for uniform grain sizes. Numerical runs will be performed with both a free surface j – e and a Reynolds stress model. Nomenclature CFD D50 E Fr Re UDVP VOF x y z DB Ks u v u*

computational fluid dynamics median bed particle size (m) constant in Eq. (8) Froude number flow Reynolds Number ultrasonic Doppler velocity profilers volume of fluid technique for solving the water surface downstream direction cross-stream direction vertical direction an expression for a rough bed in the wall function wall roughness height (m) mean streamwise velocity component (m s
1) mean cross-stream velocity component (m s
1) wall friction velocity (m s
1)

y j– e j l m sw cl F p q sij Vcell Vwater t

389

distance from the wall in Eq. (8) (m) turbulence model Von Karman’s constant molecular viscosity of the fluid flow (kg m
1 s
1) dynamic viscosity of the fluid (m2 s
1) wall shear stress (N m
2) empirically obtained constant for the eddy viscosity volume fraction of water in a computational cell pressure of the fluid flow (N m
2) density of the fluid flow (kg m
3) Reynolds stress tensor (N m
2) volume of the computational cell volume of the computational cell filled with water time

Acknowledgements This research was funded by EPSRC grant GR/ L90590. The authors are grateful to David Ashley and Mark Franklin for their help with the field and flume work, respectively. David Pelleymounter (Environment Agency, North – East Region) and Andrew Pepper (EPSRC Coastal, Estuarine and Waterway Programme Coordinator) provided valuable input and feedback during fieldwork and development of the CFD model, but all of the ideas presented in this paper are solely the responsibility of the authors. We are extremely grateful to Paul Bates (Bristol) and Paul Carling (Southampton) for their perceptive and thorough reviews. References Bates, P.D., Horritt, M.S., Smith, C.N., Mason, D., 1997. Integrating remote sensing observations of flood hydrology and hydraulic modelling. Hydrol. Processes 11, 1777 – 1795. Best, J.L., Kirkbride, A., Peakall, J., 2001. Mean flow and turbulence structure in sediment-laden gravity currents: new insights using ultrasonic Doppler velocity profiling. In: McCaffrey, W.D., Peakall, J., Kneller, B.C. (Eds.), Particulate Gravity Currents, Special Publication of the International Association of Sedimentologists, 31, 159 – 172 pp. Bradbrook, K.F., Biron, P.M., Lane, S.N., Richards, K.S., Roy, A.G., 1998. Investigations of controls on secondary circulation

390

L. Ma et al. / Geomorphology 44 (2002) 375–391

and mixing processes in a simple confluence geometry using a three-dimensional numerical model. Hydrol. Processes 12, 1371 – 1396. Bradbrook, K.F., Lane, S.N., Richards, K.S., 2000. Numerical simulation of three-dimensional time-averaged flow structure at river channel confluences. Water Resour. Res. 36 (9), 2731 – 2746. Bray, D.I., 1980. Evaluation of effective boundary roughness for gravel-bed rivers. Can. J. Civ. Eng. 7, 392 – 397. Dietrich, W.E., Smith, J.D., 1983. Influence of the point bar on flow though curved channels. Water Resour. Res. 20 (19), 1173 – 1192. Ferguson, R.I., Prestegaard, K.L., Ashworth, P.J., 1989. Influence of sand on hydraulics and gravel transport in a braided gravel bed river. Water Resour. Res. 25, 635 – 643. FLUENT User’s Guide Volume 4, 1996. Release 4.4. Fluent Europe, Sheffield, UK, 258 pp. Gibson, M.M., Rodi, W., 1989. Simulation of free surface effects on turbulence with a Reynolds stress model. J. Hydraul. Res., IAHR 27, 233 – 244. Habashi, W.G., Dompierre, J., Bourgault, Y., Yahia-A-A, D., Fortin, M., Vallet, M.-G., 2000. Anisotropic mesh adaption: towards user-independent, mesh-independent and solver-independent CFD: Part I. General principles. J. Numer. Methods, Fluids 32, 725 – 744. Hardy, R.J., Bates, P.D., Anderson, M.G., 1999. The importance of spatial resolution in hydraulic models for floodplain environments. J. Hydrol. 216, 124 – 136. Hervouet, J.M., Van Haren, L., 1996. Recent advances in numerical methods for fluid flows. In: Anderson, M.G., Walling, D.E., Bates, P.D. (Eds.), Floodplain Processes. Wiley, Chichester, pp. 183 – 215. Hey, R.D., 1979. Flow resistance in gravel-bed rivers. J. Hydraul. Div., Am. Soc. Civ. Eng. 105, 365 – 379. Hey, R.D., Thorne, C.R., 1975. Secondary flows in river channels. Area 7, 191 – 195. Hirt, C.W., Nichols, B.D., 1981. Volume of Fluid (VOF) method for the dynamics of free boundaries. J. Comput. Phys., 201 – 225. Hodskinson, A., Ferguson, R.I., 1998. Numerical modelling of separated flows in river bends: model testing and experimental investigation of geometric controls on the extent of flow separation at the concave bank. Hydrol. Processes 12, 1323 – 1338. Lane, S.N., 1998. Hydraulic modelling in hydrology and geomorphology: a review of high resolution approaches. Hydrol. Processes 12, 1131 – 1150. Lane, S.N., Bates, P.D., 2000. Introduction. In: Bates, P.D., Lane, S.N. (Eds.), High Resolution Flow Modelling in Hydrology and Geomorphology. Wiley, Chichester, pp. 1 – 12. Lane, S.N., Bradbrook, K.F., Richards, K.S., Biron, P.A., Roy, A.G., 1999. The application of Computational Fluid Dynamics to natural river channels: three-dimensional versus two-dimensional approaches. Geomorphology 29, 1 – 21. Leonard, B.P., 1979. A stable and accurate convective modelling procedure based on quadratic upstream interpolation. Comput. Methods Appl. Mech. Eng. 19, 59 – 98. Meneveau, C., Katz, J., 2000. Scale-invariance and turbulence models for large-eddy simulation. Annu. Rev. Fluid Mech. 32, 1 – 32.

Moulin, C., Ben Slama, E., 1998. The two-dimensional transport module SUBIEF. Applications to sediment transport and water quality processes. Hydrol. Processes 12, 1183 – 1195. Nallasamy, M., 1987. Turbulence models and their applications to the prediction of internal flows—a review. Comput. Fluids 15, 151 – 194. Nezu, I., Nakagawa, H., 1984. Cellular secondary currents in channel flow. J. Hydraul. Eng., Am. Soc. Civ. Eng. 110, 173 – 193. Nezu, I., Nakagawa, H., 1993. Turbulence in Open Channel Flows. Balkema, Rotterdam, The Netherlands 279 pp. Nezu, I., Rodi, W., 1986. Open-channel flow measurements with a laser Doppler anemometer. J. Hydraul. Eng., Am. Soc. Civ. Eng. 112, 335 – 355. Nezu, I., Nakagawa, H., Kawashima, N., 1988. Cellular secondary currents and sand ribbons in fluvial open channel flows. Proc. 6th Congr., APD – IAHR, vol. II, Kyoto, Japan, pp. 51 – 58. Nezu, I., Tominaga, A., Nakagawa, H., 1993. Field measurements of secondary currents in straight rivers. J. Hydraul. Eng., Am. Soc. Civ. Eng. 119, 598 – 614. Nicholas, A.P., 2001. CFD modelling of boundary roughness in gravel bed rivers: an investigation of the effects of random variability in bed elevation. Earth Surf. Processes Landforms 26 (4), 345 – 362. Nicholas, A.P., McLelland, S.J., 2000. Hydrodynamics of a floodplain recirculation zone investigated by field monitoring and numerical simulation. In: Marriott, S.B., Alexander, J. (Eds.), Floodplains: Interdisciplinary Approaches. Geol. Soc. London, Spec. Publ., vol. 163, pp. 13 – 26. Nicholas, A.P., Sambrook Smith, G.H., 1999. Numerical simulation of three-dimensional flow hydraulics in a braided channel. Hydrol. Processes 13, 913 – 929. Olsen, N.R.B., Kjellesvig, H.M., 1999. Three-dimensional numerical modelling of bed changes in a sand trap. J. Hydraul. Res. 37, 189 – 198. Patankar, S.V., 1980. Numerical Heat Transfer and Fluid Flow. Hemisphere Publishing, Washington, DC, 197 pp. Peakall, J., Ashworth, P.J., Best, J.L., 1996. Physical modelling in fluvial geomorphology: principles, applications and unresolved issues. In: Rhoads, B.L., Thorn, C.E. (Eds.), Scientific Nature of Geomorphology. Wiley, Chichester, pp. 221 – 253. Rajaratnam, N., Muralidhar, D., 1969. Boundary shear stress distribution in rectangular open channels. Houille Blanche 6, 603 – 609. Rajendran, V.P., Constantinescu, S.G., Patel, V.C., 1999. Experimental validation of numerical model of flow in pump-intake bays. J. Hydraul. Eng., Am. Soc. Civ. Eng. 125, 1119 – 1125. Rhie, C.M., Chow, W.L., 1983. Numerical study of the turbulent flow past an airfoil with trailing edge separation. J. Am. Inst. Aeronaut. Astronaut. 21, 1525 – 1532. Rogers, M.M., Moser, R.D., 1992. The three-dimensional evolution of a plane mixing layer—the Kelvin – Helmholtz rollup. J. Fluid Mech. 243, 183 – 226. Sarma, K.V.N., Lakshminarayana, P., Rao, N.S.L., 1983. Velocity distribution in smooth rectangular open channels. J. Hydraul. Eng., Am. Soc. Civ. Eng. 109, 270 – 289. Shi, J., Thomas, T.G., Williams, J.J.R., 1999. Large-eddy simulation of flow in a rectangular open channel. J. Hydraul. Res. 37, 345 – 361.

L. Ma et al. / Geomorphology 44 (2002) 375–391 Sofialidis, D., Prinos, P., 1999. Numerical study of momentum exchange in compound open channel flow. J. Hydraul. Eng., Am. Soc. Civ. Eng. 125, 152 – 165. Takeda, Y., 1991. Development of an ultrasound velocity profile monitor. Nucl. Eng. Des. 126, 277 – 284. Takeda, Y., 1993. Velocity profile measurement by ultrasonic Doppler method. In: Kelleher, M.D. (Ed.), Experimental Heat Transfer, Fluid Mechanics and Thermodynamics. Elsevier, Amsterdam, pp. 126 – 131. Thomas, T.G., Williams, J.J.R., 1995a. Large eddy simulation of turbulent flow in an asymmetric compound open channel. J. Hydraul. Eng., Am. Soc. Civ. Eng. 33, 27 – 41. Thomas, T.G., Williams, J.J.R., 1995b. Large eddy simulation of a symmetric trapezoidal channel at Reynolds number of 430,000. J. Hydraul. Eng., Am. Soc. Civ. Eng. 33, 825 – 842. Thompson, A.L., 1986. Secondary flows and the pool – riffle unit: a

391

case study of the processes of meander development. Earth Surf. Processes Landforms 11, 631 – 641. Thompson, D.M., Nelson, J.M., Wohl, E.E., 1998. Interactions between pool geometry and hydraulics. Water Resour. Res. 12, 3673 – 3681. Tominaga, A., Nezu, I., Ezaki, K., Nakagawa, H., 1989. Threedimensional turbulent structure in straight open channel flows. J. Hydraul. Res. 27, 149 – 173. Wu, W., Rodi, W., Wenka, T., 2000. 3D numerical modelling of flow and sediment transport in open channels. J. Hydraul. Eng., Am. Soc. Civ. Eng. 126, 4 – 15. Yakhot, V., Orszag, S.A., 1986. Renormalization group analysis of turbulence: I. Basic theory. J. Sci. Comput. 1, 1 – 51. Younis, B.A., 1996. Progress in turbulence modelling for openchannel flows. In: Anderson, M.G., Walling, D.E., Bates, P.D. (Eds.), Floodplain Processes. Wiley, Chichester, pp. 299 – 332.