0 Lattice Boltzmann Simulation for Shallow Water Flow Applications

Under the influence of gravity, many free-surface flows can be modelled by the well-known shallow water equations under the assumption that the vertical scale is much smaller than any typical horizontal scale and the pressure is hydrostatic. These equations can be derived from the depth-averaged incompressible Navier-Stokes equations and usually include continuity and momentum equations. Hence, the applications of depth-averaged models include a wide spectrum of phenomena in hydraulic flows such as ocean circulation modelling Salmon (1999a) and wind-driven ocean circulation Zhong et al. (2005) to name but a few. Simulation of such real-world flow problems is not trivial since the geometry can be complex and the topography irregular. Numerical methods based on the finite difference, the finite volume or the finite element methods have been applied to simulate the shallow water equations, refer to Bermúdez & Vázquez (1994); Kurganov & Levy (2002); LeVeque (1998); Stansby & Zhou (1998); Toro (1992); Vázquez-Cendón (1999); Vukovic & Sopta (2002); Xing & Shu (2006); Zhou (1995) among others. For most of these approaches, the treatment of bed slopes and friction forces often cause numerical difficulties in obtaining accurate solutions, see, for example, Bermúdez & Vázquez (1994); LeVeque (1998); Vázquez-Cendón (1999). In addition the extension of these schemes to complex geometries is not trivial, refer to Benkhaldoun et al. (2007), for example. Some of these approaches are very expensive if one considers real flows Vukovic & Sopta (2002). Since the problems are posed at a large scale it has been the aim of practitioners to develop a simple and accurate representation of the source terms in order to simulate practical shallow water flows without relying on upwind discretization or Riemann problem solvers, refer to Benkhaldoun et al. (2009; 2010); Benkhaldoun & Seaïd (2010) and references therein for these alternatives. The idea of this chapter will be to give the reader a self-contained introduction of the developments and the implementation of the shallow water lattice Boltzmann approach. In this chapter the lattice Boltzmann method will be applied to the simulation of depth-averaged models in flow hydraulics and dispersion Banda et al. (2009); Klar et al. (2008); Seaïd & Thömmes (2009); Thömmes et al. (2007). It can be pointed out that the shallow water equations referred to in this discussion are viscous and also account for the effects of bed slope, bed friction, Coriolis forces and wind stresses in two-dimensional simulations Dellar (2002); Salmon (1999a); Zhou (2002). The practical aspects 11


Introduction
Under the influence of gravity, many free-surface flows can be modelled by the well-known shallow water equations under the assumption that the vertical scale is much smaller than any typical horizontal scale and the pressure is hydrostatic.These equations can be derived from the depth-averaged incompressible Navier-Stokes equations and usually include continuity and momentum equations.Hence, the applications of depth-averaged models include a wide spectrum of phenomena in hydraulic flows such as ocean circulation modelling Salmon (1999a) and wind-driven ocean circulation Zhong et al. (2005) to name but a few.Simulation of such real-world flow problems is not trivial since the geometry can be complex and the topography irregular.Numerical methods based on the finite difference, the finite volume or the finite element methods have been applied to simulate the shallow water equations, refer to Bermúdez & Vázquez (1994); Kurganov & Levy (2002); LeVeque (1998); Stansby & Zhou (1998); Toro (1992); Vázquez-Cendón (1999); Vukovic & Sopta (2002); Xing & Shu (2006); Zhou (1995) among others.For most of these approaches, the treatment of bed slopes and friction forces often cause numerical difficulties in obtaining accurate solutions, see, for example, Bermúdez & Vázquez (1994); LeVeque (1998); Vázquez-Cendón (1999).In addition the extension of these schemes to complex geometries is not trivial, refer to Benkhaldoun et al. (2007), for example.Some of these approaches are very expensive if one considers real flows Vukovic & Sopta (2002).Since the problems are posed at a large scale it has been the aim of practitioners to develop a simple and accurate representation of the source terms in order to simulate practical shallow water flows without relying on upwind discretization or Riemann problem solvers, refer to Benkhaldoun et al. (2009;2010); Benkhaldoun & Seaïd (2010) and references therein for these alternatives.The idea of this chapter will be to give the reader a self-contained introduction of the developments and the implementation of the shallow water lattice Boltzmann approach.In this chapter the lattice Boltzmann method will be applied to the simulation of depth-averaged models in flow hydraulics and dispersion Banda et al. (2009); Klar et al. (2008); Seaïd & Thömmes (2009); Thömmes et al. (2007).It can be pointed out that the shallow water equations referred to in this discussion are viscous and also account for the effects of bed slope, bed friction, Coriolis forces and wind stresses in two-dimensional simulations Dellar (2002); Salmon (1999a); Zhou (2002).The practical aspects of the method will be emphasized.In addition the methods for coupling the shallow water flow to other mass balance equations like pollutant transport and temperature dispersion will also be discussed.Considerations will also be made for new developments in the fields and further possible extensions.
The lattice Boltzmann (LB) method, also popularly referred to as LBM, is an alternative numerical tool for simulating fluid flows Chen & Doolen (1998).The method is based on statistical physics and models the fluid flow by tracking the evolution of the distribution functions of the fluid particles in discrete phase space.The essential approach in the LB method lies in the recovery of macroscopic fluid flows from the microscopic flow behaviour of the particle movement or the mesoscopic evolution of particle distributions.The basic idea is to replace the nonlinear differential equations of macroscopic fluid dynamics by a simplified description modeled on the kinetic theory of gases.Furthermore, the LBM offers several desirable properties such as linear convection terms and nearest-neighbor stencils.On a structured mesh, the LBM can be implemented in a two-stage procedure namely, a collision operator evaluation which involves only local operations, and an advection operation where values are transported to adjacent lattice points without performing any computations.In this chapter, the dynamics of the two different but dependent models namely, (i) a depth-averaged hydrodynamic model defining the flow, and (ii) a depth-averaged advection-diffusion model defining the transport of the pollutant/temperature are solved by an LBM with two distribution functions modelling the dynamics of the hydrodynamic flow and the pollutant concentration or temperature, respectively, will be discussed.To obtain the hydrodynamic behaviour, the Chapman-Enskog expansion which exploits a small mean free path approximation to describe slowly varying solutions of the underlying kinetic equations, is undertaken Zhou (2004).The method has been proven to be effective for simulating flows in complicated geometries and implementation on parallel computer architectures Kandhai et al. (1998).Furthermore, the method has become an alternative to other numerical methods like finite difference, finite element and finite volume methods in computational fluid dynamics.
As such the LB method has found a wide range of applications in a variety of fields, which include numerical simulation of shallow water equations.The LB method has been successfully adopted to simulate shallow water equations which describe wind-driven ocean circulation Salmon (1999a); Zhong et al. (2005), to model three-dimensional planetary geostrophic equations Salmon (1999b), and to study atmospheric circulation of the northern hemisphere with ideal boundary conditions Feng et al. (2002).In Klar et al. (2008); Thömmes et al. (2007), a practical study of the LBM for shallow water flows and pollutant dispersion by tidal flow problems in complex geometry and irregular bathymetry was presented.The bathymetry is given either by an analytical function or by data points in a two-dimensional domain.In the dispersion of pollutants Banda et al. (2009) the flow characteristics and concentration profiles of dispersive species will be discussed in detail.It can be noted that all above LB methods have been mainly applied to the isothermal shallow water flows and no thermal sources have been accounted for.In Seaïd & Thömmes (2009), a presentation of shallow water equations involving thermal sources has been made.
In the next section a brief discussion of the shallow water equations will be presented.This will be followed by Section 3 in which the practical details of implementing the LBM will be presented.Section 4 will present numerical results.The chapter concludes with Section 5.
Lattice Boltzmann Simulation for Shallow Water Flow Applications 3

