## Abstract

A two-phase flow model has been developed to study three-dimensional breaking waves over complex topography, including the wave pre-breaking, overturning and post-breaking processes. The large-eddy simulation approach has been adopted in this study, where the model is based on the filtered Navier–Stokes equations with the Smagorinsky sub-grid model being used for the unresolved scales of turbulence. The governing equations have been discretized using the finite volume method, with the PISO algorithm being employed for the pressure–velocity coupling. The air–water interface has been captured using a volume of fluid method and the partial cell treatment has been implemented to deal with complex topography in the Cartesian grid. The model is first validated against available analytical solutions and experimental data for solitary wave propagation over constant water depth and three-dimensional breaking waves over a plane slope, respectively. Furthermore, the model is used to study three-dimensional overturning waves over three different bed topographies, with three-dimensional wave profiles and surface velocities being presented and discussed. The overturning jet, air entrainment and splash-up during wave breaking have been captured by the two-phase flow model, which demonstrates the capability of the model to simulate free surface flow and wave breaking problems over complex topography.

## 1. Introduction

Wave breaking plays an important role in the air–sea interactions, surf zone dynamics, near shore sediment transport, wave–structure interactions, and marine renewable energy development. Over the last three decades, significant advances have been made in the theoretical, experimental and numerical studies of the characteristics of breaking waves (see [1–4] for reviews). However, little attention has been paid to investigate three-dimensional overturning waves over complex topography due to wave nonlinearities and complexity. The challenge of modelling breaking waves is further increased in the post-breaking process, with air entrainment and jet-splash cycles [5].

There are essentially two types of mathematical models for modelling overturning waves. One is the fully nonlinear potential flow model based on Laplace's equation with inviscid and irrotational assumptions, which can simulate breaking waves and can provide insight into the kinematics and dynamics of water waves up to wave overturning. Much research has been carried out for two-dimensional overturning waves [6–12]. Few attempts have been made with potential flow models to simulate three-dimensional overturning waves in deep water in a periodic domain [13,14], and three-dimensional overturning waves in shallow water over a sand bank [15], a sloping ridge [16–18] and a non-symmetrical seabed or multiple reefs [19]. However, the potential flow models, which are very efficient in computation, usually terminate before the plunging jet touching down and cannot provide any information after wave breaking in turbulent flows.

The other mathematical model for modelling overturning waves is based on the Navier–Stokes equations with a free surface calculation. In order to track or capture the interface, several techniques [20] have been employed to study breaking waves, such as the marker and cell method [21,22], surface tracking method [23,24], volume of fluid (VOF) method [25–32], level set method [33,34], the density function method [35] and some meshless methods [36–40]. Among them, there are one-phase flow models, in which only the flow in water is considered in the computation, the pressure in the air is taken as a constant, and the boundary conditions are specified at the free surface. During wave breaking, these one-phase flow models may be inadequate to deal with the air entrainment and splash-up process. Additional complication arises in the treatment of boundary conditions at the highly distorted free surface. Thus, in order to take the air into account for wave breaking, several two-phase flow models, in which both flows in the air and water are solved, have been developed to study the details of two-dimensional breaking waves and the air entrainment during wave breaking in deep water [26,33,34], two-dimensional breaking solitary waves on sloping beaches [28,32,41] and two-dimensional breaking waves in the surf zone [42–44]. With increases in computational power and developments of numerical methods, recently some attempts have been made in the development of three-dimensional two-phase flow model to investigate breaking waves over a double reef [45], deep-water periodic breaking waves over constant water depth [30], wave–body interactions [46] and turbulence structure in breaking waves over a plane slope [47]. In addition, the coupling between the potential flow model and the Euler equations has been applied for breaking waves over a plane slope and a sloping ridge [48]. However, numerical simulations of three-dimensional breaking waves over complex topography are still limited.

