Numerical simulation of wave transformation， breaking and runup by a contravariant fully non-linear Boussinesq equations model*

F. GALLERANO， G. CANNATA， F. LASAPONARA

Department of Civil， Constructional and Environmental Engineering， Sapienza University of Rome， via Eudossiana 18， 00184 Rome， Italy， E-mail:francesco.gallerano@uniroma1.it



Numerical simulation of wave transformation， breaking and runup by a contravariant fully non-linear Boussinesq equations model*

F. GALLERANO， G. CANNATA， F. LASAPONARA

Department of Civil， Constructional and Environmental Engineering， Sapienza University of Rome， via Eudossiana 18， 00184 Rome， Italy， E-mail:francesco.gallerano@uniroma1.it

In this paper we propose a new model based on a contravariant integral form of the fully non-linear Boussinesq equations（FNBE） in order to simulate wave transformation phenomena， wave breaking， runup and nearshore currents in computational domains representing the complex morphology of real coastal regions. The above-mentioned contravariant integral form， in which Christoffel symbols are absent， is characterized by the fact that the continuity equation does not include any dispersive term. The Boussinesq equation system is numerically solved by a hybrid finite volume-finite difference scheme. A high-order upwind weighted essentially non-oscillatory （WENO） finite volume scheme that involves an exact Riemann solver is implemented. The wave breaking is represented by discontinuities of the weak solution of the integral form of the non-linear shallow water equations （NSWE）. On the basis of the shock-capturing high order WENO scheme a new procedure， for the computation of the structure of the solution of a Riemann problem associated with a wet/dry front， is proposed in order to simulate the run up hydrodynamics in swash zone. The capacity of the proposed model to correctly represent wave propagation， wave breaking， run up and wave induced currents is verified against test cases present in literature. The results obtained are compared with experimental measures， analytical solutions or alternative numerical solutions. The proposed model is applied to a real case regarding the simulation of wave fields and nearshore currents in the coastal region opposite San Mauro Cilento （Italy）.

fully non-linear Boussinesq equations， contravariant formulation， Christoffel symbols， upwind WENO scheme， wet-dry fronts

Introduction

The modelling of surface wave transformation，wave breaking， shoreline movement and nearshore currents is of fundamental importance for the simulation of hydrodynamic phenomena which occur in coastal regions. Most of these phenomena can be represented by two-dimensional Boussinesq equations.

In general， coastal engineering issues are related to nearshore currents occurring in real situations which can be morphologically complex. In coastal areas，slightly sloping and regular sea beds alternate with steep irregular bottoms and the coastlines can be characterized by articulated shapes and be interrupted by the presence of anthropic structures and/or river mouths. In order to simulate hydrodynamic phenomena over computational domains characterized by a complex boundary， two strategies can be followed. The first strategy is represented by the possibility of using unstructured grids［1-5］. The second strategy is based on the numerical integration of the governing equations on a generalized curvilinear boundary conforming grid. By using curvilinear boundary conforming grids motion equations can be written in contravariant formulation［6-8］.

In recent literature integral forms［9］or new differential conservative forms of the Boussinesq equations expressed in terms of conserved variables［10-13］in which convective terms are expressed in divergence form and the term related to total local depth is expressed in a gradient form， have been proposed. The integration of integral forms or of differential conservative forms of Boussinesq equations allows the simulation of wave propagation from deep water regions up to the coastline， including the surf and swash zone. Breaking waves can be represented by discontinuities of the weak solution of the integral form of the nonlinear shallow water equations （NSWE）， numericallysolved by a shock capturing scheme［14-16］.

The above integral or differential conservative forms of the Boussinesq equations （which in the surf zone reduce to the NSWE by switching off dispersive terms） are able to simulate wave dynamics from deep water regions up to the coastline， do not need any additional term to take into account the wave breaking energy dissipation and do not require any empirical calibration. In order to apply a shock capturing method to Boussinesq equations expressed in differential conservative form， Roeber et al.［11］adopted a strategy consisting in using a hybrid finite volume-finite difference scheme.