Shallow water equations in hydraulics with coupling to dispersion
In this section a brief discussion of the derivation of the shallow water equations will be presented.This will be extended to the modelling of pollutant and temperature dispersion.In general, modelling of fluid flow with dispersion on a free-surface requires two sets of coupled partial differential equations.The first set of equations describes the water motion on the free-surface flow while, the second set of equations models the distribution of a pollutant or temperature on the water free-surface.
The flow is governed by the depth-averaged Navier-Stokes equations involving several assumptions including (i) the domain is shallow enough to ignore the vertical effects, (ii) the pressure is hydrostatic, (iii) all the water properties are assumed to be constant with the exception that in temperature dispersion the density is temperature dependent, which is accounted for using the Boussinesq approximation, and (iv) viscous dissipation of energy is ignored and any radiative heat losses are assumed to have occurred over a time scale small compared with that which characterizes the flow motion.

Equations of depth-averaged models in hydraulics
The starting point for the discussion of depth-averaged models in hydraulic flows is the three-dimensional incompressible Navier-Stokes equations, where t is the time variable, (x, y, z) T are the space coordinates, ρ the density, (u, v, w) T the velocity field, p the pressure, Ω the Coriolis parameter defined by Ω = 2ω sin φ, with ω denoting the angular velocity of the earth and φ is the geographic latitude, g is the gravitational force, ν H and ν V are the coefficients of horizontal and vertical viscosity, respectively.In (1), ∆ = ∂ 2 ∂x 2 + ∂ 2 ∂y 2 denotes the two-dimensional Laplace operator.In most shallow water models, the ratio of vertical length scale to horizontal length scale is very small.As a consequence, the horizontal viscosity terms are typically orders of magnitude smaller than the vertical viscosity terms and their effect is normally small and obscured by numerical diffusion.Therefore, most models either neglect these terms or simply use a constant horizontal viscosity coefficient.In addition, assuming that the pressure is hydrostatic, the momentum equation in the vertical direction (1d) degenerates into the following form Integrating vertically the continuity equation (1a) from the bottom topography z = Z to the free surface z = h + Z, and using the kinematic condition at the free surface leads to the free surface equation where h(x, y, t)+Z(x, y) is the water surface elevation and Z(x, y) is the bed or bottom as depicted in Figure 1.The boundary conditions at the water free surface are specified by the prescribed wind stresses T W x and T W y The wind stresses T W x and T W y are given by a quadratic function of the wind velocity (W x , W y ) T as 5) where C W is the coefficient of wind friction defined by Bermúdez & Vázquez (1994) as for example.Here, ρ a denotes the density of ambient air.The boundary conditions at the bottom are given by expressing the bottom stress in terms of the velocity components taken from values of the layer adjacent to the sediment-water interface.The bottom stress can be related to the turbulent law of the wall, a drag coefficient associated with quadratic velocity or using a Manning-Chezy formula such as where T b x and T b y are the bed shear stresses defined by the depth-averaged velocities as Hydrodynamics -Theory and Model www.intechopen.com The coefficient C z = h 1/6 /η is the Chezy friction coefficient, and η denotes the Manning roughness coefficient at the bed.Thus, using the free surface equation (3) and the boundary conditions ( 4) and ( 6), and after standard approximations on convective terms, the two-dimensional vertically averaged system of shallow water equations rewritten in conservative form is obtained as 8) where U and V are the depth-averaged horizontal velocities in x-andy-direction given by Note that the system (8) has been widely used in the literature to model physical phenomena of water flows such as flood waves, dam breaks, tidal flows in an estuary and coastal water regions, and bore wave propagation in rivers, among others.