The objective of this study is to investigate the kinematics and dynamics of three-dimensional breaking waves over complex topography, including the wave pre-breaking, overturning and post-breaking processes. In this study, the two-phase flow model [32], which solves the flow in the air and water simultaneously, has been further developed and extended for modelling three-dimensional breaking waves and providing detailed phenomena during wave overturning. The key aspect of the model is that in order to deal with complex topography, the partial cell treatment has been extended in three dimensions and introduced in the finite volume discretization together with an implicit time stepping, which prevents the instability in small cut cells without the commonly used cell-merging technique in the literature. The organization of this paper is as follows. The description of the mathematical model for the two-phase flow is described in §2. The numerical method and implementation are presented in §3. The applications of the three-dimensional two-phase flow model are shown in §4. The numerical model is first validated by simulating solitary wave propagation over constant water depth, in which numerical results are compared with the analytical solutions and the convergence property is checked. Breaking waves over a plane slope are then calculated and both quantitatively and qualitatively compared with the experimental measurements. Finally, three-dimensional overturning waves over three different bed topographies are presented. The paper ends with the conclusion in §5.

## 2. Mathematical model

### (a) Governing equations

The large-eddy simulation (LES) approach is adopted in this study, and the governing equations used for incompressible two-phase flow are based on the filtered Navier–Stokes equations, given as
*x*,*y*,*z*), *t* is the time, ** g** is the gravitational acceleration vector,

*ρ*and

*μ*are the density and dynamic viscosity of the fluid.

The term ** I** is the unit tensor and

*μ*

_{t}is the turbulent eddy viscosity defined as

*C*

_{s}=0.1 used in this study. The symbol for spatial filtering ‘

^{−}’ is dropped hereinafter for simplicity.

The momentum equation is closed with the constitutive relations for the density and dynamic viscosity of the fluid as given by
*F* is the volume fraction defined as

The air–water interface is then within the cells where 0<*F*<1. A particle on the surface stays on the surface and the volume fraction *F* has a zero material derivative

### (b) Initial and boundary conditions

In order to completely describe the mathematical model, it is necessary to define the boundary conditions in a computational domain. On the solid boundary, the wall model [50] is used for the near-wall treatment in LES modelling. For the outlet, the zero-gradient or radiation boundary conditions are applied for the flow. As both fluids in air and water are solved simultaneously in the present two-phase flow model, the kinematic and dynamic free surface boundary conditions are already implemented and they do not need to be specified explicitly at the air–water interface. The inlet boundary conditions are case dependent. In general, the time history of the velocity field and the volume fraction at the inlet are obtained from an analytical solution of water waves.

In the computation, the initial flow field at *t*=0 has to be prescribed. For calculations with the fluids initially at rest, the flow field is initialized with zero velocity and hydrostatic pressure, and the volume fraction is computed from the initial water depth. When the wave is initialized in the computational domain, the water velocities and water surface are specified using the corresponding analytical solution for water waves. The velocity in the air is initialized as zero in this case as little is known about the flow in the air and the pressure distribution in the whole domain is hydrostatic.

## 3. Numerical method

In order to study three-dimensional breaking waves over complex topography, partial cell treatment without the instability issue has been implemented to deal with complex geometries in the Cartesian grid. In this study, the finite volume method (FVM) is used to discretize the governing equations on a staggered Cartesian grid, which has the advantage of strong coupling between the velocity and the pressure. Figure 1*a* shows a typical variable arrangement in a three-dimensional Cartesian grid, in which the velocities are located on the face centre of the control volume; the pressure, all other scalar variables and the volume fraction *F* are stored at the cell centre. Figure 1*b* shows a typical control volume used in this study, in which P is the present node, the upper-case letter E, W, N, S, B and R denote neighbouring nodes on the east, west, north, south, back and front with respect to the central node P. The lowercase e, w, n, s, b and r denote the corresponding face of the control volume, whereas c denotes the centre of the control volume.

### (a) Finite volume discretization

In the FVM, also known as the control volume method, the whole domain is divided into a number of control volumes, such that there is a control volume surrounding each grid point. The differential equation is integrated over each control volume in order to derive the algebraic equation containing the grid-point values of *ϕ*, where *ϕ* is the considered variable. The discretized equation expresses the conservation principle for a finite control volume, just as the differential equation expresses it for an infinitesimal control volume. The FVM is conservative and can deal with complex geometries [51,52], thus it is especially suitable for modelling free surface flows due to the requirement of mass conservation and the deformed interface, therefore it is adopted in this study.

Consider a volume of fluid *Ω* which has an arbitrary domain, the surface of the control volume is *S* and the unit outward normal vector to the face f is ** n**. All the governing equations can be recast into a general integral formulation as below