In this paper a new model based on an integral form of the Fully Non-linear Boussinesq Equations in contravariant formulation， in which terms up to（withand， in whichis the wave amplitude） and second order vertical vorticity terms are included， is presented. In this contravariant integral form the continuity equation does not contain any dispersive term. The proposed Boussinesq equations are cast in a contravariant integral form in which Christoffel symbols are absent. The Boussinesq equation system is numerically solved by a hybrid finite volume-finite difference scheme: the convective terms and the terms related to the gradient of the square of the total local water depth are discretized by an high-order upwind WENO shock-capturing finite volume scheme， based on an exact Riemann solver， the dispersive terms and the terms related to the approximation to the second order of the vertical vorticity are discretized by a cell-centred finite difference scheme. The wave breaking is represented by discontinuities of the weak solution of the integral form of the NSWE. As a consequence no additional term to take into account the wave breaking energy is added in the equations. On the basis of the shockcapturing high order WENO scheme a new procedure，for the computation of the structure of the solution of a Riemann problem associated with a wet/dry front， is proposed in order to simulate the run up hydrodynamics in swash zone. The capacity of the proposed model to correctly represent wave propagation， wave breaking， run up and wave induced currents is verified against test cases presented in literature. The results obtained are compared with experimental measures，analytical solutions or alternative numerical solutions. The proposed model is used for the application on a real engineering case regarding the simulation of wave fields and nearshore currents in the coastal region opposite San Mauro Cilento （Italy）.

1. Contravariant form of the fully non-linear boussinesq equations

that represents the second order term in depth power expansion of the velocity vector， in whichis the two-dimensional differential operator defined byin a Cartesian reference system. Letbe the depth averaged value ofwhich is

The transformation relations of the genericvector from the Cartesian coordinate system to its covariant and contravariant componentsin the curvilinear coordinate system， as well as the covariant derivate and Christoffel symbols are given by:

Hereinafter the comma followed by a subscript position index indicates the derivative covariant operation.

In the curvilinear coordinate system， the Boussinesq equations expressed in contravariant formulation can be written as

In this paper we choose， as conserved variables，the total local depthand the contravariant quantity

a form of momentum equation， expressed in contravariant formulation， is given by

in which

the term （18） can be written in the form:

In order to obtain the contravariant integral form of the FNBE devoid of Christoffel symbols， the strategy devised by Gallerano and Cannata［15］is followed.

By integrating Eq.（13） and applying Green’s theorem to the second term on the left hand side， the continuity equation becomes

Equations （20） and （21） represent an integral form of the FNBE expressed in a contravariant formulation in which Christoffel symbols are absent. These equations are accurate up toin dispersive terms and retain the conservation of potential vorticity up to

2. Numerical scheme

The numerical integration of Eqs.（20） and （21） is carried out by a high order Upwind-WENO scheme，and for time integration of Eqs.（20） and （21） the strategy devised by Gallerano and Cannata［15］is followed.

In numerical integration of NSWE a particular treatment of the wet and dry front in the swash zone，to model a moving shoreline， is required［12］. In order to simulate the up rush and the backwash dynamics of the wet and dry front in the swash zone， the following original procedure is proposed.

For the sake of brevity the procedure of the up rush and the backwash dynamics of the wet and dry front is exposed referenced to a shoreline which is parallel to the curvilinear coordinate line. At the centre of the segments which separate the dry cellfrom the wet cell， point values of the unknown variables are reconstructed， by means of an asymmetric WENO reconstruction defined on the wet cell. For example， at the centre of the segment which is the interface between dry celland wet cell， WENO reconstructions defined on thecells lead to the evaluation of the variables. The advancing in time of the aforementioned variables is carried out by means of the exact solution of an apposite Riemann problem， with initial data given by the pair of point-values computed by the WENO reconstruction. It must be noted that the point values of the unknown variablesandare equal to zero because they belong to the dry cell

