Computation of free surface waves in coastal waters with SWASH on unstructured grids

Computation of free surface waves in coastal waters with SWASH on unstructured grids

Zijlema, Marcel

DOI

10.1016/j.compfluid.2020.104751

Publication date

2020

Document Version

Final published version

Published in

Computers and Fluids

Citation (APA)

Zijlema, M. (2020). Computation of free surface waves in coastal waters with SWASH on unstructured grids.

Computers and Fluids, 213, [104751]. https://doi.org/10.1016/j.compfluid.2020.104751

Important note

To cite this publication, please use the final published version (if applicable).

Please check the document version above.

Copyright

Other than for strictly personal use, it is not permitted to download, forward or distribute the text or part of it, without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license such as Creative Commons. Takedown policy

Please contact us and provide details if you believe this document breaches copyrights. We will remove access to the work immediately and investigate your claim.

This work is downloaded from Delft University of Technology.

ContentslistsavailableatScienceDirect

Computers

and

Fluids

journalhomepage:www.elsevier.com/locate/compfluid

Computation

of

free

surface

waves

in

coastal

waters

with

SWASH

on

unstructured

grids

Marcel

Zijlema

Delft University of Technology, Department of Hydraulic Engineering, P.O. Box 5048, 2600 GA Delft, The Netherlands

a

r

t

i

c

l

e

i

n

f

o

Article history:

Received 10 April 2020 Revised 30 July 2020 Accepted 30 September 2020 Available online 6 October 2020 Keywords:

Shallow water flow Nearshore waves Triangular mesh Finite difference Covolume discretization SWASH

a

b

s

t

r

a

c

t

Thispaperaimstopresenttheextensionofthenon-hydrostaticwave-flowmodelSWASHwiththe co-volumemethodtobuilddiscretizationschemesonunstructuredtriangulargrids.Centraltothismethod thatisfreeofspuriouspressure modes,istheuseofdualpairs ofmeshes thataremutually orthogo-nal,suchastheDelaunay-Voronoimeshsystems. Theapproximantssoughtarethecomponentsofthe flowvelocityvectornormaltothecellfacesoftheprimalmesh.Inadditiontothecovolumeapproach, anovelupwinddifferenceschemeforthehorizontaladvectiontermsinthemomentumequationis pro-posed.ThisschemeobeystheRankine-Hugoniotjumprelations andpreventstheodd-evendecoupling ofthevelocityfieldaccordingly.Moreover,caseswithflowdiscontinuities,suchassteadyboresand bro-kenwaves,areproperlytreated.Inspiteofthelow-orderaccuracyoftheproposedmethod,unstructured mesheseasilyallowforlocalrefinementinawaythatretainsthedesiredaccuracy.Theunstructured-grid versionofSWASHisapplicabletoawiderangeof2DHwave-flowproblemstoinvestigatethenonlinear dynamicsoffreesurfacewaves overvaryingbathymetries.Itsefficiency and robustnessistested ona numberoftheseproblemsemployingunstructuredtriangularmeshes.