*ϕ*denotes the dependent variable,

*Γ*is the viscosity and

*Q*

^{S}

_{ϕ}is the source term in the control volume.

Table 1 shows the various values of *ϕ*, *Γ* and *Q*^{S}_{ϕ} in the general integral formulation to represent the filtered Navier–Stokes equations. It is noted that the final form of the continuity equation (2.1) used here is obtained under the assumption that the fluid is incompressible.

### (b) The complex geometry treatment in Cartesian grids

To deal with complex geometries in engineering applications, overlapping grids, boundary-fitted grids and unstructured grids can be used. Unstructured grids provide great flexibility to conform onto complex stationary or moving boundaries, and can easily refine or coarsen meshes in a region of the domain depending on the flow feature. However, they require additional computational efforts and further complicate the algorithm implementation as there is no pre-define order of the control volumes and their geometric layouts need to be calculated. Furthermore, generating high-quality boundary-fitted or unstructured grids is usually very cumbersome [53]. Cartesian grid methods which can simulate flow with complex geometries on Cartesian grids, avoid these problems. The most popular methods are the immersed boundary method [53–56] and Cartesian cut cell method [57–62]. The primary advantage of the Cartesian grid method is that only little modification of the program on Cartesian grids is needed to account for the complex geometries. It also has the advantage of simplified grid generation and simulating flow with moving boundaries due to the use of stationary, non-deforming grids. The drawback of these methods is that implementing boundary conditions is not straightforward and instability problems may occur in small cells when explicit schemes are used. Thus, the cell-merging technique [63] and using slightly different control volumes [64] are developed to avoid this instability, both of which effectively increase the size of the cut cell.

As most previous studies investigated breaking waves over a simple geometry, the main objective of this study is to investigate the effect of complex topography on breaking waves. In this study, the partial cell treatment is used and a *θ* function in a control volume, arising from the Fractional Area–Volume Obstacle Representation method [65], is extended to three dimensions and introduced in the finite volume discretization. The *θ* function is defined whose value is 1 for a point accessible to fluid and 0 for a point inside an obstacle. The average of this function over a control volume or cell face is the fraction of the volume or area available to the flow. Figure 2 shows a typical cut cell in a three-dimensional Cartesian grid. It is worth mentioning that, to prevent the instability in small cells, an implicit scheme for time integration is employed for the governing equations together with the finite volume discretization in the present method, which is different from the existing literature mentioned above.

In contrast to a full cell, the convective and diffusive fluxes at cell faces are modified in a cut cell, which will be presented in the following discretization.

### (c) Spatial discretization

#### (i) Convection term

The finite volume discretization of the convection term is obtained as
*A* is the area of the face and *m*=*ρ*** u**⋅

*n**θA*is the mass flux through the face.

In cut cells, the mass flux has also to be modified by the *θ* function on the boundary. If *θ*=0, there is no mass flux through the face and the convective flux is obtained as *m*_{f}=0.

The mass flux at the faces of the momentum control volume can be obtained by the interpolation of values of *ρ* and ** u**, such as

*m*

_{f}=

*ρ*

_{f}

*u*_{f}⋅

*n**θ*

_{f}

*A*

_{f}, however, the mass conservation in the momentum control volume can be only guaranteed to the accuracy of the interpolation procedure [51]. Thus, in this study, the

*m*

_{f}is obtained from the interpolation of the mass fluxes, which is already available at the faces of the continuity control volumes. The face value

*ϕ*

_{f}can be obtained from different schemes and a high-resolution scheme [52], which is based on the first-order upwind scheme with high order terms through flux limiting, is used in this study. It combines the high order accuracy with monotonicity and more details can be found in [32].

#### (ii) Diffusion term

The finite volume discretization of the diffusion term is obtained as
_{e}=|eE|/|PE|. Analogous expressions can be derived for all other faces ( f=w, n, s, b, r) by making appropriate index substitutions and will not be shown here.

The gradient at the face is calculated by the finite difference approach as
_{Pnb} is the distance from the present point P to the neighbouring point nb.

When the control volume is a cut cell, special attention has to be paid to the spatial discretization. When the face of a momentum control volume is on the wall, the diffusion flux is obtained as
*ϕ*/∂*n* is calculated by the finite difference approach in equation (3.5) and *τ*_{w} is the shear stress of the wall on the face of the control volume.