Generally speaking， the Riemann problem in a curvilinear coordinate system is more difficult to solvethan the Riemann problem for the same set of hyperbolic equations in an orthonormal frame. We solve all Riemann problems in a locally valid orthonormal basis.

For the sake of brevity， we define the normal and tangential to the coordinate linedepth-averaged velocity components asand， respectively. Asandare the unit vectors which are normal and tangential， respectively， to the coordinate lineand recalling relation （4）， the following transformation relations are obtained

For example， in the point of coordinatesbelonging to the segment that lies on coordinate line， which is the interface of cellsandthe WENO reconstructions lead to the definition of the point values of dependent variablesand. Letbe the propagation velocity of the wet and dry front at the advanced time levelLet，andbe the solution of the wet and dry Riemann problem defined by the hyperbolic homogenous system of the shallow water equations， written in the locally valid orthonormal basis. The exact solution of this Riemann problem on the interface between the wet celland the dry cell， gives

3. Test cases

The ability of the proposed model to accurately represent wave propagation， wave breaking， run up and wave induced currents is verified by means of test cases presented in literature. The tests are performed by using both Cartesian grids and highly distorted curvilinear grids. The results obtained are compared with experimental measurements， analytical solutions or alternative numerical solutions.

Fig.1 Run-up of a solitary wave on a conical island. Schematic plot of the computational domain

3.1 Solitary wave run-up on a conical island

In this section we simulate the run-up of a solitary wave onto a conical island. To this aim we numerically reproduce a laboratory test of Liu et al.［19］. A definition sketch for the computational domain used for the simulation is shown in Fig.1. The outer circle shows the base of the island， which is centred at（x，and has a radius， the middle dashed circle represents the initial still water shoreline（radius）， the inner circle represents the island top （radius）， the island height is 0.625 m.

Fig.2 Run-up of a solitary wave on a conical island. Maximum run-up around the conical island

Fig.3 Run-up of a solitary wave on a conical island. Snapshots of the computed free surface elevation at time t = 7.2 s（a）， t = 9.7 s （b）， t = 11.0 s （c）

As initial condition a rightward propagating solitary wave is imposed to the left boundary of the domain， on an otherwise calm free surface. The following expressions［20］are used for the free surface elevationand the depth averaged velocity componentin thedirection:

Fig.4 Run-up of a solitary wave on a conical island. Time series from the measurement locations depicted in Fig.1

In Fig.2 the comparison between the maximum computed run-up around the island and that measured by Liu et al.［19］is shown. The computed values of the maximum run-up around the island are in good agreement with the experimental data. It can be notedthat the run-up on the back side of the island， caused by the collision of edge waves circling the island from both sides， is well simulated by the present model. In Fig.3 snapshots of the computed free surface elevation at three different instants of the simulation are shown.

In Fig.4 the computed and measured time series from the measurement locations depicted in Fig.1 are shown. These four measurement locations have been chosen in order to represent the free surface elevation to the front， side， and rear of the island. From the comparison between the computed and measured values it can be seen that the proposed model is able to simulate the run-up at each measurement location around the island. Some secondary oscillation of the free surface elevation that has been observed during the laboratory experiments are slightly underestimated in the numerical simulation.

The good agreement between the computed and observed values of the free surface elevation shows the ability of the proposed model to adequately simulate the large run-up heights produced by a tsunami wave on the lee side of small islands.

3.2 Simulation of wave train propagation on a highly distorted grid

The FNBE proposed by Wei et al.［17］and Chen et al.［18］represent the evolution of the previous Boussinesq models because they improved the dispersive characteristics of the equations retaining terms up toin the variable depth power expansion and using as dependent variable of the equations the horizontal velocity from a certain depth. Furthermore the FNBE improve the non linear properties of the equation retaining terms up to， these terms represents the non linear effects in dispersion terms. Retaining these terms in the Boussinesq Equations leads to improvement in the model capability to represent the steepening wave profile just seaward of the surf zone. Then Wei et al.［17］demonstrate that FNBE extend a lot the validity envelope of the equations about dispersion and nonlinearities.