Equations for free-surface flow with temperature distribution
The starting point for the derivation of the free-surface flow model with temperature distribution is the three- where the variables have the same meaning as in (1).In (9), the force term F is given according to the Boussinesq approximation as where α is the thermal expansion coefficient and T ∞ is the reference temperature.In addition, assuming that the pressure is hydrostatic, the momentum equation in the vertical direction (9d) reduces to the following form Integrating vertically the continuity equation (9a) and using the kinematic condition at the free surface leads to the free-surface equation as presented above in (3).

259
Lattice Boltzmann Simulation for Shallow Water Flow Applications www.intechopen.comThus, using the free surface equation (3) and the boundary conditions ( 5) and ( 6), and after standard approximations on convective terms, the two-dimensional vertically averaged system of shallow water equations rewritten in conservative form is obtained as where g ′ = g (1 + αT ∞ ), Θ is the depth-averaged temperature, U and V are the depth-averaged horizontal velocities in x-andy-direction given by Numerical treatment of the equations ( 8) or ( 12) often presents difficulties due to their nonlinear form, hyperbolic nature of the homogeneous system, and presence of complex source terms particularly the differential terms involving irregular topography.Therefore, the treatment of topography and friction source terms is of major importance in many practical applications of shallow water models.These could be a source of numerical instability and may produce nonphysical oscillations mainly because dicretizations of the flux and source terms are not well balanced in their reconstruction Bermúdez & Vázquez (1994); Toro (2001); Vázquez-Cendón (1999); Vukovic & Sopta (2002).
The shallow water equations (8) or (12) have to be solved in a bounded spatial domain with smooth boundaries endowed with given initial and boundary conditions along with a prescribed bed elevation.In practice, these conditions are problem dependent and their discussion is postponed until Section 4 where numerical examples are discussed.

Equations for pollutant or temperature dispersion
To model solute or heat transport by water flows, equations ( 8) or ( 12) are coupled to the depth-averaged convection-diffusion equation of the form where Θ is the depth-averaged pollution concentration or temperature, Q is the depth-averaged source, and ν C is the diffusion coefficient.In practical situations the eddy viscosity ν C and the eddy thermal diffusivity coefficients depend on water temperature, water salinity, water depth, flow velocity, bottom roughness and wind, compare Bartzokas (1985); LaCasce & Mahadevan (2006) for more discussions.For the purpose of the work presented in this Chapter, the problem of the evaluation of eddy diffusion coefficients is not considered.

The Lattice Boltzmann Method (LBM)
The central idea of the LBM is the discretization of the kinetic equation formulated for a two-dimensional geometry as where f i is the particle distribution function which denotes the number of particles at the lattice node x =( x, y) T and time t moving in direction i with velocity e i along the lattice ∆x = ∆y = e i ∆t connecting the nearest neighbors and N is the total number of directions in a lattice.In ( 14), J i represents the collision term and F i includes the effects of external forces.
Using the single time relaxation of the Bhatnagar-Gross-Krook (BGK) approach Bhatnagar et al. (1954), the discrete Boltzmann equation takes the form where τ f is the relaxation time and f eq i is the equilibrium distribution function.In the current work the D2Q9 square lattice model Qian et al. (1992), as depicted in Figure 2 is considered.The nine velocities e i in the D2Q9 lattice are defined by where c = ∆x/∆t = ∆y/∆t.Here, ∆t is chosen such that the particles travel one lattice

261
Lattice Boltzmann Simulation for Shallow Water Flow Applications www.intechopen.comspacing during the time step.The corresponding weights w i to the above velocities are The choice of relaxation time, τ f , and equilibrium distribution function, f eq i , in (15) depends on the macroscopic equations under study.Next, the formulation of these parameters for the shallow water equations ( 8) or ( 12) and the convection-diffusion equation ( 13) are described.

Lattice Boltzmann discretization for the shallow water equations
For the shallow water equations ( 8) or ( 12), the equilibrium distribution function f eq i depends on the water depth h and the velocity field u =(U, V) T which are recovered by For the D2Q9 lattice, the equilibrium function f eq i in ( 15) is defined as Dellar (2002); Salmon (1999a) with the weight factors w i in ( 16).It is easy to verify that the local equilibrium function satisfies the following conditions where I denotes the 2 × 2 identity matrix.The central idea in the LBM lies essentially in the recovery of the macroscopic flow behaviour from the mesoscopic flow picture of the particle movement Salmon (1999a).
After discretization, equation ( 15) can be written as where τ f is the relaxation time for the flow simulation and F represents the force term in the shallow water equations ( 8) or ( 12), for example 262 Hydrodynamics -Theory and Model www.intechopen.com By applying Taylor expansion and the Chapman-Enskog procedure to equation ( 20), it can be shown that the solution of the discrete lattice Boltzmann equation ( 20) with the equilibrium distribution ( 18) results in the solution of the shallow water equations ( 8) or ( 12).For details on this multi-scale expansion, the reader is referred to Dellar (2002); Salmon (1999a); Zhong et al. (2005).
In this LBM implementation, the relaxation time is determined by the physical viscosity in ( 8) and the time step through the formula In the lattice Boltzmann method, equation ( 20) is solved in two steps: collision and streaming.
In the collision step, the equations for each direction are relaxed toward equilibrium distributions with forcing Then, at the streaming step, the distributions move to the neighboring nodes

Lattice Boltzmann discretization of the convection-diffusion equation
The LBM for convection-diffusion equation ( 13) is derived using a similar approach as the one used for the shallow water equations ( 8) or ( 12).Hence, starting from equation ( 15) and using the D2Q9 lattice from Figure 2, a lattice Boltzmann discretization of the convection-diffusion equation is where g i is the distribution function, τ g is the relaxation time and Q i is the source term associated with the convection-diffusion equation ( 13).In (23), g eq i is an equilibrium distribution function satisfying the following conditions 8 To process equation ( 23), a relaxation time and equilibrium function are required.For the convection-diffusion equation, the equilibrium function is given by where the lattice weights w i are defined in ( 16).For this selection, the source term in ( 23) is set to It should be noted that the convection-diffusion equation ( 13) can be obtained from equation ( 23) using the Chapman-Enskog expansion.Details on these derivations were given in Banda et al. (2009); Klar et al. (2008); Thömmes et al. (2007).