#### (iii) Source term

The finite volume discretization of the source term is obtained as
*ρ*_{c} and *θ*_{c}) is obtained by the volume averaging of two values on the face of the control volume and the pressure gradient is calculated as

### (d) Temporal discretization

A backward finite difference is used for the time derivative, which leads to an implicit scheme for the governing equations
*t* is the time step and the superscripts *n*+1 and *n* mean the value in current and previous time step, respectively. The implicit scheme has the advantage of unconditional stability and thus can prevent the instability problem in small cut cells.

### (e) General form of the discretization

Substituting all the above discretized terms into equation (3.1) and subtracting the continuity equation ∂*ρ*/∂*t*+**∇**⋅(*ρ*** u**)=0 multiplied by

*a*

^{ϕ}is the coefficient, the subscripts P and nb=E,W,N,S,B,R denote the variables in the present and neighbouring cells, respectively, and

The coefficients depend on the approximations used and the first-order upwind scheme is used in this study as the basis of the formulation; the high-resolution scheme is implemented using the deferred correction method via source term *Q*^{H}_{ϕ} [51]. For example, the coefficients for the momentum equation *ϕ*=*u*,*v*,*w* are

The algebraic equations are solved by the Strongly Implicit Procedure method or Bi-Conjugate Gradients Stabilized method in this study.

### (f) Pressure–velocity coupling

In the incompressible Navier–Stokes equations, pressure and velocity are decoupled as the pressure term does not appear in the continuity equation. For some numerical discretizations, this may cause convergence problems. However, when a staggered mesh is used, as in this work, coupling occurs as a result of the discretization, as velocity updates on cell faces contain pressure terms. In this study, the PISO algorithm [67] is employed for the pressure–velocity coupling, and more details can be found in [68]. The PISO algorithm is used to calculate the corrected pressure twice, and after solving the pressure correction equations the updated pressure and velocity are added by the pressure and velocity correction terms, respectively.

### (g) Volume of fluid scheme for interface capturing

A key requirement for modelling two-phase flow is a method of calculating the shape of the interface. Numerous methods have been proposed and used for the simulation of interfacial flows. However, the VOF method for tracking the interface is one of the most popular approaches due to its advantages: mass conservation, computational efficiency and easy implementation. From a general point of view, there are two classes of algorithms to solve the *F* transport equation (2.8): algebraic and geometric computation [69]. Excellent reviews on the VOF methods can be found in [20,69].

Considering the advantages of the VOF method and efficiency in algebraic computation, the high-resolution VOF scheme Compressive Interface Capturing Scheme for Arbitrary Meshes (CICSAM) is employed in this study to capture the air–water interface for breaking waves. CICSAM is a high-resolution scheme based on the normalized variable diagram used by Leonard [70]. It contains two high-resolution schemes and the weighting factor is based on the angle between the interface and the direction of motion. An outline of CICSAM is given below. Refer to [71] for the details.

The normalized variable *V*_{f} is the volumetric flux. However, the Hyper-C scheme is inadequate to preserve the shape of an interface which lies tangentially to the flow direction. Thus CICSAM switches to the ULTIMATE-QUICKEST (UQ) scheme [70]

Thus, depending on the angle between the interface and the flow, CICSAM combines these two schemes, then
*k*_{γ} is a constant introduced to control the dominance of the different schemes and the recommended value is *k*_{γ}=1, *α*_{γ} is the angle between the vector normal to the interface and the vector which convects the centres of donor and acceptor cells.

The final expression for the face value of *F* is
*β*_{f} is obtained by

It is noted that the normalized variable in equation (3.12) will be divided by zero if the volume fraction *F* has the same value in the acceptor and upwind cell. In the numerical implementation, the numerator and denominator of the weighting factor in equation (3.18) are multiplied by (*F*_{A}−*F*_{U}), resulting in a modified expression of the normalized variable on the face (not shown here), to avoid the singularity in the computation.

## 4. Results and discussion

### (a) Solitary wave propagation over constant water depth

In order to validate the three-dimensional two-phase flow model, solitary wave propagation over constant water depth [25] is investigated in this section, which has been considered as a classical benchmark problem for a wave simulation model. The analytical solution for the solitary wave profile is given as
*H* is the wave height, *D* is the constant water depth, *x*_{L} is the initial location of the wave.