From a theoretical point of view， Wei et al.［17］underline that standard Boussinesq equations （SBE）have a very strict range of validity: the SBE are based on the assumption thatandafter which terms ofare neglected. Assuming that the validity limit for the standard approximation corresponds approximately to0.2， the value ofis around 0.04 （as shown in Fig.1 of Wei et al.［17］）， the fully non linear model extends the validity envelope of the Boussinesq model since the curve

In this paper a conservative and contravariant form of fully non lin near Boussinesq equations， in which terms up toand second order vertical vorticity terms are included， is presented.

In this section we demonstrate that the proposed conservative and contravariant form of FNBE， integrated with the proposed numerical scheme， is able to simulate waves with high nonlinearitiesalso on curvilinear highly distorted grids.

Fig.5 Simulation of wave train propagation on a highly distorted grid. Grid distortion

In order to demonstrate the model capability to correctly simulate the propagation of strongly dispersive and highly nonlinear wave trains on curvilinear highly distorted grid， we choose to simulate three wave trains with different characteristics by using a calculation grid with the highly distorted zone shown in Fig.5. Two sponge zones are laterally situated to correctly dissipate the wave energy. In order to analyze if the wave profile shows any deformation during the simulation， a flat bottom， situate at 0.85 m from the still water level， is used.

Fig.6 Wave profile. Wave train A:，，and

The first wave train （A） has an amplitude of， a period of， a non linearity ofand a strong dispersion of. In Fig.6 is shown the wave profile after 50 s of simulation. It is possible to notice how the proposed model can correctly represent the high dispersion of the wave train: the profile does not present any disturbs or modification of the wave characteristic and it is in good agreement with the data extrapolated from the analytical solution. In Fig.7 the free surface elevation of thiswave train is shown. From this figure it can be seen that the wave form is preserved and is not affected by alterations that could disturb the wave profile.

Fig.7 Instantaneous wave field. Wave train A:，and

Fig.8 Wave profile. Wave train B:，，and

The second wave train （B） has an amplitude of， a period of， a good non linearity ofand a good dispersion of0.25. In Fig.8 are shown the wave profile after 50 s of simulation. It can be seen （Fig.8） how the numerical results and analytical solution are in good agreement and the high distortion of the grid does not affect the model capability to represent non linear and dispersive waves.

Fig.9 Wave profile. Wave train C:，and

In Fig.9 is shown the wave profile related to train wave C with an amplitude of， a period of， a highly non linearity ofand a good dispersion of. It is shown how the wave profile， after 50 s of simulation， does not have any deformation and it is in good agreement with the analytical solution.

In conclusion， these tests demonstrate how the proposed model， based on a contravariant and conservative form of the FNBE that retains the term up toand second order vertical vorticity and integrated with a shock capturing high order WENO scheme， is able to represent strongly dispersive and highly non linear waves also on highly distorted calculation grid.

3.3 Simulation of wave field and currents in the coastal region opposite San Mauro Cilento （Italy）

In this section the proposed model is applied to the coastal region opposite San Mauro Cilento （Italy）. This coastal region is characterized by the presence of several inlets bounded by promontories and rocky protrusions that represent the boundaries of physiographic units. The proposed model is used to study the hydrodynamic fields in a coastal region opposite a portion of coastline， washed by the Tyrrhenian sea，that is bounded by two promontories that delimited a physiographic unit. We simulate the wave field and the nearshore currents produced by a wave train with wave heights of 3 m and a wave period of 10 s coming from 225oN. We are interested in the hydrodynamic fields produced by the interaction between the abovementioned incoming wave train and a system of emerged shore-parallel barriers.

Fig.10 Simulation of wave fields and nearshore currents in the coastal region opposite San Mauro Cilento （Italy）. Instantaneous wave field in presence of a system of shore parallel barriers