263
Lattice Boltzmann Simulation for Shallow Water Flow Applications www.intechopen.com The relaxation time is defined by the diffusion coefficient in (13) as well as the time step Notice that conditions ( 22) and ( 27) establish a relation between the diffusion coefficient and the relaxation time used in the LBM simulations.

Implementation and boundary conditions
The computational domain is discretized by a square lattice using the D2Q9 model with 9 velocities as shown in Figure 2. The scheme to advance the solution from the time t n to the next time t n+1 can be implemented based on the following steps: Step 1.
Flow boundary conditions for the height, h, and/or the velocities, (U, V), are needed at the inlet and the outlet of computational domains.When the height h l is prescribed at the left boundary, the three distributions f 1 , f 5 and f 8 are unknown.The techniques described in Zhou (2002); Zou & He (2002) for flat interfaces to implement these boundary conditions in 264 Hydrodynamics -Theory and Model www.intechopen.comthe framework of LBM can be used.Assuming that V = 0, the velocity in x-direction can be recovered from the relation and the unknown distributions are defined as Neumann boundary conditions are implemented by imposing the equilibrium distribution corresponding to the prescribed height, h l , and the velocity of the nearest neighbor in direction of the normal, (U n , V n ) Dirichlet boundary conditions for a prescribed concentration/temperature Θ 0 can be imposed by the equilibrium for the unknown populations Neumann boundary conditions in convection-diffusion problems are implemented in a similar way by prescribing the concentration of the neighbour node Θ n at the boundary For more details on implementation of general boundary conditions in LBM, the reader is referred to Gallivan et al. (1997); Klar et al. (2008); Zou & He (2002) and further references therein.General details on the implementation of an LB method for irregular domains can also be found in der Sman & Ernst (2000); Mei et al. (1999) among others.

Numerical examples and results
In this section a practical study of the LB method to shallow water problems on complex geometry and irregular bathymetry is presented.The bathymetry is given either by an analytical function or by data points in a two-dimensional domain.The aim of this section is to test the accuracy, efficiency and study challenges for the LB approach for practical situations.Numerical results for several test cases will be presented.To verify this approach, the problem of mean flow in the Strait of Gibraltar has also been used as a test example.This latter example presents a challenge to numerical schemes because of its irregular geometry, complex bathymetry and presence of bottom friction and wind stresses.The results obtained are competitive in comparison with other approaches that solve the macroscopic equations using direct discretization methods.They are obtained without consideration for well-balancing, or adaptive grids and other technical details as is the case with other approaches.The term well-balanced schemes refers to those methods that require special treatment of the source terms such that the discretization of the flux gradients is balanced with the one used for the source terms.For more details on well-balanced schemes for shallow water equations

265
Lattice Boltzmann Simulation for Shallow Water Flow Applications www.intechopen.comthe reader is referred to LeVeque (1998); Vázquez-Cendón (1999); Xing & Shu (2006), while references on adaptive methods for solving shallow water equations can be found in Ambrosi (1999), among others.The findings in this section inform applied scientists to consider the LB method as an alternative practical numerical scheme for solving flow problems modelled by the shallow water equations.For all the results presented in this section, the gravity acceleration is set to g = 9.81 m/s 2 , the relaxation times τ f and τ g are fixed and time steps are selected according to conditions ( 22) and ( 27).

Verification of the method using the hydraulic model
To verify the performance and accuracy of the LBM, the one-dimensional shallow water equations with known analytical solutions are considered.In the current simulations, the bed friction, Coriolis forces and wind stresses are neglected in equations ( 8).It should be pointed out that a two-dimensional LBM code has been used to reproduce numerical solutions for the one-dimensional problems.Therefore, boundary conditions in the y-direction have to be supplied for the two-dimensional code.For these test examples, the dimension in y-direction is fixed to 50 lattice points and periodic boundary conditions are assumed on the upper and lower walls.Thereafter tests will be performed on flow through the Strait of Gibraltar.

Lake-at-rest example
The benchmark problem of a lake at rest proposed in Bermúdez & Vázquez (1994) to test the conservation property of numerical methods for shallow water equations is solved.The lake bed is irregular, so this test example is a good illustration of the significance of the source term treatment for practical applications to natural water-courses.It is expected that the water free-surface remains constant and the water velocity should be zero at all times.The LBM is executed using τ f = 0.6, c = 200 m/s and the results are displayed at the time t = 10800 s.
Figure 3 shows the free water surface along with the lake bed.In this figure a close-up of the free surface is included for better insight.As can be seen, small perturbations appear on the free surface.The amplitude of these perturbations decreases as the number of lattice points increases.A comparison of the relative errors on a sequence of grids reveals nearly first-order convergence of the water level h + z (see Table 1).The relative error is computed using the maximum norm where H ref = 16 m is the reference height of the water surface in this case.On a mesh with 320 lattice points along the channel length, the amplitude of perturbations is of the order of 0.1 m.
The error of 0.3 m on the grid with 80 nodes seen in Figure 3 corresponds to a relative error of 2%, and one can argue if this is sufficiently small.In general, on the one hand, the LBM is a simple scheme that is easily implemented but one needs a finer grid, while more sophisticated schemes can use a coarser grid to achieve the same accuracy.On the other hand, one wishes to simulate an application with a complex geometry where more elaborate schemes are difficult to implement, while this is not a significant problem for LBM.It can be mentioned that the performance of the LBM approach is very attractive since the computed solution remains stable and accurate even when coarse lattices are used without solving Riemann problems or reconstructing upwind fluxes or requiring complicated techniques to balance the source terms and flux gradients as those reported in Vázquez-Cendón (1999).Table 1.Relative error of the free surface h + z for the lake at rest on different meshes at time t=10800s.The convergence order obtained from a least squares fit is p = 1.09.Here δx = ∆x/L is the cell size ∆x relative to the domain length L = 1500 m.