© 2020TheAuthor.PublishedbyElsevierLtd. ThisisanopenaccessarticleundertheCCBYlicense(http://creativecommons.org/licenses/by/4.0/)

1. Introduction

A commonly encountered coastal engineering application

in-volvesthesimulationofdispersivewavesincoastalwaters.Waves

are usually presentalmost everywhere inthe domainofinterest,

while thewavedynamicscan bedominantina smallpartofthe

entiredomain,e.g.thesurfzoneduetowavebreakingand

gener-atingcurrents.Formanylarge-scalewaveproblems,themeshonly

needstobeconcentratedinparticularareaswithstrongwave

fea-turesandlargecurrentgradients. Localmeshrefinementis

there-foreasuitableapproachtooffersufficientgridresolutionin

other-wise under-resolvedareas(e.g.breakingzones,swashzones,

wet-lands,coastalflooddefences).

The developmentto increase the flexibility of structured grid

methods to allow variable grid resolution andlocal mesh

refine-menthasalonghistory.Traditionaloceanandcoastalmodels

usu-ally employ curvilinear grids with domain decomposition

tech-niques. Such gridsare typicallyorthogonal so thatthe coordinate

transformation offers advantagesin solution algorithm efficiency,

whereas undesirable numerical artefacts arising from grid

non-orthogonality,highcellaspectratioandskewnessareavoided.

E-mail address: m.zijlema@tudelft.nl

Structured grids greatly facilitate the use of high-order

dis-cretization methods. Such methods provide a higher spatial

ac-curacy andhencereduce thenumber ofgrid cellsrequiredfor a

givenlevelofsolutionaccuracy.Nevertheless,we arguethat

low-orderdiscretizationsonlocallyfinegrids mightbemoreefficient.

Low-orderupwind-biasedschemesoftenexhibitdesirablestability

properties butmay sufferfrom numerical damping,especially at

relatively coarse grids, which can easily be counteracted by grid

refinements. Furthermore, they only require small discretization

stencils and often deliver a good scaling capability for high

per-formance computing (HPC). In contrast, high-order discretization

methodsusuallylackrobustnessasthey havethepotential to

de-stroythestabilityoftheunderlyingnumericalapproach.In

partic-ular,centralandhigh-orderupwinddiscretizationstendtodisplay

non-physicaloscillationsandcommonlyresultininstabilities.

On the other hand, when local mesh refinement is applied

to flow problems in water of finite depth, the accuracy is

of-ten limited by the error in topographic and forcing data rather

than thediscretization error[1].Hence, grid refinement

substan-tiatesthebenefitoflow-orderdiscretizations.Despitethefactthat

low-order(upwind)methods aretypicallydeprecated inthe

liter-ature,we believe that such methods areconsidered adequate for

predictingwaveprocesses includingpropagation, shoaling,

refrac-tionandwavebreaking,providedthatsufficientlyfinemeshesare

https://doi.org/10.1016/j.compfluid.2020.104751

employed. While such processes are governed by mass and

mo-mentum balances, thismaynot be trueforthe transportof

con-taminants insurfacewatersand variabledensitysolute transport

(e.g.salinity,temperature,sediment).Thatistosay,highresolution

schemes aredesiredtoavoida detrimentallossoftracer

concen-trations[2].

Adrawbackofstructuredmeshesistheinherentdifficulties

as-sociated withthe generation ofappropriate boundary-fitted grids

and excessive grid refinement. The useof unstructured grids

of-fersaremedytotheseproblemsastheyaremoreflexibleand

eas-ilyprovide highlyvariablegridresolution incomplexgeometries.

Inthisregard,thedomainiscommonlydiscretizedintotriangular

ortetrahedralelementswithnoimpliedconnectivity.Unstructured

mesheshavereceivedgreatattentionforthepastseveralyearsin

the oceanandcoastalmodellingcommunity.Ageneralreviewon

thistopiccanbefoundin[3].Oneprinciplefeatureofan

unstruc-tured gridis thatadvanced gridrefinement techniquesare easier

to implement. Unstructured grid methods are therefore

particu-larly suitedforwave problems withcomplex geometriesthat

re-quireareasonablynon-uniformgridresolution.However, theydo

notallowforeaseofimplementationofhigh-orderdiscretizations.

Also,high-resolutionschemesonunstructuredgridsusuallydonot

takethefulladvantageofhigherorderaccuracythatcaneasilybe

achievedonstructuredgrids[2].Infact,low-ordermethodsfor

un-structuredmeshesarestillinusetodayinestablishedcodes[1].In

spiteofthislowerorderaccuracy,wepremisethatahigherrateof

convergence isachieved by localgrid refinementrather than use

ofhigh-order(central)discretizations.

Acommonfeatureofthetraditionalnumericalmodels

employ-ing unstructured meshes is the application of the finite volume

method,whichisa widelyusedapproachintheCFDcommunity.

Thismethoddirectlyutilizestheintegral formulationof

conserva-tion lawsin divergence form which are subsequently discretized

bymeansoffluxapproximationsappliedonnon-overlapping

con-trolvolumesineitherprimalordualgrids.Standardfinitevolume

schemesthushavethebenefittodiscretelyconserveunknownsof

PDEsonarbitrarymeshes.Inthecontextofthesolutionofshallow

waterproblems,suchschemestypicallyinvolvethevertex-centred

or cell-centred discretization of the conserved variables, i.e. the

waterdepthandthedepth-integratedvelocityvector,oncolocated

andsemi-staggeredgrids.

Yet,low-ordervertex-centredfinitevolumeschemesareknown

toberatherinaccurateoreveninconsistentonirregulargrids(for

details see, e.g.[4–6]). Such schemesrequirepolygon-shaped

du-alsthattypicallyprovokehighlydistortedvolumes.Forinstance,a

detailed analysisof Svärd etal. [4]reveal that vertex-centred

fi-nite volumeapproximationsyieldorderofaccuracylessthanone

fornon-smoothmixedgridsoftrianglesandquadrilaterals.In

par-ticular, thenon-uniformity in thesize ofdual mesh polygons

re-sultsinalossindiscretizationaccuracyofoneorder.Thecauseis

the differencesin fluxeson opposite sidesof thecontrol volume

that donot cancelin pairs when the mesh is non-uniform.It is

likely that their capacityto achieve desired physicalaccuracy on

irregulargridsisratherrestricted.Ontheother hand,thesolution

qualityofacell-centredfinitevolumescheme,operatingonprimal

cells, can easily be greater than that ofa vertex-centred scheme.

However,cell-centredschemesincurlargeroverheadsthan

vertex-centredschemesastheyinvolveapproximatelytwiceasmany

un-knowns.Therefore,wearguethatlow-orderfinitevolumemethods

usingunstructuredmeshesislesstractableforrobustandefficient

large-scalesimulationsofwavedynamicsrequiringhighaccuracy.

Another class of numerical methods involves the useof

stag-gered grids. In such grids the velocity components are staggered

withrespecttothepressure.Thepressureisplacedinthecentreof

thegridcells,whilethevelocitiesarestoredatthecellfaces.

Con-trarytocolocatedandsemi-staggeredgridschemes,staggeredgrid

methods prevent the well-known odd-even decoupling between

velocityandpressure whichistypicallythe keyingredient tothe

developmentof robust andefficient incompressible flow models.

Staggered finitedifference methods were first developedby

Har-lowandWelch[7]forincompressibleflowsonCartesiangrids1_and

hasalonghistoryofsuccessfulachievementsincoastalandocean

modelling(see,e.g.[9–11]).

Afruitfulextension oftheclassical staggeredgridapproach to

unstructured triangular and tetrahedral meshes is the covolume

method [12,13]. This methodis designed as a low-order method

that does not encounter spurious pressure modes. The covolume

method makes use of two dual orthogonal grids to approximate

the flow velocity and pressure.In this respect, the Delaunay

tri-angulation (primal grid) and the corresponding Voronoi diagram

(dual grid) arethe naturalchoice.Onlythe normalvelocity

com-ponents are defined on the faces of the primal mesh, whereas

the pressure is stored at the vertices of the dual grid, i.e. the

circumcentres of the primal cells. The accuracy of the covolume

scheme is generally first order in space but can be second

or-der accuratefora regular mesh[14].However, moreimportantly,

covolume techniques usually respect underlying physical

princi-ples dueto a range offundamental properties,such astopology,

conservationandsymmetry, leading tothe developmentof

high-fidelity discretization methods for incompressible Navier-Stokes

andMaxwell’sequations[15,16].The effectofdiscretizationerror

originatedfromsuchdiscretizationapproximationsisthuslimited,

in that the numericalsolution is merely influenced by the mesh

resolutionandmeshquality.

Over the past two decades, the covolume method, frequently

referred to as the unstructured C-grid discretization, has been

successfullyemployed forthe modellingof large-scaleocean and

coastalflows.See,forinstance,PerotandNallapati[17],Stuhneand

Peltier[18],KramerandStelling[19]andKleptsovaetal.[20].

Ad-ditionally,regionalcoastalandestuarineoceanmodelslikeUnTRIM

[21],SUNTANS[22]andD-FlowFM[23]arewidelyusedbymany

researchers.BothUnTRIMandSUNTANScanalsobeappliedto

in-cludenon-hydrostatic dynamicsforenvironmental flow problems

(e.g.internalwavesandoceanicfronts).

SWASH is a non-hydrostatic wave-flow model developed at

DelftUniversityandexploitsthefinitedifferencemethodon

stag-gered Cartesian grids [24]. This phase-resolving model has been

subjectedtoextensivebenchmarkingandhasbeensuccessfully

ap-pliedtoa numerouswave andflowapplications.Recentexamples

can be found in [25–28]. The present work describesthe

exten-sion ofSWASH to unstructuredtriangular grids. To thispurpose,

we adopt the covolume strategy because of the wish to

maxi-mizethephysicalaccuracy,stabilityandefficiencyofthenumerical

wavesimulation.Thefocusofthepresentstudyisontailor-made

discretizationssuitableforthesimulationofwave dynamics.They

ensurethecorrecthandlingofdiscontinuitiesandshocksas

mani-festedinbores,hydraulicjumpsandbrokenwaves.

In what follows, the shallow water equations including

non-hydrostaticpressureareintroducedinSection2.Inorderto

facil-itatetransparencyondiscretization schemes,thegoverning

equa-tions are presented in a depth-averaged form. Next, the concept

of the covolumeapproach and its numerical implementation are

treated in Section 3. Special attention is paid to the

discretiza-tionofmomentumadvectionbasedonagenuinelyupwind-biased

approachthat guaranteesconservation ofmomentumflux locally.

The key contribution of the present studyis to demonstrate the

efficiencyandrobustnessoftheproposedmethodforwave-related

simulations.InSection4anumberoftest casesare presentedfor

1 Historically, the use of a staggered grid for the solution of the shallow water

validationpurposesandtheresultsarediscussed.Conclusionsare

providedinSection5.

2. Governingequations

The governing two-dimensional, primitive variable equations

forthedepth-averaged,non-hydrostatic,freesurfaceflowofan

in-compressiblefluidoverabedtopographyd(x,y)aregivenby

∂ζ

∂

t +

∇

· q=0 (1)

∂

hu

∂

t +

∇

·

(

qu

)

+gh

∇

ζ

=−  _ζ −d

∇

pdz− cfu



u



(2)

withthecoordinatedirectionsx,yandzaligningintheeast,north,

andvertical directions,respectively. The bedleveld(x,y) is

mea-suredfromthereferencelevelz=0(positivedownwards),whereas

ζ

(x, y, t) is the surface elevation with respect to the reference

level (positive upwards). The water depth is givenby h

(

x,y,t

)

=

ζ

(

x,y,t

)

+d

(

x,y

)

.Furthermore,u=

(

u,

v

)

istheflowvelocitywith

the depth-averaged components u(x,y, t) andv(x,y, t) along the

x and y coordinates, respectively, q=huis the mass flux, p(x,y,

z,t)isthenon-hydrostaticpressure(normalisedbythedensity),g

isthegravitationalacceleration,andc_fisthedimensionlessbottom

frictioncoefficient.Theshallowwaterequations(1)and(2)are

de-rived fromintegratingthemassconservationandthemomentum

balanceoverthedepth,respectively,whereas thetotalpressureis

decomposed intoits hydrostatic andnon-hydrostaticcomponents

(see, e.g. [21,22,24]). In the case of a hydrostatic pressure

distri-bution, i.e.p ≡ 0, these equations can be reformulatedas a set

ofnonlinear hyperbolicequations,andmaythus generate

discon-tinuous solutions featuring shockwaves [29]. Such solutions can

readilybeunderstoodasweaksolutionsinthevariationalcontext.

Using Leibniz’ rule,a conservative expression forthe gradient

ofnon-hydrostaticpressureisobtained[30]

 _ζ

−d

∇

pdz=

1

2

∇

(

hpb

)

− pb

∇

d

with pb the non-hydrostaticpressure atthe bed.This pressureis

associatedwiththeverticalmotionthatisgovernedbythe

follow-ingequation

∂

ws

∂

t = 2pb h −

∂

wb

∂

t

wherewsisthevelocityinthez−directionatthefreesurface.This

equationisderivedusingtheKeller-boxmethodtofurtherimprove

thedispersivebehaviourofthewaves[24,30].Theverticalvelocity

atthe bedwb canbe foundby meansofthefollowing kinematic

condition

wb=−u·

∇

d atz=−d

Finally, the system of equations is complete with the following

equation

∇

· u +ws− wb

h =0 (3)

whichensuresconservationoflocalmass.

The momentum equation (2)is symbolically written invector

formwhichisnotthemostnaturalformtodescribethe

underly-ing physics in wave-dominatedregimes.Instead, we considerthe

followingmomentumbalancecharacterizingtheflowbelowafree

surface in water that moves in a nominaldownstream direction

definedbytheunitvectornf,

∂

huf

∂

t +[

∇

·

(

qu

)

]· nf+ghnf·

∇ζ

=−1

2nf·

∇

(

hpb

)

+pbnf·

∇

d− cfuf



u



(4)

withu_f=u· nf thelocal streamwise,depth-averaged flow

veloc-ity. Vector nf is assumed to be slowly varied such that pressure

gradientsredirectmomentumthroughbendswhereastheflow

ve-locityu_f isalignedwiththe directionofthepressure forces.

Fur-thermore, the lateral variability in the free surface slope is

neg-ligible on the scale offlow induced motion,especially in coastal

regions[31].

The depth-integrated momentum equation (4) essentially

de-scribesalocalbalancebetweenflowacceleration,frictional

decel-erationandthepressureforceactingonanelementoffluidalong

streamlines,andgovernsmomentum perunit volume normalised

bythedensity,i.e.uf.Themassfluxq isthetransportvelocityof

momentumandthesecond termofEq.(4)representsthe

advec-tive acceleration of the flow field witha nonlinear spatialeffect

onuf.Forinstance,itcauses thechangeinthe momentumin an

open channel dueto expansion orcontraction, i.e.slowing down

orspeedingupthefluidlocally.Suchflowtransitionsaretypically

rapid and frequently arise in hydraulic jumps, dam breaks, tidal

boresandflowsoveraweir.

3. Discretizationonunstructuredtriangularmeshes

This paperis not intended to provide an outline of the

over-all numerical methodology of SWASH, but rather to explain the

essential steps of the discretization of the governing equations

on unstructured staggered grids since thistopic differs from the

structured-grid version of SWASH as described in Zijlema et al.

[24].Therefore,thediscussionontemporaldiscretizationsandthe

solution procedure is not included in the present study asthey

havebeenpreviouslytreatedindepth(in,e.g.[24]).

The starting point is that the flow velocities are segregated

from the surfaceelevation and non-hydrostatic pressure and are

assignedto differentgrid locationsbased ona staggered grid

ar-rangement.Thisisfurther discussedinSection3.1.Notealsothat

thesolutionofnon-hydrostaticpressureisnotpursuedherein.The

interestedreader is referred to the aforementioned referencefor

detaileddescriptionofthepressurecorrectiontechniqueaimingto

enforce localmass conservation.Section 3.2 outlines the concept

ofthecovolumeschemeforspatialdiscretization.The

implementa-tionofthismethodinSWASHiselaboratedindetailinSection3.3.

Next,thediscretizationoftheadvectiveaccelerationisdealtwith

inSection 3.4.Finally,a simplewet-dry schemeistreatedbriefly inSection3.5.

3.1. Staggeredgrid

Staggered grid methods typically employ velocity components

astheprimary discreteunknownsin thespatialdiscretization.In

thisrespect,theflowvelocity componentufisa properone.This

physical variableis located at the face midpoint inthe direction

nfnormaltothecellface,whereasthephysicalscalars(e.g.water

depth,pressure) are positioned at thecell circumcentre.As a

re-sult, theflow velocity normalthroughthe faceofthe gridcell is

influenced directlyby the waterdepth andnon-hydrostatic

pres-sureatthecircumcentersoneitherside,becauseitisthepressure

gradient that drives the flow. In addition, the change in volume

cell, confinedby thesurface elevation, is dueto the massfluxes

hu_fat thefaces ofthe gridcell. Thisphysicaldualityexpresses a

naturalpropertyofthestaggeredmeshapproach.

Staggered schemes thus require the existence of the cell

cir-cumcentressincethelinesegmentconnectingthecircumcentresof

thetwo cellsintersects withtheirshared face,whilethey are

or-thogonaltoeachother.Themostcommonpolygonswithaunique

circumcenter are the triangles and rectangles. Furthermore, any

twonon-colinearsidesofsuchcellsspanabasisinR2_.(Notethat

sides.)Thisimpliesthatthevelocityvectoru(momentumperunit

volumenormalised bythedensity)atanarbitrarylocationwithin

the cell can be described by a linear combinationof the normal

components u_f at the cell faces; see also next section. However,

inthecaseofrectangles,bothindependentvelocitycomponents,u

andv,arerequiredtocharacterizemomentuminthetwo

dimen-sionalspaceduetotheirmutualorthogonality.

3.2. Thecovolumemethod

The mainmeritoftheabovestaggered treatmentofthe

prim-itive variables is the proven robustness and efficiency as these

variables are located in optimalpositions on the mesh withtwo

unique consequences. First, it prevents the odd-even decoupling

betweenthepressureandthenormalvelocitycomponent.Second,

it yields excellent conservationproperties.Forexample,both

pri-mary (mass and momentum) and secondary (e.g. kinetic energy

and enstrophy) quantities are conserved by Cartesian staggered

grid methodsfortheincompressibleNavier-Stokesequations(see,

e.g.[14,32,33]).Theconservativebehaviourofstaggeredschemesis

dueto thediscrete divergenceandgradient operators thatmimic

thephysicalpropertiesoftheircontinuouscounterparts.Asa

con-sequence, thephysicalaccuracyis preservedatthediscrete level.

In this context, the staggered discretization is classified as the

mimeticapproach[15].

The covolumemethod[12,13]isoneofthemanymimetic

dis-cretizationsthathaveappearedinliterature,ofwhichanexcellent

overview is given by Lipnikov et al. [15]. The mimetic character

ofdiscretecovolumeoperatorsprovidesconservationpropertiesin

theresultingnumericalmethods.Forinstance,thelocalandglobal

conservationpropertiesofthecovolumeschemeappliedtothe

in-compressibleNavier-Stokesequationshavebeenextensively

exam-inedbyPerot[14]andZhangetal.[34].

The covolumemethodusuallyrequirestwoorthogonal meshes

to compute the change in mass and momentum ofshallow

wa-ter flows. The Delaunaymesh and its dual, the Voronoi

tessella-tion,aretheobviouschoice.TheverticesoftheVoronoimeshare

the circumcentres ofthe Delaunay cells. Like the Cartesian

stag-geredapproach,onlythenormalvectorcomponentisemployedat

the face midpoint of the primal grid and the pressure is

conse-quently defined atthe cell circumcentres,i.e. the vertices ofthe

dualgrid.Incontrast,forinstance,themethodimplementedin

EL-CIRC [35] is not based on the covolume approach as it requires

thesolutionofthevelocity componenttangentialtothecellface.

Suchasemi-staggeredschemecangiverisetoerroneous

pressure-velocity oscillations[36].Ingeneral,preservationofthekey

phys-ical variables, i.e. thepressure at the cell centreand the normal

velocitycomponentu_fatthecellfaceratherthanthefull

momen-tum vector u, in distinct discrete structures such as primal and

dualmeshesiscrucialtopreventnon-physicalartefacts.(SeeTonti

[16] for details on the central role of the relationship between

physical variables andmesh topologies, i.e.points, lines,surfaces

and volumes,in computational physics.)Inaddition, vector fields

canbereconstructedbymeansofaninterpolationmethod.

Exam-plesrelatedtotriangularmeshescanbefoundin[37–40].

It isimportanttonote thattheuseofcovolumeorC-grid

dis-cretizationsmightbeproblematicforthesimulationofgeophysical

flows[41].Forinstance,spuriousgeostrophicmodesmayoccuron

both structured andunstructured meshes in the presence ofthe

Coriolis forceandvarying bathymetry. Aproper reconstruction of

the tangential velocity is the key solution to this problem [39].

Furthermore, three-dimensional hydrostatic models for resolving

large-scalebarocliniccurrentstendtoexhibitcheckerboardpattern

inthedivergenceofthehorizontalvelocityfield,andhenceinthe

vertical velocity(cf.Eq.(3)),duetothecombined useof

triangu-largeometry,notablyequilateraltriangles,andvariablestaggering

[41,42]. Some remedieshave beenput forwardin alleviatingthis

problem such as divergenceaveraging, increaseddiffusion,

filter-ingandmimeticdiscretizations(see,e.g.[42–44]),thoughthe

un-wanted modesareusually not entirelyabsent.Nevertheless,such

checkerboardmodeshavenotbeenencountered inthenumerical

wavesimulations aspresentedinthisstudy.Therefore,nofurther

attentionwillbepaidtothisissueinthispaper.

Althoughthecovolumemethodisknowntobefirstorder

accu-rateforgeneralunstructuredmeshesandsecondorderforaclass

ofregulargrids2_{,ityieldsphysicallyconsistentsolutiontothe}

gov-erning PDEs. Furthermore, discretizations should reflect accurate

changesintheresponseofthewaterwavemotiontoadvective

ac-celerationandpressuregradientand,inturn,bathymetricforcing.

From this perspective, we postulate that low-order mimetic

dis-cretizationsareaviablealternativetohigh-orderfinite-volume

dis-cretizationschemes, includingcolocatedandsemi-staggered ones

employing the full momentum vector as the primary unknown.

Despite their higher order accuracy, numerical methods that do

not meet certain physical and mathematical properties can

pro-duce unacceptably large errors due, principally, to the presence

ofnonphysicalparasiticmodesandthe nonlinearamplificationof

such modes and discretization errors [32]. The above reasoning

motivatestheimplementationofthecovolumeschemeinSWASH

whichispresentedinthenextsection.

3.3. Covolumediscretization

Weconsideratwo-dimensionalgridconsistingofDelaunay

tri-angles.Notethatthefacesinthedualmeshareorthogonaltothe

cell faces in the primal mesh. Letindices c andf enumerate

tri-angularcellsandfaces,respectively.Inthediscretizationapproach

we denoteSc thesetoffaces forming theboundaryofcell cand

vectorn_c,f theoutward-pointing normaltofacefofcell c.Rather

thanemployingtheoutwardnormal,theunitnormaltofacef

co-incidedwiththestreamwisedirectionnfisadopted.Sincethis

nor-malcaneitherbeoutwardorinwardwithrespecttocellc,we

se-lectoneofits directiontobefixedatallfaces inthegrid,sothat

uf· nf=uf withufthevelocityvectoratfacef.Themutual

orien-tationoftheoutwardnormalnc,fandtheuniquenormalnfatface

fofcellcisgivenby(cf.Fig.1) nc, f =

α

c, fnf,

α

c, f =nc, f· nf=±1

The continuity equation (1)is integrated over a cell c as

fol-lows  c



∂ζ

∂

t +

∇

· q



dA=0

Application of the midpoint rule and the divergence theorem of

Gaussyields Ac d

ζ

c dt +  ∂c q· ncdl=0

with

ζ

cthewaterlevelstoredatthecellcircumcentre,Acthearea

ofcell candnc the outwardunit normalvector tocell boundary

∂

c.Theintegralisapproximatedbytakingthesumoverthefaces

ofthecell  ∂c q· ncdl≈  f∈Sc qf· nf



nc, f· nf



lf=  f∈Sc

α

c, fQf

whereQ_f=q_f· nflf isthe mass flux integratedalong cell facef

withlf thecorrespondinglength. Notethat

α

c,fQf isnegativeifqf

2 Formally speaking, the discretization error for the pressure gradient term is sec-

ond order in space if the circumcentres of two adjacent cells are equally spaced from the shared face.

Fig. 1. The layout of the triangular mesh for the discretization of the momentum equation at face f for u f > 0. This face with length l f is shared by two triangle cells

c L and c R with distance s f between the associated circumcentres. Also depicted

are the unit outward-pointing normal n c_L ,k to face k of cell c L and the unit normal

vector n k indicating the flow direction. The normals are constant along each face.

Definitions of various geometric quantities and variables are provided in the text.

directsintocellc,otherwise,itispositive.(SeeSection3.5for

fur-therdetailsonthemassflux.)Theresultingsemi-discretescheme

forEq.(1)is Ac d

ζ

c dt +  f∈Sc

α

c, fQf=0 (5)

Next,themomentumequation(4)issolvedatthefacesofeach

interiorcell.Theapproachtofollowisessentiallyafinitedifference

methodasthespatialdiscretizationwillbecarriedoutdirectlyfor

the equationinsteadofits integral formexpressingthe

conserva-tionbalanceofmomentum.

Fig.1showscellfacefanditstwoadjacentcellscLandcR,with

the subscriptsindicating theposition ofthe cellswithrespectto

theface.

Thedistancebetweentheassociatedcircumcentresisgivenby



sf=



sc

L, f+



scR, f

with



sc,f the distance fromthe circumcentreof cell c to face f.

Wecomputetherateofchangeofthestreamwisemomentumper

unit facelengthbetweentwocircumcentresadjacenttofacef,i.e.



sfhfuf,withufthevelocitycomponentnormaltofacefandhfthe

waterdepth atface f.This isdetermined by,amongstothers,the

pressure forcesacting onthe area



s_f (perunit face length)and

the advective flow acceleration altering spatially the momentum

intheupstreampartofthearea,i.e.



sc,fhfuf(perunitfacelength),

eitherwithc₌c_L(ifu_f>0)orc₌c_R (ifu_f<0).Thisupwindbias

ensuresthestabilityofthepresentmethodsinceitinheritscertain

physicalconditions,aswewillseeinthenextsection.

Byvirtueofthemeshorthogonalityproperty,weobtainthe

fol-lowingfinitedifferenceform



sf d˜hfuf dt +



sc, faf+ghf

ζ

c_R−

ζ

c_L

=−1 2

hc_Rpc_R− hc_Lpc_L

+p˜f

dc_R− dc_L

− cf



sfuf



u˜f



(6)

whereafisthecomponentofmomentumadvectioninthe

stream-wisedirection.ThiswillbetreatedinSection3.4.Furthermore,pc,

hc anddc denotethepressure,thewaterdepth andthebedlevel

locatedatthecircumcentreofcellc,respectively.Finally,h˜_f,h_f,p˜_f

andu˜faretheinterpolatedwaterdepths,non-hydrostaticpressure

andvelocityvectoratfacef,respectively;seebelow.Notethat

sub-scriptbhasbeendropped fromthepressurevariables.Itisworth

noting that, owing to the mesh orthogonality, the discrete

pres-suregradients ofEq.(6)and thediscrete divergenceofthe mass

flux, Eq.(5), are mutuallyadjoint.This dualitycontributesto the

robustnessofthecovolumemethod(see,e.g.[15]fordetails).

Wenowprovideappropriateinterpolationsforhf,pfanduf,

re-spectively.Letusconsider aflow acrossan abruptbedchangeat

cellfacef.Aconditionthatmustbemetisthehydrostaticbalance

betweenthe(hydrostatic)pressurefluxandthetopographyterm

1 2g[h]

2

f=ghf[d]f

wheretheoperator[· ]f denotesthe jumpacrossface f.This

con-ditionpreservesthequiescentwaterovervarying bathymetry,i.e.

u=0and

ζ

=constant.Hence,Eq.(6)reducesto

ghf

ζ

c_R−

ζ

c_L

=0 Bydefining hf = 1 2

hc L+hcR

wehave ghf

ζ

c_R−

ζ

c_L

= 1 2g

hc_L+hc_R

hc_R− hc_L+dc_L− dc_R

=0 (7) sothat 1 2g

h2 c R− h 2 c L

=ghf

dc_R− dc_L

whichisthehydrostaticbalanceatthediscretelevel.

Next,h˜f in thetime-derivative term ofEq. (6)represents the

volumeperunitareaatfacef.Hence,forconsistency,wedefine



sfh˜f=



sc_L, fhc_L+



sc_R, fhc_R

Likewise,thenon-hydrostaticpressureatfacefisapproximatedby

takingthe volume-weighted average of the pressures of the two

cellsadjacenttotheface

˜ pf=



sc_L, f



sf pc_L+



sc_R, f



sf pc_R

Since the covolume scheme employs the face normal

veloc-ity components uf, a velocity vector u=

(

u,

v

)

needs to be

re-constructedfromthesecomponents.Onereconstructionmethodis

dueto Perot[14]and computesthe cell-basedvelocity vector uc

bymeansofthedivergencetheorem.Theresultisgivenby

uc= 1 Ac  f∈Sc lf



sc, fufnf (8)

which expresses the cell velocity vector in terms of the normal

velocity components at the faces of the cell. Perot [14] argued

that computing thecell-based velocity vector by means of

inter-polation(8)impliesdiscrete conservationofmomentum(see also

AppendixA). Thisinterpolationisfirst orderaccurate andsecond

orderifthegridisregular.Analternativeisthemethodproposed

byVidovic[38]usingaleastsquarestechniquewhichprovides

bet-teraccuracy,butisratherinvolvedandcomputationallydemanded.

We therefore prefer the use of Perot’s interpolation scheme(see

thevolume-weightedaverageofthecell-basedvectorsofthetwo

cellsadjacenttotheface,asfollows

˜ uf=



sc_L, f



sf uc L+



sc_R, f



sf uc R

3.4. Discretizationofmomentumadvection

In [45],the incorporationof the Rankine-Hugoniotjump

rela-tions in anumericalscheme isdevised topreclude theodd-even

decoupling in the velocity field. In addition, fulfillment of these

jumpconditionsisdesiredtoaccuratelycaptureflowswith

discon-tinuities.Inthisregard,theadvectiveflowaccelerationisprojected

onthedirectionnfatfacefrepresentingalocalshock,

af=ac· nf

withac the upstream cell-basedadvectionvector. The manner in

which this nonvolumetric vector is to be determined is closely

linked to the fulfillment of the jump conditions atthe

disconti-nuity. Morespecifically, considersteady, frictionlessflow acrossa

jumpatfacef,undertheassumptionofhydrostaticpressure,then

Eq.(6)becomes



sc, faf+ghf

ζ

c R−

ζ

cL

=0 (9)

and ac is constructed such that Eq.(9) reduces to the following

Rankine-Hugoniotrelation



qu+1 2gh 2



f= ghf[d]f (10)

withqthemassfluxperunitwidth.Thiswillimprovephysical

ac-curacy atthediscretelevel sincethebalance betweenthe

advec-tive accelerationandthepressuregradient isrespected.

Addition-ally, such a schemeremains accurate on coarse meshes, whereas

any consistent scheme that does not adhere to jump conditions

explicitly converges relatively slower to a weak solution.

Exam-pleswheresuchalatterschemehasbeenemployedare[14,22,23].

Thus,fulfillingtheRankine-Hugoniotcondition(10)hasthe

poten-tialtofurtherenhancetherobustnessandefficiencyofthe

numer-icalmethod[45].

To achieve the desired jump conditionat face f, we use

one-sideddiscretizationsformomentumadvection[45].First,we

con-sider acellcupwindofface fwherebyvelocityuf ispointingout

ofthecell.LetchoosethiscellcL,i.e.uf>0(cf.Fig.1).Vectorac_L

is thecell-basedadvectionvector oftheupwindcellandis

com-putedusingGauss’theorem

acL = 1 Ac_L cL

∇

·

(

qu

)

dA= 1 Ac_L  ∂cLq· ncLudl ≈ 1 AcL  k∈ScLq· nk



ncL,k· nk



ˆ uklk = 1 AcL  k∈ScL

α

cL,kQkuˆk

The transportedflow velocityvectorû_katthefaces kofcellcL is

obtainedbymeansofaone-sidedinterpolation(see,e.g.[19,20])

ˆ uk≈ ucL,k with uc

L,k the cell-basedvelocity vector of theupwind cellcL of

face k (see Fig. 1; keep in mind that uf > 0). In turn, this

cell-basedvelocityvectorisinterpolatedfromfacenormalcomponents

usingEq.(8).Finally,wearriveatthefollowingexpressionforthe

facenormalcomponentofmomentumadvection

af= 1 Ac L  k∈Sc L

α

c L,kQkucL,k· nf (11)

Approximationinthecaseofnegativeflow,i.e.uf<0,isdone

sim-ilarly. Thiscompletes ourdiscretization ofthemomentum

advec-tion.

We now demonstrate that jump condition(10) is indeed

ful-filled.SubstitutingEq.(11)inEq.(9)andusingEq.(7)gives



sc L, f Ac_L  k∈Sc L

α

c_L,kQkuc_L ,k· nf+ 1 2g

h2 c R− h 2 c L

=ghf

dc_R− dc_L

Letk1,k2 andfdenotethefacesoftheupwindcellcL,seeFig.1.

Expansionresultsin −



scL, f Ac_L Qk1uc_L ,k1 · nf−



sc_L, f Ac_L Qk2uc_L ,k2 · nf +



scL, f Ac_L Qfuf+ 1 2g

h2 c R− h 2 c L

=ghf

dc_R− dc_L

sinceuc

L, f· nf≈ uf· nf=uf.Furthermore,letdesignate

qk=



sc_L, f

Ac_L

Qk

therateofvolumeflowper unitwidthoftheupwindcell.Hence,

wehave qfuf−

qk1uc_L ,k1 +qk2uc_L ,k2

· nf+ 1 2g

h2c_R− h2c_L

=ghf

dc R− dcL

which is the Rankine-Hugoniot jump condition at the discrete

level.Notethatthisconditionholdsonfacefrepresentingthe

dis-continuity,whereasthefirsttwotermsexpressthejumpupstream

ofthatface.

Forreasonsofefficiency,Eq.(6)isnotsolved asitstands,but

the solutionof an equation forflow velocity u_f is considered

in-stead. For the purpose of deriving this equation, the right hand

side of Eq.(6) is set to zero. Furthermore,we considerthe case

ofpositiveflowatfacef,i.e.u_f>0.Weobtain



sfh˜f duf dt +



scL, f dhc_L dt uf+



scL, faf+ghf

ζ

c R−

ζ

cL

=0

astheoutflowatfacefonlyaltersthewaterdepthupstream.Next,

substitutingEq.(5)yields



sfh˜f duf dt −



scL, f  k∈Sc L

α

cL,kQk Ac L uf+



sc L, faf +ghf

ζ

c R−

ζ

cL

=0

andsubsequently, thecomponentofmomentumadvectionaf.We

thenhave



sfh˜f duf dt +



scL, f  k∈Sc L

α

cL,kQk

uc L,k· nf− uf

Ac_L +ghf

ζ

c_R−

ζ

c_L

=0

The outgoingflux atface fcan be omitted withoutchanging the

amount of momentum because uc

L, f· nf=uf. Hence, including

theright handside ofEq.(6),the final semi-discretized

duf dt +



sc_L, f



sfh˜fAc_L



Qk1

uf− uc L,k1 · nf

+Qk2

uf− uc L,k2 · nf



+ghf ˜ hf

ζ

c R−

ζ

cL



sf =− 1 2˜hf hc RpcR− hcLpcL



sf +p˜f ˜ hf dc R− dcL



sf −cf uf



u˜f



˜ hf (12)

Asimilar semi-discreteequation canbe derivedforthecasewith

negativeflow(u_f< 0).Notethatthecurrentdiscretizationdiffers

from the one proposed by Kramer and Stelling [19] in that it is

abletoresolveflowdiscontinuitiesmoreaccurately.Thisis

demon-stratedinSection4.2.

Anattractivecharacteristicoftheresultingfinitedifference

ap-proximation,apartfrompreservingphysicalaccuracy,isthatit

en-sures total momentum is conserved at the discrete level. Such a

mathematicalpropertyisdesiredtokeeptheoveralldiscretization

errorlow byeliminating theconservationerrors.Aproof isgiven

inAppendixA.

3.5. Wettinganddrying

ThemassfluxintegratedalongacellfaceinEq.(5)isevaluated

with

Qf =qf· nflf =lfhˆfuf

wherehˆf isthewaterdepthapproximatedatfacef.Sincethe

wa-terdepthisstoredatthecircumcentresanupwind-biased

interpo-lationisapplied,asfollows

ˆ hf=

⎧

⎪

⎪

⎪

⎪

⎪

⎪

⎨

⎪

⎪

⎪

⎪

⎪

⎪

⎩

ζ

c L+min

dc L,dcR

, ifuf>0

ζ

c R+min

dc L,dcR

, ifuf<0 max

ζ

c L,

ζ

cR

+min

dc L,dcR

, ifuf=0 (13)

so that thewater depth inthe outgoing mass flux ofany cell is

that of the cell itself.As a consequence,a sufficient conditionto

guaranteeanon-negativewaterdepthincellccanbederived.

De-tails can be found in Kramer andStelling [19]. This condition is

givenby



tlf

|

uf

|

≤ Ac

andsoitforcesatmostonecellpertimesteptobeflooded.

If a wetcell becomesdry, thewater depthin thecell can be

arbitrarysmall,whilethereisnooutgoingmassflux.Theflow

ve-locity atface fis thenset to zeroifhˆf<

δ

with

δ

the threshold

depth(

δ

₌10−10minthepresentstudy),whilemomentum

equa-tion(12)isnotsolved.Thiscompletesourdescriptionofthe

wet-dryalgorithm.

4. Testcases

4.1. Introduction

This section presents results for four test cases with an

in-creasing degree ofdifficulty inthe unstructured mesh

configura-tion. With the exception of the first case, the selection of these

cases was motivated by the wish to investigate the main

fea-tures ofthe nearshorewave dynamics, suchas wave

transforma-tion due to shoaling, refraction, diffraction and runup, and

non-linearprocesses,inparticulartheresonant three-waveinteraction

and depth-induced breaking. The objective is to assess the

abil-ityofthepresentunstructuredstaggeredmeshapproachto

repro-duce thesewave processeswithan emphasis onpredictive

accu-racy. The predictive performance of SWASH for simulating these

wave-relatedcasesonrectangulargridsiswelldocumentedin,e.g.

[24,46,47]. The first academic test caseis added to thisstudy to

verifythe numericalaccuracywith whichthe momentum

advec-tionisapproximated.

Some brief remarks are given here in relation to marching

in time and imposing boundary conditions. The semi-discrete

Eqs. (5) and (12) are integrated in time using a second order

explicit leapfrog scheme in conjunction with the explicit Euler

methodfortheadvectiontermandtheimplicitEulermethodfor

the non-hydrostatic pressure term[24]. For stability reasons the

time stepisrestricted tocomplywiththeCFL conditionthat

de-pendsonthewavecelerity.Thisconditionisgivenby

Cf=



t



gh˜f+

|

uf

|





sf ≤ 1

with



t the time step andCf the Courant number evaluated at

a face midpoint. A variable time step is employed so that the

Courantnumberislessthanaprescribednumber(<1)inallwet

faces. The balance of experience derived from ourvarious

wave-relatedstudies suggeststhatthisroute hasemerged asproviding

thebestcompromiseintermsofaccuracyandstability.Forfurther

detailsonthetimeintegrationseeZijlemaetal.[24].

Furthermore, boundary conditions considered here pertain

mainlyto the generationofincident free shortwaves andbound

infragravity waves. Theyare imposed as weakly reflective at the

offshore boundary.Details ontheir implementation can be found

in[24,47].

As a final note, thestaggered meshmethod described inthis

paperwill soon be released in a futureversion of SWASH (http:

//swash.sourceforge.net),thoughthesourcecodeandadetailed

de-scriptionoftheimplementation canbe obtainedfromtheauthor

uponrequest.

4.2. Dambreakoverwetbed

The aimis to demonstratethe accuracy withwhich flow

dis-continuities are reproduced using the proposed upwind scheme

outlined in Section 3.4 and the scheme of Kramer and Stelling

[19]forthemomentumadvection.Weconsider anunsteady flow

withpropagationofadambreakwaveinahorizontalfrictionless

channel.Thepressureisassumedtobehydrostatic.Thelengthand

thewidthofthebasinwas100mand10m,respectively,whereas

thedamwasinitiallylocatedinthecentreofthebasin.Theinitial

upstreamwaterdepthwas1m,andthedownstreamwaterdepth

wasset to0.1 m.The velocity wasinitially zero. The rather

uni-formtriangularmeshemployedhasanaveragedgridsizeof0.4m,

whereasthetimestepwassetto0.01s.

Attimet=7s afterdam failure,themodelresultsofthetwo

schemesare depictedinFig.2.Alsoshownis theanalytical

solu-tion.

The results obtained with the proposed scheme are in better

agreementwiththeanalyticalsolutionthantheschemeofKramer

and Stelling [19]. In particular, the simulation showsthat a

cor-rect speed and height of thebore was virtually rendered by the

proposedschemethatsatisfiestheRankine-Hugoniotjump

condi-tions.

4.3. Wavedeformationbyasubmergedshoal

Inthistestcase,theUnSWASHmodelisappliedtostudywave

Fig. 2. Comparison of free surface (top panel) and velocity (bottom panel) profiles at t = 7 s for dam break over wet bed obtained with different schemes and analytical solution.

experimental arrangement conducted by Berkhoff et al. [48] has

servedasaclassicaltestcaseforthisstudy.Thesimulationis

con-sideredinabasinof20mwideand35mlongwithaplaneslope

of1/50onwhichanellipticshoalisrested.Amonochromatic,

uni-directionalwave withaperiodof1sandaheightof4.64cm

en-terstheareaandallowstopropagatethroughoutthedomain.

De-tailedinformationonthebathymetryandthewaveconditionscan

befoundelsewhere(e.g.[30,48]).

Thedomainwasmeshedataresolutioninbetween0.03mand

0.07m,withatotalofapproximately370,000triangularelements.

Themeshisquiteuniformasthedominantwavesarepresent

ev-erywhere inthedomain. Furthermore,itisrelatively coarse with

20to50gridcellsperincidentwavelength,whilethewavelength

becomes smaller on the slope and shorter waves are generated

over the shoal. The time step wastakeninitially as 0.005s. The

maximumCourantnumberwassetto0.8andthesimulationtime

equaled30s.

First,afewwordsconcerningthecomputationalefficiency.The

performance oftheunstructured-meshversionofSWASH is

com-paredto thestructured-grid version. Theserialmode simulations

were run on an Ubuntu 16.04 desktopwith a 64-bit Intel Xeon

processor (2.3GHz).The executiontime ofUnSWASH neededfor

thecurrentcaseappearedtobe about108CPUminutes,withthe

highestlevel ofresourcerequiredby thesolution ofthepressure

Poisson equation. Furthermore, the total CPU time per grid cell

per time step required for UnSWASH was approximately 2.7 μs,

whereasfortheregular-gridSWASHwasabout1.8μs.Thisincrease

in thecomputationtime is probablydueto theoverhead related

tounstructuredmeshdatastructureaccess,causinglowcache

ef-ficiency.Itisclearthathigh-resolutionUnSWASHmustutilizeHPC

resourcesinordertobescalableandefficient.Aparallelcodewill

be build usingmessage-passing paradigms andwill be testedon

a commodity computer cluster. This will be reported ina future

paper.

Attentionis turnednext to thecomparisonbetweenthe

com-putedandobservednormalisedwaveheightsshowninFig.3.

Wave heights were measured along eight transects atregular

intervals(see [48]fordetails),of whichfourinteresting ones are

consideredinthepresentstudy.Thewavefocussinginthewakeof

theshoal(transects2and4),thewave shoalingandrefractionon

theslope(transects6and7),andtheinterferencepatterncaused

bydiffraction(transect4)arewellcapturedbythemodel,despite

theuseoflow-order discretizations.Inaddition, theseresults are

ingoodagreementwithpreviously computedresultsdiscussedby

StellingandZijlema[30].

4.4. Breakingwavesoverabarredtopography

ThelaboratoryflumetestofBoers[49]isconsideredwith

ran-domwavespropagatingoverabarredsandybeach(seeFig.4).

The case 1C is discussed with waves generated at the

wave-maker based on a target JONSWAP spectrum with a significant

waveheightof0.103mandapeakperiodof3.33s.Thiscasehasa

relativelylowwave steepness,whilewavesshoaloverthesloping

bottomuntilthey breakinthesurfzone. Fordetailedanalysison

thistestcaseincludingtheinfragravitywave dynamics,see

Rijns-dorpetal.[47].

The computational domain is 32 m long and 1 m wide. The

non-uniformunstructuredgrid employed inthe presenttest case

consistsofapproximately12,000triangles withthecellsize

vary-ing in between 0.2 m offshore and 0.01 m near the coast. This

mesh resolution should provide an economical representation of

the bathymetricvariability in the considered area, while the

nu-merical approach must remain accurate in the presence of

rela-tively large depth gradients. An initial time step of 0.004 s was

taken with a maximum Courant number of 0.5. This time step

turnedout tobe unchangedthroughoutthesimulation.Since the

Fig. 3. Computed (solid line) and measured (circles) relative wave heights along four transects for the wave over submerged shoal. H 0 is the incident wave height.

Fig. 4. Bottom topography and location of the selected wave gauges of the experiment of Boers [49] .

Fig. 5. Computed and measured significant wave heights (left panel) and mean wave periods (right panel) along the flume for Boers 1C. Bullets: experimental data; solid line: UnSWASH.

front approximation(HFA) ofSmit etal.[46] isemployed.

Imple-mentationdetailscanbefoundinAppendixB.Furthermore,a

fric-tioncoefficient c_f=0_.001wasadopted.Attheoffshoreboundary,

a weakly reflective condition is imposed based on the recorded

time seriesofsurfaceelevationatthemostseawardlocatedwave

gauge.

Fig. 5 compares the computed and observed significant wave

heights Hm0=4√m0 and mean wave periods Tm01=m0/m1 of

shortwaves.

The momentsmn= fnE

(

f

)

df are computedfromtheenergy

density spectra E(f) (with f frequency) of the surface elevation.

Clearly, thetrendoftheintegral waveparameters throughoutthe

domainiswellreproducedbythemodel.

Next, we compare the computed wave spectra with the

ob-servedonesinFig.6.

Thelocationoftheselectedwave gaugesisindicatedinFig.4.

Boththegenerationofhigherharmonicsandthetransferofenergy

towards low frequencies are properly predicted. Also, the energy

lossatthemid-tohigh-frequency rangeduetowave breakingis

resolved quiteaccurately.In general, themodelis ableto predict

the primary characteristics ofthe spectral evolution, both in the

shoalingregionandthesurfzone.

Though the potential of the present low-order unstructured

staggeredmeshmethodisclearlydemonstrated,theresultscanbe

Fig. 6. Observed (thick line) and predicted (thin line) energy density spectra at selected gauges of shoreward propagating waves for case 1C of Boers [49] .

Fig. 7. Schematic views of the conical island experiment. Plane view of the wave basin and the island (top panel) and side view of the island along section A - A (bottom panel).

unstructured-meshversionofSWASHtoincludeverticalresolution

willbedoneinthefuturework.

4.5. Runupofwavesonaconicalisland

Inthislasttestcasetherunupofasolitarywavearounda

con-ical island is discussed. The experimental configuration ofBriggs

etal.[50] isselected andisdepictedinFig. 7.Surface elevations

wererecordedusingwavegaugesateightdifferentlocationsas

in-Fig. 8. Non-uniform mesh used for the conical island test case with local refine- ment towards the island.

dicatedinthefigure.Thisbenchmarkischallenging,becauseofthe

wettinganddryingcycleswhilewavesrunupanddownonthe

is-land.

The grid used to obtain accurate wave runup andinundation

onthe islandcontainsapproximately171,000triangularcellsand

Fig. 9. Comparison with the data of Briggs et al. (1995) for the runup on a conical island: H = 0 . 045 d (left panel), H = 0 . 096 d (middle panel), and H = 0 . 181 d (right panel). Time histories of surface elevation at selected gauges around the island are shown. Dashed line: experimental data; solid line: UnSWASH.

Thelargestgridsizeawayfromtheislandis0.4m,whereasthe

smallestsizeisaround0.01m.Subsequently,asolitarywavewith

givenheightHwasimposedatthewesternboundary.

Three wave conditions are considered in the present study:

H=0.045d,H=0.096d andH=0.181d withd=0.32mthe still

waterdepthofthebasin(cf.Fig.7).Thecorrespondingcaseswere

simulated with UnSWASH using the following settings. The

ini-tial time step was 0.001 s with the maximum CFL of 0.8 for

the first two cases, whereas for the last severe case they were

0.0005 s and 0.5, respectively. The simulation period wasset to

25 s forall cases. Since wave breaking wasobserved on the lee

side of the island for the case H₌0_.181d_, the HFA is applied

accordingly.

Fig.9providesacomparisonbetweenthemodelresultsandthe

measureddataattheselectedgauges.Theobservedsurface

eleva-tionsareverywellpredictedbytheUnSWASHmodel.

In particular, thecomputed arrival time andheight of the

in-coming wave agrees with the measured data.There is, however,

a phase shift observed in the peak at gauge 22, especially for

the case H=0.181d. This known issue can be resolved by

in-cluding two vertical layers inthe model.Finally, itis worth

not-ingthat themodelresultsareconsistentwithpreviouslyobtained

results for the structured mesh case [24]. The associated mesh

is uniform with a grid size of 0.05 m and the number of grid

cells is approximately twice as large asthat of the unstructured

grid.

5. Conclusions

An extension of thewave-flowmodel SWASH tounstructured

triangular meshes has been discussed. The main motivation for

theapplicationofsuchmeshesforthesimulationofwave

dynam-icsistheeaseoflocalgridrefinement.Thecovolumemethodhas

beenadoptedforthespatialdiscretizationofthedepth-integrated

shallowwaterequationswiththeprimitivevariables.The velocity

componentsnormaltothecellfacesareemployedastheprimary

unknownsinthediscretization.Toaccountforwavedispersionthe

non-hydrostaticpressureisincluded.

Thecovolumemethodhasthegreatbenefitofallowingto

con-struct finite difference discretizations that mimic desirable

prop-erties of PDEs (e.g. topology, conservation, jump relation) while

keepingthe computational stencilcompact. Thissignificantly

im-provesboth therobustnessandtheefficiencyofSWASH.In

addi-tion,mimeticdiscretizationsroutinely enablephysically

meaning-fulresultstobeobtainedonrelativelycoarsemeshes.

Still,theapplicationofthecovolumeschemeispractically

lim-itedtoDelaunay-Voronoimeshes,whichmayimpedetheuser

flex-ibilitytogenerateadequategridscomprisingthenecessary

refine-ments.However,theassociatedorthogonalityrequirementisfound

nottobealimitingfactorinwave-relatedapplications.Thisis

ex-plainedbythefactthattherequiredmeshresolutionisreasonably

non-uniformforthescaleofwave dynamicsacrosstheentire

do-main.

A robust and efficient upwind-biased scheme has been

pro-posed.TheschemecomplieswiththeRankine-Hugoniotjump

re-lationsandis specificallydesignedwitha view topreservingthe

local momentum flux. This iscrucial to the simulation of

break-ingwavesandunsteadybores.Thishasbeendemonstratedbythe

dambreaktestcase.

Theproposed methodisgenerallyfirstorderaccurateon

non-uniform meshes. Despite this fact, the model is capable of

re-producing essential wave processes, including shoaling,

refrac-tion, diffraction, nonlinear interactions, depth-induced breaking

and wave runup, owing to the mimetic nature of applied

dis-cretizationsandtheuseoflocallyfinegrids.Thishasbeenverified

DeclarationofCompetingInterest

Theauthordeclaresthat hehasnoknowncompetingfinancial

interests orpersonalrelationshipsthatcouldhaveappearedto

in-fluencetheworkreportedinthispaper.

CRediTauthorshipcontributionstatement

MarcelZijlema:Conceptualization,Methodology,Software,

Val-idation,Formalanalysis,Writing-originaldraft.

Acknowledgements

The authorthanktheanonymousreviewersfortheir

construc-tive comments and suggestions. The author also thank Marien

Boerswho provideduswiththedataofthelaboratory flume

ex-perimentsandHaiyangCuiforprovidingtheunstructuredmeshof

theconicalislandtests.

AppendixA. Proofofmomentumconservation

Theobjectiveofthisappendixistodemonstratethatrecurrent

relation (12) does not produce a momentum conservation error,

givenauniformbedandazerobedfriction.UsingEq.(7),we

ob-tainthefollowingfinitedifferenceform



sf duf dt +



sc L, f Ac_L



Qk1 ˜ hf

uf− uc_L ,k1

+Qk2 ˜ hf

uf− uc_L ,k2



· nf =−

Pc_R− Pc_L

(14) where Pc= 1 2h˜f



gh2 c+hcpc



is thedepth-averaged totalpressureatthe circumcentre.Discrete

momentum conservation can be expressed as therate ofchange

inthe totalamountofmomentum huwithin a controlvolume is

due onlyto the net flux throughthe edges of the volume.Since

thenormalfacevelocityuf istheprimaryunknown,interpolation

must be employed toobtain the momentum at a singlelocation

within themesh.Tothisend, themeshcellistreatedasthe

con-trolvolume.Subsequently,wederiveanequationforthecell-based

momentum vector using Perot’sinterpolation scheme(8). Finally,

summation ofits flux contributionsover thecell faces mustlead

to adiscreteequivalentofthemomentumequation indivergence

form, whichcompletes the proof. Thisformal procedureof proof

utilizes thegeometricpropertiesofthe meshandisalso

applica-bletoCartesianstaggeredschemes(see,e.g.[14,33]).

First,Eq.(14)isrewrittenas





scL, f+



scR, f



duf dt +





scL, fccL+



scR, fccR



· nf =−

(

PcR− PcL

)

with cc= 1 Ac  k∈Sc

α

c,kQ_˜k hf



ˆ uk− uf



(15)

the cell-based discretization of the advection term evaluated for

the two cellsadjacent toface f.Note that the transported

veloc-ityuˆ_k isinterpolatedfromfacenormalcomponentsfromthecell

upwindoffacek.Asaconsequence,cc_R· nf=0ifuf>0,and

like-wise, cc_L· nf=0ifuf <0.Next,thediscretized equation is

mul-tiplied by thenormal ofthe facenf andits length lf,and

subse-quentlysummedoverthefacesofcellc

 f∈Sc nflf



sc_L, f+



sc_R, f

_duf dt +  f∈Sc lf



sc_L, fcc L+



scR, fccR

· nfnf=−  f∈Sc nflf

Pc R− PcL

Wenowdemonstrateconservationofmomentumby

transform-ing the above equation into an equation for the cell-based

mo-mentum vector u. Each interior face is shared by two triangular

cellsofwhichonecontributestothecellunderconsideration.

Fur-thermore,attheboundaryfaceafluxofmomentumisprescribed.

First,therateofchangeinmomentumcanberecastedas

 f∈Sc nflf



sc, f duf dt = d dt  f∈Sc nflf



sc, fuf=Ac duc dt

followingfromEq.(8).Thisholdsforallcellsinthecomputational

mesh.

Next, the advective acceleration termcan be rearranged with

theaidofthefollowinggeometricidentity



f∈Sc

lf



sc, fnfnf=AcI

with I the identity matrix. This identity follows from the

diver-gencetheorem [14]. Then,using Eq.(15),theadvection termcan

beexpandedas  f∈Sc lf



sc, fcc· nfnf=cc·  f∈Sc lf



sc, fnfnf=cc· IAc=Accc = k∈Sc

α

c,kQ_˜k hf



ˆ uk− uf



which conservesmomentum in the considered cell since

α

_c

R,k=

−

α

c_L,k at each interiorface k. Thisimplies that thesum of

con-tributions to theleft andright cellscancelsthe advectiveflux at

face k. This holds for all cells except the boundary cells where

thechange inmomentum isdueto the fluxesacrossthe

bound-aryfaces.

Finally,the pressureterminthe interiorcell canbe rewritten

as  f∈Sc nflfPc=−Pc  f∈Sc nc, flf=0

by notingthat the pressuregradient isaligned withtheflow

di-rection and the cell has a closed surface. In case the cell under

considerationisaboundarycell,then

 f∈Sc nflfPc=  f∈Sc nc, flfPf

withP_fthecontributiontotheprescribedmomentumflux.

Inshort,byrewritingthefinitedifferenceform(14)intoa

dis-crete equation in divergence formwithout introducing anyother

approximation,we haveshownthat itdoesnot createordestroy

momentumineachindividualmeshcellofthecomputational

do-main.Thisimpliesthat thecomputed amountofmomentum can

onlychangeasaresultofanon-zeronetmomentumfluxoverthe

boundary of the domain, a non-uniform bed, or a non-zero bed

friction.

AppendixB.Hydrostaticfrontapproximation

Dueto coarseverticalresolutions some measuresare required

toproperlyapproximatewavebreaking.Theyareintroducedwith

theHFA techniqueof Smitet al.[46]. Inprinciple, toinitiate the

onsetofthebreakingprocessatlowresolutions,asteepbore-like

wavefrontistrackedandsubsequentlyahydrostaticpressure

dis-tributionisimposedlocally.Furthermore,theHFArequiresthe

in-clusion of the diffusion termin momentum equation (4),as

fol-lows

∂

huf

∂

t +· · · =· · · +