The proposed model takes into account the complexity of the hydrodynamic phenomena taking placein the above mentioned coastal region， consisting in non-linear wave-wave interactions between the incoming wave motion and the waves reflected and diffracted around the structures and coastal circulations induced by the wave breaking. The equations at the basis of this model are discretized on a boundary conforming curvilinear grid.

Fig.11 Simulation of wave fields and nearshore currents in the coastal region opposite San Mauro Cilento （Italy）. Details of the time averaged velocity field. A vector out of four is represented

In Fig.10（a） plane view of the simulated instantaneous wave field is shown. From this figure it is possible to see the reduction of the wave height caused by the wave breaking in the proximity of the coastline and the effects produced on the wave fronts by the emerged shore-parallel barriers. The figure shows the reflection of the incoming waves against the barriers and the bending of the wave fronts around the edge of the barriers. From the same figure， a drastic reduction on the wave heights landward of the barriers can be noted. In Figs.11（a）-11（c） the timeaveraged velocity field related to the wave field depicted in Fig.10 is shown.

From Fig.11（a） it is evident that part of the wave induced longshore current coming from the south is driven to the area between the barrier and the coastline，where its magnitude is not reduced. The aforementioned current flows in the alongshore direction up towhere it changes its direction producing a rip current， in the proximity of the coastline several recirculation cells are generated. The velocity field represented in Fig.11（b） shows also that this barrier drives a longshore current coming from the south in the area between the barrier and the coastline: this current flows in a rip current at. Figure 11（c） shows that the presence of the two adjacent barriers induce the formation of recirculating cells of opposite sign in the area between the barriers and the coastline: it is also evident the presence of a rip current flowing across the gap between the two barriers. The longshore current flowing from North to South in part is driven offshore， in part is intercepted by the barriers.

The analysis of the hydrodynamic field provides useful information for a qualitative assessment about the way in which the barriers， by changing substantially the wave field and the wave-induced velocity field，can modify the sediment transport conditions that are responsible of seabed modifications.

4. Conclusions

In this paper a new model based on a contravariant integral form of FNBE in order to simulate wave transformation phenomena， wave breaking， run up and nearshore currents in computational domains representing the complex morphology of real coastal regions has been presented. In the above-mentioned contravariant integral the Christoffel symbols are absent.

The wave breaking is represented by discontinuities of the weak solution of the integral form of the NSWE. A new high order shock capturing scheme has been proposed， in which an exact Riemann solver is involved. On the basis of the shock-capturing high order WENO scheme a new procedure， for the computation of the structure of the solution of a Riemannproblem associated to a wet/dry front， is proposed in order to simulate the run up hydrodynamics in swash zone.

It has been demonstrated that new model based on a contravariant integral form of FNBE proposed in this work， in which the Christoffel symbols are absent and which is numerically integrated by the new Upwind WENO scheme， allows the simulation of wave propagation from deep water regions up to the coast line including the surf and swash zone， and the simulation of longshore currents induced by the wave breaking on generalized curvilinear grids， even in the presence of highly distorted calculation cells.

It has been demonstrated that the presented Boussinesq model can be used for the simulation of wave fields and nearshore currents in the coastal region characterized by morphologically complex coastal lines and irregular seabed and by the presence of maritime infrastructures， such as in the real case dealt with in this paper， represented by the coastal region opposite San Mauro Cilento （Italy）.

References

［1］ SUN Zhong-bin， LIU Shu-xue and LI Jin-xuan. Numerical study of multidirectional focusing wave run-up on a vertical surface-piercing cylinder［J］. Journal of Hydrodynamics， 2012， 24（1）: 86-96.

［2］ WANG Ben-long， ZHU Yuan-qing and SONG Zhi-ping et al. Boussinesq-type modeling in surf zone using mesh-less least-square-based finite difference method［J］. Journal of Hydrodynamics， Ser. B， 2006， 18（3）: 89-92.

［3］ CASONATO M.， GALLERANO F. A finite difference self-adaptive mesh solution of flow in a sedimentation tank［J］. International Journal for Numerical Methods in Fluids， 1990， 10（6）: 697-711.