Tidal wave flow
Secondly, the problem of a tidal wave flow in a frictionless (C b = 0) channel with length, L = 14Km is considered.The bottom topography is analytically defined by The initial conditions for the water height and velocity are h(x,0)=60.5− Z(x), u(x,0)=0.
At the channel inflow and outflow, respectively, the following Following Bermúdez & Vázquez (1994), an asymptotic analytical solution for this example can be developed as This asymptotic analytical solution is used to quantify the results obtained by the LB method.
The relative L ∞ -, L 1 -andL 2 -error norms are defined as , where e n ij = u n ij − u(x i , y j , t n ) is the error between the numerical solution, u n ij , and the analytical solution, u(x i , y j , t n ),a tt i m et n and lattice point (x i , y j ).F o rt h eL Bm e t h o d ,τ f = 0.6, c = 200m/s is used and the results are displayed at time t = 9117.5s.For this test example the ratio U/c = 0.0009.In Figure 4 the error norms for the velocity solution using four uniform lattices with sizes ∆x = ∆y = 56m, 28m, 14m and 7m are plotted.Logarithmic scales are used on the x-andy-axis.It is easy to verify that decreasing the lattice size results in a decrease of all error norms.As expected the LB method shows a first-order accuracy for this test example.The velocity values corresponding to the considered lattices are plotted along the analytical solution as shown in Figure 5. Grid convergence is clearly observed in this figure .Only a small difference between the LB solutions obtained with lattice resolution ∆x = ∆y = 7m and the asymptotic analytical solution is observed.Figure 6 presents the numerical and analytical solutions for the free surface at the simulation time t = 9117.5susing ∆x = ∆y = 7m.There is an excellent agreement between the numerical results obtained by the LB method and the asymptotic analytical solution.

Mean flow in the Strait of Gibraltar
The next application is the problem of mean flow in the Strait of Gibraltar.The schematic description of the Strait of Gibraltar is given in Figure 7.The system is bounded to the North and South by the Iberian and African continental forelands, respectively, and to the West and East by the Atlantic Ocean and the Mediterranean sea, respectively.This test problem is chosen because it presents a true practical test of lattice Boltzmann shallow water flow for two major reasons.Firstly, the Strait of Gibraltar's domain is a large-scale domain including high gradients of the bathymetry and well-defined shelf regions.Secondly, the Strait contains complex fully two-dimensional flow structures, which present a challenge in the shallow water modelling.The Strait of Gibraltar has also been the subject of numerous investigations such as water circulation, hydrodynamic processes and tides, compare Almazán et al. (1988); González & Sánchez-Arcilla (1995); Lafuente et al. (1990); Tejedor et al. (1999) among others.In all these references, the simulation domain is restricted by the Tangier-Barbate axis from the Atlantic Ocean and the Ceuta-Algeciras axis from the Mediterranean sea, see Figure 7.A schematic map of the Strait of Gibraltar is depicted in Figure 7 along with the main locations and ports.In geographical coordinates, the Strait is 35 o 45 ′ to 36 o 15 ′ N latitude and 5 o 15 ′ to 6 o 05 ′ W longitude.This domain is taken in numerical simulations mainly because measured data is usually provided by stations located on the above mentioned cities.The main objective in this numerical example is to test the capability of the LB method to handle complex geometry and irregular topography.The main astronomical tidal constituents in the Strait of Gibraltar are the semidiurnal M 2 , S 2 and N 2 tides, and the diurnal K 1 tide.Thus, in this computational study, the boundary conditions on open boundaries are prescribed as where A k is the wave amplitude, ω k the angular frequency and ϕ k the tide phase for the tide k, k = M 2 , S 2 , N 2 or K 1 .In ( 29), h 0 is the averaged water elevation set to 3m.
Initially, the flow was at rest and two weeks of real time were simulated.At the end of the simulation time the velocity fields were sampled for each tidal simulation at four different times t = 0, t = T/4, t = T/2 and t = 3T/4, where T represents the period of the considered tidal wave.
First the lattice dependence of the solutions was examined.To this end, the LB code using the M 2 tidal conditions on three different meshes with lattice sizes ∆x = ∆y = 500m, 250m and 125m was executed.In Figure 9 the cross sections of the water height at mid-width of the Strait at times t = T/4 and t = 3T/4 are shown.It is evident that, for this flow regime, the results obtained on the coarse lattice of 500m show differences to those obtained on the fine lattice of 125m.These differences noticeably decrease for the lattice of 250m.For instance, the discrepancies in the maximum water height on the lattices with sizes 250m and 125m are less than 1.92% and 2.34% at t = T/4 and t = 3T/4, respectively.Similar results, not reported here, were obtained for the water velocity and for the other tidal waves.Therefore, bearing in mind the relatively small differences on the results from a lattice with size 125m and 250m at the expense of rather significant increase in the computational costs, the lattice with size 250m  was believed to be adequate to obtain reasonable results subject to minimal lattice effects.Hence, the results presented herein are based on the mesh with lattice size ∆x = ∆y = 250m.
The computed velocity fields using the parameters of the semidiurnal M 2 , S 2 and N 2 tidal waves are presented in Figure 10, Figure 11 and Figure 12, respectively.The results for the diurnal K 1 tidal wave are presented in Thömmes et al. (2007).The results at four different times using the corresponding time period of each tide are displayed.Once the period is completed, dynamics of the water flow are repeated reproducing analogous velocity fields.The results show different aspects in the flow generated using tidal conditions for the semidiurnal M 2 , S 2 and N 2 tidal waves.Using the conditions for the semidiurnal tides, the flow exhibits a recirculating zone of different magnitudes near the Craminal Sill.At later time, before the period is completed, the flow generated by semidiurnal M 2 , S 2 and N 2 tidal waves changes the direction pointing towards the Atlantic Ocean.A recirculating flow region is also detected on the top eastern exits of the Strait near Algeciras.Similar features have also been reported in Almazán et al. (1988); González & Sánchez-Arcilla (1995).The lattice Boltzmann shallow water model performs well for this test problem since it does not diffuse the moving fronts and no spurious oscillations have been observed near steep gradients of the flow field in the computational domain.It can be clearly seen that the complicated flow structures on the Caraminal Sill and near Tarifa narrows and Tangier basin are being captured by the LB method.In addition, the presented results clearly indicate that the method is suitable for the prediction of mean flow in the Strait of Gibraltar.
Finally, computational cost, in terms of CPU seconds per time step, is 0.22s for each simulation using the M 2 , S 2 , N 2 and K 1 tidal waves.Approximately 2 × 10 6 time steps were needed to reach the real time of two weeks in a solution.All the computations were performed on a Pentium IV 2.66 GHz having 1Gb of RAM.Considering the computational cost and the accuracy achieved, the LB algorithm can be considered as a competitive alternative to the 273 Lattice Boltzmann Simulation for Shallow Water Flow Applications www.intechopen.comfinite volume methods widely used in the literature to perform numerical studies on shallow water flows, in terms of both numerical accuracy and computational cost.