In the numerical simulation, the constant water depth *D*=1 m is used, with the ratio of wave height to still water depth, *H*/*D*=0.1 [25], being considered in this study. The computational domain is 100 m long, 1 m wide and 1.4 m high. It is discretized by a uniform grid including 1000×3×140 points in the streamwise, spanwise and vertical directions, respectively. The turbulence model is switched off and the dynamic viscosity is set to zero, so that the potential flow is solved by the model, which is the assumption of the analytical solution for the solitary wave. At the inlet, the solitary wave is generated by giving the water surface profile with equation (4.1) and the water particle velocities based on the analytical solution [72], whereas the radiation boundary condition is applied at the outlet and free-slip elsewhere. At *t*=0 s, the velocities in water are initialized from the analytical solution and the velocities in air are initialised as zero. The simulation is run for *t*=30 s with a constant time step Δ*t*=0.0025 s.

Figure 3*a* show the comparison between the predicted wave profiles and the analytical solutions at different times along the centreplane. It can be seen that the numerical results are very close to the analytical solutions, although there is a slight phase shift after the solitary wave propagating about 100 water depth. Figure 3*b* show the predicted results of the time histories of normalized mass and energy in water, including both the kinetic and potential energy. When the solitary wave propagates into the domain, the mass of the water is nearly kept constant, demonstrating good conservation property for the interface capturing method used in the present model. It can also be found that the energy is conserved during the simulation, as there is no turbulence and viscous effects in this case. From *t*=12 s to *t*=24 s when the solitary wave is completely within the computational domain, the error for the mass is less than 0.01% and the error for the total energy is 0.51%.

It is shown from the above that good agreement is obtained between the predicted surface profiles and analytical ones. Nevertheless, it does not mean that the internal kinematics of water waves are correctly simulated in the model. In order to demonstrate this, the comparison between the predicted velocity field and the analytical solution for the solitary wave at the same location is shown in figure 4. Although the solitary wave propagating about 80 water depth, the predicted velocity field agrees well with the analytical solution in terms of the magnitude and the direction of the velocity vectors, indicating accurate internal kinematics during wave propagation are obtained by the present model.

In order to check the convergence property related to the time step, additional simulations have been performed for the same spatial resolution as mentioned above, but with three different time steps, namely Δ*t*= 0.02, 0.01 and 0.005 s. Figure 5*a* shows the comparison between the predicted wave profiles obtained using four different time steps and the analytical solution at *t*=24 s. It is shown that the numerical results generally agree with the analytical solution. The differences between the numerical results for different time steps are small provided the grid is fine enough. The numerical results for Δ*t*=0.005 and 0.0025 s are quite close to each other, indicating convergent solutions have been obtained in this study. Figure 5*c* shows the temporal convergence and it can be seen that the solution is converging on the smallest time step. There is nearly no change if we reduce the time step and the slope is smaller than the first-order temporal discretization because other discretization truncation errors are dominant at this stage.

In order to check the convergence property related to the grid size, additional simulations have also been performed using the same time step Δ*t*=0.0025 s for three different grid sizes, namely Δ*x*=0.8 m and Δ*z*=0.02 m, Δ*x*=0.4 m and Δ*z*=0.015 m, and Δ*x*=0.2 m and Δ*z*=0.0125 m. As the solitary wave is two-dimensional, the number of grid points in the spanwise direction is kept constant in all cases. Figure 5*b* shows the comparison between the predicted wave profiles for four different grids and the analytical solution at *t*=24 s. It is shown that better agreement with the analytical solution is obtained when the grid is refined. The coarsest grid, corresponding to 20 points in one wave length and five points in whole wave height, is not sufficient to resolve the wave profile during wave propagation. However, when the finest grid is used, close agreement is obtained between the predicted result and the analytical result, indicating the numerical solution nearly converges to the exact solution. Figure 5*d* shows the spatial convergence and it can be seen that it is second-order accurate at large grid sizes, corresponding to the high-resolution scheme we used in the present method. When the mesh is refined, the result is nearly converged in space and other discretization truncation errors are dominant at this stage.

### (b) Breaking waves over a plane slope

Non-breaking solitary wave propagation over constant water depth has been presented in §4a, in which good agreement with the analytical solutions has been obtained. In this section, the three-dimensional two-phase flow model is used to study breaking waves, and compare both quantitatively and qualitatively with the experiment [73] for a breaking solitary wave splash-up on a 1 : 15 sloping beach.

In the simulation, the computational set-up is set up to replicate the laboratory studies undertaken and reported by Li [73]. The schematic of the run-up of a breaking solitary wave is shown in figure 6, where the origin of the coordinate system is on the still water level above the toe of the beach along the centreline, *x*, *y* and *z* are the streamwise, spanwise and vertical coordinates, respectively. The slope of the beach is *D*=0.3048 m. Two cases for the incident solitary wave with the ratio of wave height to still water depth, namely *H*/*D*= 0.45 and 0.40, are considered in this study. The computational domain starts from the toe of the beach and extends to the location beyond the maximum run-up point 18.75*D*. The width and height of the computational domain are 1.3*D* and 1.75*D*, respectively. In order to get good results and save the computational effort, based on our experience for a two-dimensional simulation study [32], the computational domain is discretized by a 900×90 non-uniform grid in the streamwise and vertical directions, with minimum meshes Δ*x*/*D*=0.01 and Δ*z*/*D*=0.01 in the breaking region. A uniform distribution with 40 points is used in the spanwise direction. At the inlet, the solitary wave is generated by giving the water surface profile and the water particle velocities based on the analytical solution [72]. Zero-gradient boundary conditions are applied at the top and the outlet of the computational domain. The sloping beach is considered as a solid wall, whereas sidewalls are free-slip.