［4］ ZHANG Jing-xin， SUKHODOLOV A. N. and LIU Hua. Non-hydrostatic versus hydrostatic modelings of free surface flow［J］. Journal of Hydrodynamics， 2014， 26（4）: 512-522.

［5］ ZHOU W.， OUYANG J. and ZHANG L. et al. Development of new finite volume schemes on unstructured triangular grid for simulating the gas-liquid two phase flow［J］. International Journal for Numerical Methods in Fluids，2016， 81（1）: 45-67.

［6］ XU D.， DENG X. and CHEN Y. et al. On the freestream preservation of the finite volume method in curvilinear coordinaes ［J］. Computer and Fluids， 2016， 129: 20-32.

［7］ FANG Ke-zhao， ZOU Zhi-li and LIU Zhong-bo et al. Boussinesq modelling of nearshore waves under body fitted coordinate［J］. Journal of Hydrodynamics， 2012，24（2）: 235-243.

［8］ GALLERANO F.， CANNATA G. Compatibility of reservoir sediment flushing and river protection［J］. Journal of Hydraulic Engineering， ASCE， 2011， 137（10）: 1111-1125.

［9］ GALLERANO F.， CANNATA G. and VILLANI M. An integral contravariant formulation of the fully non-linear Boussinesq equations［J］. Coastal Engineering， 2014，83（1）: 119-136.

［10］ FANG Ke-zhao， ZHANG Zhe and ZOU Zhi-li. Modelling of 2-D extended Boussinesq equations using a hybrid numerical scheme［J］. Journal of Hydrodynamics， 2014，26（2）: 187-198.

［11］ ROEBER V.， CHEUNG K. F. and KOBAYASHI M. H. Shock-capturing Boussinesq-type model for nearshore wave processes［J］. Coastal Engineering， 2010， 57（4）: 407-423.

［12］ SHI F.， KIRBY J. T. and HARRIS J. C. et al. A high-order adaptive time-stepping TVD solver for Boussinesq modeling of breaking waves and coastal inundation［J］. Ocean Modelling， 2012， 43-44: 36-51.

［13］ FANG K.， ZOU Z. and DONG P. et al. An efficient shock capturing algorithm to the extended Boussinesq wave equations［J］. Applied Ocean Research， 2013， 43: 11-20.

［14］ DONG Jie， WANG Ben-long and LIU Hua. Run-up of non-breaking double solitary waves with equal wave heights on a plane beach［J］. Journal of Hydrodynamics，2014， 26（6）: 939-950.

［15］ GALLERANO F.， CANNATA G. Central WENO scheme for the integral form of contravariant shallow-water equations［J］. International Journal for Numerical Methods in Fluids， 2011， 67（8）: 939-959.

［16］ GALLERANO F.， CANNATA G. and TAMBURRINO M. Upwind WENO scheme for shallow water equations in contravariant formulation［J］. Computers and Fluids，2012， 62: 1-12.

［17］ WEI G.， KIRBY J. T. and GRILLI S. T. et al. A fully nonlinear Boussinesq model for surface waves. Part 1. Highly non linear unsteady waves［J］. Journal of Fluid Mechanics， 1995， 294: 71-92.

［18］ CHEN Q.， KIRBY J. T. and DALRYMPE R. A. et al. Boussinesq modeling of longshore currents［J］. Journal of Geophysical Research， 2003， 108（C11）: 2601-2618.

［19］ LIU P. L. F.， CHO Y. S. and BRIGGS M. J. et al. Runup of solitary waves on a circular island［J］. Journal of Fluid Mechanics， 1995， 302: 259-285.

［20］ NIKOLOS I. K.， DELIS A. I. An unstructured nodecentered finite volume scheme for shallow water flows with wet/dry fronts over complex topography［J］. Computer Methods in Applied Mechanics and Engineering，2009， 198（47）: 3723-2750.

March 7， 2015， Revised June 22， 2015）

* Biography: F. GALLERANO （1953-）， Male，

Ph. D.， Professor