Application to pollutant transport
Recently, the authors in Banda et al. (2009) extended this method to pollutant transport by the shallow water flows.The behaviour of the pollutant is investigated especially in connection with a non-flat topography and the surface stress originated by the shear of blowing winds.Firstly, the accuracy and convergence features of the LB method are verified.Finally, the LB method is applied to the simulation of a contamination event taking place in the Strait of Gibraltar.

Pollutant transport in a squared cavity
A problem of convection-diffusion of a pollutant transport in a 9000 m × 9000 m squared cavity with bottom slopes given by Komatsu et al. (1997), uniform flow velocities u 1 = u 2 = 0.5 m/s are imposed as well as uniform flow water depth (h + Z) as initial condition.The initial condition for the pollutant concentration is given by the superposition of two Gaussian pulses centered, respectively, at (x 1 = 1400 m, y 1 = 1400 m) and (x 2 = 2400 m, y 2 = 2400 m),

275
Lattice Boltzmann Simulation for Shallow Water Flow Applications www.intechopen.comwhere C 1 = 10, C 2 = 6.5 and σ 1 = σ 2 = 264.For this example, the pollutant concentration is a wave that moves along the diagonal cross-section x = y with the constant speed u 1 = u 2 = 0.5 m/s.Here, the wind effects are neglected (W 1 = W 2 = 0) in the hydraulic equations, no source (Q = 0) is considered in the pollutant transport equation, and a diffusion coefficient ν C = 100 m/s 2 is used in all LB simulations.Neumann boundary conditions are used for both hydraulic variables and pollutant concentration on all walls of the cavity.We used τ = 0.01 and simulations were stopped at time t = 4600 s.At this time, the pollutant concentration reaches the end corner of the cavity.
Figure 13 shows the initial concentration and the numerical result using a uniform mesh with lattice size ∆x = ∆y = 50 m.The corresponding contour plots are presented in Figure 14.It is clear that the LB method preserves the expected transport trajectory and captures the correct dynamics.It can be remarked that due to the diffusion present in the equation of pollutant transport, the two initial pulses merge into one concentration pulse during the time process.
In addition, the obtained solutions are completely free of spurious oscillations and the moving fronts are well resolved by the LB method.
In order to check the grid dependence of the LB method for this test example, in Figure 15 cross sections of the pollutant concentration in the main diagonal (x = y)attimet = 4600 s using different meshes are presented.Four meshes with ∆x = ∆y = 200 m, 100 m,50m and 25 m. are considered.For the selected pollutant conditions, a large difference in the concentration profile is detected for the coarse mesh with ∆x = ∆y = 200 m compared with finer meshes.This difference becomes smaller as the mesh is refined.For instance, the discrepancy in the pollutant concentration on meshes with ∆x = ∆y = 50 m and ∆x = ∆y = 25 m is less than 1%.Similar trend was observed in the hydraulic variables.Therefore, bearing in mind the slight change in the results from a mesh with ∆x = ∆y = 50 m and ∆x = ∆y = 25 m at the expense of rather significant increase in computation time, the mesh ∆x = ∆y = 50 m is believed to be adequate to obtain computational results (shown in Figure 13 and Figure 15) free of grid effects for the considered pollutant transport problem.