Figure 7*a* shows snapshots of the water surface profiles during wave overturning on a sloping beach for *H*/*D*=0.45 at *b* shows the quantitative comparison between the predicted surface profiles along the centreplane and the experimental data [73]. At *V* <*C* before wave breaking. At *C*. As a consequence, the wave has passed the breaking point, which is defined as when the front of the wave becomes vertical, and starts to overturn. At *a* that the surface velocities in most parts of the wave is less than the wave phase speed and only the surface velocities in the breaking region are greater than the wave phase speed. It can be seen from figure 7*b* that the computational results agree with the experimental measurements in terms of the wave shape and location before wave curling down, and there is only a slight difference in the size of the cavity enclosed by the plunging jet. The slight discrepancy may be caused by the solitary wave at the inlet generated from the analytical solution differing slightly from the experiment [74].

It is noted that the numerical results are only presented up to the wave overturning process in figure 7 because the measured water surface profiles are not available after the wave touching down for the case (*H*/*D*=0.45). In order to show the capability of the present model in simulating wave post-breaking and subsequent splash-up processes, the case for the breaking solitary wave with *H*/*D*=0.4 is computed with the same above-mentioned computational set-up, and qualitatively compared with the laboratory photographs [73,75] in figure 8. Although we cannot use the same viewpoint and location as used in the experimental photos, it can be seen that the development of plunging jet, jet impingement and generation of splash-up are generally reproduced in the three-dimensional simulation, which reasonably agree with the experimental measurements. However, there is a discrepancy between the simulation and experiment for the enclosed cavity and the secondary jet generated during the splash-up process. This is attributed to the aeration and compressibility of the fluid, which have been neglected in the present model, and may play important role during the post-breaking process. Thus, the present model has limitations to simulate the strong two-phase flow mixing and the generation of small bubbles during the splash-up process, which is a complex multiphase flow phenomena.

### (c) Breaking waves over complex topography

The three-dimensional two-phase flow model has been validated against analytical solutions and experimental measurements for water waves over constant water depth and a sloping beach, respectively, in §4a and b. This section aims at investigating characteristics of three-dimensional overturning waves over complex topography. Three bed topographies have been considered in this study, namely the three-dimensional sloping ridge [16,19], the three-dimensional sloping trough and the three-dimensional sloping asymmetry, which can be considered as one or a combination of examples of convex and concave coastlines. The schematic of breaking waves over complex topography is shown in figure 9, where the origin of the coordinate system is on the still water level above the toe of the bed along the centreline, *x*, *y* and *z* are the streamwise, spanwise and vertical coordinates, respectively. The bed topography is expressed as
*s*_{b}=1/15 is the slope of the bed along the centreline and sech^{2}(*k*_{b}*y*/*D*) is the spanwise modulation of the bed topography with coefficient *k*_{b}=0.5.