Pollutant transport in the Strait of Gibraltar
The Strait of Gibraltar is used heavily for shipping traffic and oil cargo.As a consequence, the Strait is considered as one of the most chronically contaminated regions, see Gómez (2003).In this example the LB method will be applied to simulate a contamination event in the Strait of Gibraltar accounting for all the hydraulic effects such as friction sources, wind stresses, Coriolis forces and horizontal eddy viscosity.
Initially, the flow is assumed to be at rest and no pollutant is present i.e., where t release is the release time and D release is the release region to be located in the Strait of Gibraltar.For the continuous release, t release corresponds to the final simulation time while for the instantaneous release, t release is set to 3 hours.In this sense, the simulations are schematic, since the number, the arrangement, and the capacities of pollution sources in the Strait of Gibraltar only partially correspond to the real situation.
The flow is initialised as in Section 4.1.3.The shallow water equations ( 8) are solved without pollutant release for two weeks of real time to obtain a well-developed flow.The obtained results are taken as the real initial conditions and the pollutant is injected at this stage of computation.Depending on the wind direction, three cases are simulated namely, calm situation, eastern wind and western wind Banda et al. (2009).At the end of the simulation time the velocity fields and concentration of pollutant are displayed after 1, 3 and 6 hours from the injection of pollutant.A mesh with lattice size ∆x = ∆y = 250 m is used for all the results presented in this section.A Neumann boundary condition is used for the pollutant concentration at the open boundaries and zero concentration is imposed at the coastlines of the Strait.
First a calm situation corresponding to W x = W y = 0 m/s is simulated.The pollutant is injected in the middle of the Strait of Gibraltar and located at (5 o 53 ′ W,35 o 56 ′ N).T h e simulated results are presented in Figure 16 at three separated instants from the injection of pollutant.A simple inspection of this figure reveals that the velocity field changes the direction during the time process according to the period of the considered tides.The decrease and increase of the strengths of velocities with time can also be seen in the figure.Obviously, the spread of contaminant patch on the water free-surface is very slow for both continuous and instantaneous releases.This fact can be attributed to the small velocities generated by the tidal waves and also to the periodic character of these tides.As expected, a wider spread of contaminant patch is observed for the instantaneous release than the continuous case.
Next, a pollutant transport subject to blowing wind from the east with W x = −1 m/s and W y = 0 m/s is considered.In contrast to the previous test example, the present pollutant transport is solved with extra velocity field due to wind.The results of the simulation for the continuous release as well as those obtained for the instantaneous release are shown in Figure 17.It is clear that the pollutant transport is influenced significantly by the action of the wind.The figure shows that the proposed LB method accurately reproduces concentration fronts.Moreover, the steep gradient in the shallow water flow and the high concentration in the convection-diffusion equation highlights the good stability and capability of the LB model to resolve pollutant transport by tidal flows.
In summary, the pollutant transport is captured accurately, the flow field is resolved reasonably well, and the concentration front is shape preserving.All these features illustrate the robustness of the LB method.
Fig. 16.Flow field (first row), concentration contours for continuous release (second row) and for instantaneous release (third row) using calm conditions at three times after the release.

Free-surface temperature in the Strait of Gibraltar
Temperature can strongly interact with hydraulics in many situations of engineering interest and neglecting its effects may have significant consequences in the overall predictions.For a discussion on the thermal effects on hydraulic flows reference is made to  Polyak et al. (1996).
Fig. 17.Flow field (first row), concentration contours for continuous release (second row) and for instantaneous release (third row) using eastern wind at three times after the release.
In the simulations the flow parameters are set as in Section 4.1.3.As discussed above a mesh with lattice size ∆x = ∆y = 250 m is used for all the results presented in this section.Depending on the wind conditions, two situations are presented namely: (i) Calm situation corresponding to (ω x = 0 m/s, ω y = 0 m/s) , (ii) Wind blowing from the west corresponding to (ω x = 1 m/s, ω y = 0 m/s).
Initially, the simulated flow has been at warm rest i.e., where Θ h = 23  sea-surface temperature is included at this stage of the simulation.At the end of the simulation time the velocity fields and temperature contours are displayed after 12, 18 and 24 hours from the inclusion of sea-surface temperature.
In Figure 18 numerical results obtained using calm wind conditions are presented.Those obtained for the wind blowing from the west are displayed in Figure 19.In these figures, the velocity field is shown and 10 equi-distributed contours between Θ c and Θ h of the temperature at the instants t = 12, 18 and 24 hours.It is clear that using the conditions for the tidal waves and the considered wind situations, the flow exhibits a recirculating zone with different order of magnitudes near the Caraminal Sill (i.e. the interface separating the water bodies between the Mediterranean sea and the Atlantic Ocean).At the beginning of simulation time, the water flow enters the Strait from the eastern boundary and flows towards the eastern exit of the Strait.At later time, due to tidal waves, the water flow changes the direction pointing towards the Atlantic Ocean.A recirculating flow region is also detected on the top eastern exits of the Strait near Algeciras.Similar flow behaviours have been also reported in Almazán et al. (1988); González & Sánchez-Arcilla (1995); Thömmes et al. (2007).
The effects of wind conditions are observed in the temperature distributions presented in Figure 18 and Figure 19.A boundary layer of high sea-surface temperatures has been detected on the Spanish coastal lines.For the considered tides and wind conditions, the buoyancy force has been seen to play a weak role in influencing the sea-surface tamperature in the Strait of Gibraltar which results in thinner mixing layers.In Figure 20 the time evolution of the water free-surface elevation at the Tarifa narrows for a time period of two weeks is displayed.As expected, the time series show two tidal periods with different amplitude and frequencies.They are in good agreement with those previously computed in Castro et al. (2004); Tejedor et al. (1999).Similar results not presented here, have been obtained at other locations in the Strait of Gibraltar.
It can be clearly seen that the complicated flow structures on the Caraminal Sill and near Tarifa narrows and Tangier basin are being captured by the LB method.In addition, the presented results clearly indicate that the method is suited for prediction of sea-surface temperature dispersion in the Strait of Gibraltar.It should be stressed that ideally, results from the temperature dispersion model should be compared with observations of real sea-surface temperatures in the Strait of Gibraltar.However, there are no available data until now to 282 Hydrodynamics -Theory and Model www.intechopen.comcarry out this work.Thus, only some hypothetical simulations have been undertaken simply to show that LB results are logical and consistent.

Conclusion
The most common two-dimensional lattice Boltzmann method using nine particle speeds arranged on a D2Q9 squared lattice was used to approximate numerical solutions to the shallow water equations.The model is simple, accurate, easy to implement, and can be used to solve both steady and unsteady shallow water problems.The method also provides a straightforward treatment of source terms without relying on complicated discretization techniques.Other source terms such as wind stresses or bed shear stresses can naturally be added to the lattice Boltzmann equation as force terms without special treatment.In this chapter the main focus is to demonstrate the ability of the lattice Boltzmann method to solve practical shallow water flows on non-flat beds with irregular bathymetry.
The efficiency of the method for predicting shallow water flows was assessed in the benchmark problems such as the tidal wave flow and steady flow over a hump.The results clearly indicate that the method captures the correct flow structures and reproduces results which satisfactorily agree with those available in the literature for the same test problems.To demonstrate the ability of the lattice Boltzmann method on complex practical shallow water problems, the method has been applied to the mean flow in the Strait of Gibraltar.The numerical results show correct physics in different test regimes.The influence of different spatial resolutions on the numerical results has also been discussed.Refined spatial models in which a larger number of total particles is used in the simulation can resolve more small-scale effects at the expense of long computational times.Nevertheless, flows in such complex domains can be computed, providing correct physics without the need for generating adaptive grids or complicated reconstruction of numerical fluxes using exact or approximate Riemann solvers.Overall the method shows reasonable accuracy while ensuring the required properties of the shallow water flows.
Furthermore the lattice Boltzmann method has been extended and tested for pollutant dispersion by shallow water flows.The mass, momentum and transport equations are obtained from the nine-velocity distributions of hydraulic flow and pollutant concentration variables.Two types of distribution functions have been developed for the hydraulic variables and pollutant concentration.Although lattice Boltzmann methods are very promising, they are still in the early stage of development and validation.More investigations are needed to explore the capability of lattice Boltzmann methods for more practical engineering applications in environmental fluid flows and species transport.For instance, extension of the current solver to complex geometries involving turbulent effects can also be of interest both for hydraulics and pollutant transport.

Fig. 1 .
Fig. 1.Vertical section of the hydraulic domain and notations.

Fig. 3 .
Fig.3.The free-surface for the lake at rest on different meshes at time t = 10800 s.

Fig. 4 .
Fig. 4. Grid convergence for the tidal wave flow at time t = 9117.5s.

Fig. 7 .Fig. 8 .
Fig. 7.The Strait of Gibraltar map in geographical coordinates along with the main locations.The considered computational domain is marked by dashed lines.

Fig. 15 .
Fig. 15.Cross sections of the pollutant concentration for different meshes.

)
Parameters of the flow are set as in the previous subsection.A diffusion coefficient of ν C = 100 m/s 2 is used for the pollutant transport.The contaminant source is implemented as an 277 Lattice Boltzmann Simulation for Shallow Water Flow Applications www.intechopen.comindicator function of the form

Fig. 19 .
Fig.19.The same as in Figure18but for a wind blowing from the west.

Fig. 20 .
Fig. 20.Time evolution of water elevation at the Tarifa narrows.
Using the concentration (and velocity from shallow water equations) at time t n , compute from (25) the equilibrium function g Zou & He (2002) (1997) 1.a.Using the water depth and velocity at time t n , compute from (18) the equilibrium function f eq i , i = 0,1,...,8.1.b.eq i , i = 0, 1, ..., 8.Step 2. Distribution functions:2.a.Calculate the distribution function f i , i = 0, 1, ..., 8, using the lattice Boltzmann equation (20) with an appropriate relaxation time τ f and impose the corresponding boundary conditions.2.b.Calculate the distribution function g i , i = 0,1,...,8, using the lattice Boltzmann equation (23) with relaxation time τ g and impose the corresponding boundary conditions.Step 3. Solution reconstructions:3.a.Update the water depth and velocity using the equations (17).3.b.Update the concentration using equations (24).Step 4. Change the time t n −→ t n+1 ,got oStep 2 and repeat until the stopping criterion is reached.The time evolution is stopped in Step 4 either when a fixed time is reached for instationary problems, or by comparing the deviation between two consecutive solutions for steady problems.The implementation of boundary conditions in the LBM has a crucial impact on the accuracy and stability of the method, seeGallivan et al. (1997);Zou & He (2002)for more discussions.When no-slip boundary conditions are imposed at walls the bounce-back rule is usually used in the lattice Boltzmann algorithm.At a boundary point x b , populations f i of links e i which intersect the boundary and point out of the fluid domain are simply reflected (bounce-back) since they cannot participate in the normal propagation step

Table 2 .
Parameters of tidal waves in the stations considered in the present study.
Millán et al. (1995)rcilla (1995)al.(1999)Polyaketal. (1996);Samelson et al. (2006);Vargas et al. (1999)and further references can be found therein.The basic circulation in the Strait of Gibraltar consists of an upper layer of cold, fresh surface Atlantic water and an opposite deep current of warmer, salty Mediterranean outflowing water, compareAlmazán et al. (1988);González & Sánchez-Arcilla (1995).The sea-surface temperatures in the Strait of Gibraltar are maxima in summer (August-September) with average values of 23-24 • C and minima in winter (January-February) with averages of 11-12 • C. The north Atlantic water is about 5-6 • C colder than the Mediterranean water, elaborate details are available inMillán et al. (1995); • C is the Mediterranean temperature and the western temperature boundary of the Strait is fixed to the Ocean temperature Θ c = 17 • C. The shallow water equations (12) are solved without temperature dispersion for two weeks of real time to obtain a well-developed flow.The obtained results are taken as the real initial conditions and the Fig. 18.Flow field (first row) and temperature contours (second row) for calm situation at three different times.From left to right t = 12, 18 and 24 hours.