In the simulation, the constant water depth *D*=1 m is used, with the ratio of wave height to still water depth, *H*/*D*=0.6, and the solitary wave crest is initialized at *x*/*D*=0. The computational domain of 18*D*×4*D*×2*D* is discretized using a grid of 575×40×70 points in the streamwise, spanwise and vertical directions, respectively. The grid is uniform in the spanwise direction and non-uniform in the streamwise and vertical directions, with minimum meshes Δ*x*/*D*=0.02 and Δ*z*/*D*=0.02 in the breaking region. The boundary conditions are the same as used in §4b.

Figure 10 show a sequence of snapshots of the predicted water surface profiles for three-dimensional breaking waves over complex topography, in which the water surfaces are coloured based on the local values of the normalized velocity magnitude *V*/*C*.

Before wave breaking at *y*/*D*=2 in the asymmetry case, respectively. At

Compared to the numerical results presented in §4a,b, it is shown that three-dimensional wave profiles are developed during wave overturning over complex topography, which would better represent the actual wave breaking process seen in nature. The three-dimensional breaking waves can break first at some points and continue with spreading laterally, depending on the local topography. It can be seen from figure 10 that the effect of seabed changes the wave shape and hydrodynamics associated during the wave breaking process. Air entrainment and jet-splash cycles can be observed during wave overturning and post-breaking processes, respectively. Both symmetric and asymmetric three-dimensional overturning waves have been presented here, demonstrating the capability of the present model in simulating breaking waves over complex topography.

## 5. Conclusion

In this paper, a two-phase flow model has been developed for investigating three-dimensional breaking waves over complex topography, focusing on the wave overturning and post-breaking processes. The LES approach has been adopted in this study, where the model is based on the space filtered Navier–Stokes equations with the Smagorinsky sub-grid model being employed for the unresolved scales of turbulence. The FVM has been used to discretize spatial derivatives with the PISO algorithm being employed for the pressure–velocity coupling, and a backward finite difference discretization has been used for the time derivative, which lead to an implicit scheme for the governing equations. The air–water interface of breaking waves has been captured using the high-resolution VOF scheme CICSAM, with the partial cell treatment being implemented to deal with complex three-dimensional topography.

In order to validate the model, solitary wave propagation over constant water depth was first computed, with the water surface profiles and internal kinematics being compared with the analytical solutions. It is shown that good agreement between numerical and analytical solutions, as well as good conservation property for the mass and energy were obtained in this study, with the property of convergence being assessed. Furthermore, breaking waves over a plane slope were calculated, where both quantitative and qualitative comparisons were made between numerical simulations and corresponding experimental data. The main features during wave pre-breaking, overturning and post-breaking processes were generally reproduced in the simulation, although there was a discrepancy between the simulation and experiment for the enclosed cavity and the secondary jet generated during the jet-splash cycles. Finally, the three-dimensional two-phase flow model was used to study three-dimensional breaking waves over three different bed topographies, namely a sloping ridge, trough and asymmetry. It was found that the waves developed different three-dimensional shapes during wave overturning, depending on the local bed topography. The three-dimensional wave profiles with surface velocities during wave breaking were presented in the simulation, including the wave post-breaking process, which was ignored in previous studies with potential flow models.

This study demonstrates the capability of the present two-phase flow model to predict three-dimensional free surface flow and breaking waves over complex topography. The model can act as a complementary approach to experimental investigations to gain further insight into the kinematics and dynamics of three-dimensional breaking waves. It is worth noting that the VOF scheme can also be switched off in the present model when considering the rigid lid approximation in open-channel flows [76–78]. Further work will include the surface tension effect for bubble generation and adaptive mesh [79], which can reduce computational efforts without sacrificing accuracy.

## Data accessibility

There are no data to report in this manuscript.

## Competing interests

We declare we have no competing interests.

## Funding

The author acknowledges the financial support from the Marie Curie EST fellowship for his PhD study at the University of Leeds. Z.X. is funded by the Marie Curie EST fellowship.

- Received February 14, 2015.
- Accepted June 18, 2015.

- © 2015 The Author(s)

Published by the Royal Society. All rights reserved.