The role of the Rankine-Hugoniot relations in staggered finite difference schemes for the shallow water equations

Delft University of Technology

The role of the Rankine-Hugoniot relations in staggered finite difference schemes for the

shallow water equations

Zijlema, Marcel

DOI

10.1016/j.compfluid.2019.104274

Publication date

2019

Document Version

Final published version

Published in

Computers & Fluids

Citation (APA)

Zijlema, M. (2019). The role of the Rankine-Hugoniot relations in staggered finite difference schemes for the

shallow water equations. Computers & Fluids, 192, [104274].

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

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.

‘You share, we take care!’ – Taverne project

https://www.openaccess.nl/en/you-share-we-take-care

Otherwise as indicated in the copyright section: the publisher

is the copyright holder of this work and the author uses the

Dutch legislation to make this work public.

ContentslistsavailableatScienceDirect

Computers

and

Fluids

journalhomepage:www.elsevier.com/locate/compfluid

The

role

of

the

Rankine-Hugoniot

relations

in

staggered

finite

difference

schemes

for

the

shallow

water

equations

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 2 April 2019 Revised 14 August 2019 Accepted 23 August 2019 Available online 23 August 2019

Keywords:

Shallow water equations Finite difference schemes Rankine-Hugoniot relations Staggered grids

Odd-even decoupling

a

b

s

t

r

a

c

t

ThepurposeofthisworkistopointouttherelevanceoftheRankine-Hugoniotjumprelationsregarding thenumericalsolutionoftheinviscidshallowwaterequations.Toarriveatphysicallyrelevantsolutions inrapidlyvariedflow,itisofcrucialimportancethatcontinuityofmassfluxandmomentumfluxacrossa steadydiscontinuityisfulfilledatthediscretelevel.Byadoptingthisviewpoint,finitedifferenceschemes canbestudiedthatmaybewellsuitedtosolveshallow waterflowproblemsinvolvingdiscontinuities, whiletheyarenotbasedonacharacteristicdecompositionofthegovernedhyperbolicequations.Three schemesonstaggeredgridswitheither thewater level orthe waterdepthatthecell centreand the flowvelocityorthedepth-integratedvelocityatthecellinterfaceareexamined.Theydifferin(1)the characterofthe transport velocitytobias thediscretization ofthe advectiveaccelerationterm inthe upwinddirection,and(2)thedeterminationofthewaterdepthatthecellfacewithwhichthe depth-integratedvelocity must belinked tothe flow velocity.A detailed analysisis provided and aimed at highlighting the necessityof fulfillingthe Rankine-Hugoniotjump conditions forpreventing the odd-evendecouplingproblem.Theaccuracyandrobustnessofthreeselectedschemesisassessedbymeansof convergencetests,threeidealized1Dtestproblemswithexactsolutionsanda1Dlaboratoryexperiment ofthebreaking,runupandrundownofasolitarywaveonaslopingbeach.Numericalresultsrevealthat schemessatisfying exactlythejump conditionsdisplayimprovedperformanceoverschemes whichdo notsharethisproperty.Also,theseresultssupportstrongevidenceonthelinkbetweennotfulfillingthe jumpconditionsandtheappearanceofodd-evenoscillations.

© 2019ElsevierLtd.Allrightsreserved.

1. Introduction

The celebrated nonlinear shallow water equations have ob-tained widespread acceptancefor modelinga rich varietyof free surface flows in rivers, lakes, floodplains, coastal areas and con-tinental shelves. Such equations may lead to the formation and propagation ofhydraulic discontinuitiessuch as hydraulic jumps, tidalbores,andbreakingwavesinshallowwater[1].Widelyused numerical methods for solving inviscid shallow water equations commonly respect the property of (strict) hyperbolicity of the equations, exploit the features of associated wave characteristics, andrelyupon the fullyconservativeformofthe governing equa-tions withthe waterdepth andthe depth-integrated velocity as conserved variables.These methodsare usually based ona time-explicit,finitevolumediscretizationemployingacell-centred colo-cated grid arrangement with the aim of computing the inviscid flux atthemidpoint ofacell face.Such methodsare appropriate forcapturing discontinuitiesinrapidlyvaried flows. Inparticular,

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

theGodunov methodemploying approximateRiemann solvers of Roe [2]and HLL [3]has proven successful. A survey of this vast subjectisgiveninLeVeque[4]andToro[5].

Althoughhyperbolicsystemsadmitdiscontinuoussolutions, hy-perbolicityis usually the source of problemsthat arise in devel-oping a numerical scheme forthe solution of the shallow water equations.Onedrawbackofthispropertyisrelatedtothe numeri-calimbalancebetweenthepressurefluxandthetopographyterm thatmayappearinthecaseofdiscontinuousbedtopography.This imbalanceleadstoartificialflowsacrosssharpbottomgradientsin waterinitiallyatrest.Thisnumericalartifact maybeovercomeby constructinga well-balancedscheme that preservescertain kinds ofsteadystate solutions [6–8].A varietyof numericaltechniques toachievethewell-balancedpropertyintheRoeandHLLmethods isproposedin[9–13].However,thedevelopmentofwell-balanced schemescanberatherdifficultandsuch schemesareoftenfound tolackrobustnessinreal-lifeapplications.Anotherdifficultyisthat stricthyperbolicityislostifadrybedispresent[14].Especiallythe Roe’s Riemann solver isnot able to cope withproblems that in-volvewet-dryfronts.Inaddition,dedicatedmeasuresarerequired toavoid non-negativityof thewater depthin the casea wet

re-https://doi.org/10.1016/j.compfluid.2019.104274

gionfallingdry.Theytypicallymodifythecellvariablesatthecell faceandsubsequentlytheinterfacefluxbymeansofahydrostatic reconstruction[7].

The numerical methods mentioned above target to exploit the hyperbolicity property of the governing equations featuring the propagation of wave-like disturbances along the character-istics with a finite speed. Where these characteristics intersect corresponds to a shock or discontinuity, which must satisfy the Rankine-Hugoniot jump relations [1]. They describe the relation betweentheflowstatesonbothsidesofthediscontinuitythat fol-lowfromconservationofmassandmomentum.The jump condi-tionsusuallyexpresscontinuityofmassfluxandmomentumflux acrossasteadydiscontinuity(i.e.theshockspeediszero).Inthis paperwewillshowthatcomplyingwiththesejumpconditionsat discretizationlevelwarrantsthe useofnumericalschemeswhich are suitable to solve problems involving (steady) discontinuities, whereas they do not explicitly utilize the hyperbolicity property oftheunderlyingequations.Inthiscontext,theflowvariablesare assumedtobeconstantineachgridcellorcontrolvolume,while treatedasdiscontinuousacrosstheedgeofthecell/controlvolume. Yetitisofcrucialimportancethattheconstancyofmassand mo-mentumfluxesacross thediscontinuityis fulfilledatthediscrete level.

Motivatedbythisobservation,anyfinitedifferenceschemecan, inprinciple,beemployedforthemodelingofrapidlyvariedflows providedthatthealgebraicRankine-Hugoniotrelations,expressing the continuity of mass flux and momentum flux across the vol-ume edge, are fulfilled. The present study deals with finite dif-ference schemes on staggered Cartesian grids. Such schemes are an attractivealternative to colocatedfinite volumeGodunov-type schemes.Thereasonsarethreefold. First,staggeredschemestreat theinviscid flux componentwise and the hyperbolicityfeature is ignored.Second,theassociatedspatialdiscretizationprovides sim-pleandaccurateapproximationsandcaneasilybecombinedwith the implicitADI-type splitting schemes [15,16] andsemi-implicit

θ

-methods[17,18],whicharethewidelyemployedtimeintegration techniquesduetotheir largestabilityregionanddelicatebalance betweenaccuracyandlowcomputationalcost.Third,thestaggered approachissuccessfullyemployed incoastal andoceanmodeling of,amongothers,large-scalebarotropicflows,stormsurges, strati-fiedflows,floodwavesinrivers,inundationbytsunamiwaves,and morphodynamicsincoastalareasin2Dand3D(see,e.g.[18–22]).

Staggered schemesarenot onlyapplicable tosubcriticalflows, butalsosuitableforrapidlyvariedflows.Forexample,thevarious staggeredgrid approachesof ZhouandStansby[23],Stelling and Duinmeijer[24],Madsen etal.[25],KramerandStelling[26],Cui etal.[22], andDoyen andGunawan [27] to solve shallow water equationsyield accurateresultsforflowsathighFroudenumbers includinghydraulicjumpsanddambreakproblems.Inthecurrent workitisexplainedindetailwhythesemethodsareadequatefor suchflows.

The objective of this paper is twofold. Firstly, it addresses a numberoffinitedifferenceschemesonstaggeredgridsforthe so-lution of the one-dimensional, depth-integrated, inviscid shallow waterequations.These 1Dequations servedasa basisfor identi-fyinga number ofnumerical issues to be discussed. The current work aims to analyzein-depth the spatialdiscretization ofthree staggered schemes, each with a different set of dependent vari-ables, either primitive or conserved, in relation to the Rankine-Hugoniotjumpconditionsholdingacrossdiscontinuities. Thiswill shedlightonconservationpropertiesofthepresentedschemes.In thisrespect,itwillbeshownthattheincompressibilityconstraint playsakeyroleinachievingenhancedaccuracyandrobustness.

The second goal of the paper is to demonstrate how the Rankine-Hugoniot jump conditions can help suppress odd-even spuriousmodesthatwouldotherwiseoccurinstaggeredschemes.

Inparticular, theaimis toshow thattheseconditions are neces-sary forpreventing theodd-even decoupling problem. This turns out to be beneficial for the numerical accuracy and stability of the solution.This forms the basis of ourclaim that fulfillingthe Rankine-Hugoniotrelations is to be preferredover exploitingthe hyperbolicity property of the governing equations with which a goodfinitedifferenceschemeshouldcomply.

Therestofthispaperisorganizedasfollows.Section 2briefly reviewstheshallowwaterequationswewishtosolve.Inthis con-text,itisimportanttoemphasizethatincompressibilityisakey el-ementintheformulationanddiscretizationofthegoverning equa-tions.Inparticular,theanalogywithcompressibleEulerequations ofgasdynamicsiscriticallyreviewed.The spatialdiscretizationis treatedfirstfora1Dgeneric continuityequation inSection 3and then,in Section4,forthe 1Dshallowwaterequations,by means of astaggered approach, including additionalanalysis concerning the Rankine-Hugoniotrelations. Section 5 outlines the aspects of time integration relevant forthe scope ofthe presentstudy. The performance ofthreeselectedschemesis assessedfirstby means ofconvergencetests andthenby threeidealized 1D tests involv-ing flow discontinuities, for which the conservationof mass flux andmomentumfluxiscrucial.Inaddition,a1Dlaboratory exper-imentofthe breaking,runupandrundownofasolitary wave on a sloping beach is usedin the assessment.Numerical results are presentedtodisplaythecomparisonbetweenthevariousschemes. These assessmentsand outcomesare dealtwithin Section 6. Fi-nally,inSection7,thefindingsofthisstudyaresummarized, fol-lowedby a few practical comments upon theimplementation of thediscussedschemesintwodimensions.

2. Governingequations

Weconsiderthe shallowwaterflow withafree surfaceunder gravity in a Cartesian coordinatesystem x=

(

x,y,z

)

. Thisflow is confined by the surface elevationz₌

ζ

(

x_,y_,t

)

(positive upwards with respect to the datum level z=0 with z-axis in the direc-tionoppositetothatofgravity)andthebedelevationz=−d

(

x,y

)

(positivedownwardswithrespecttothedatumlevel),andis uni-formindepth.We restrictourselvestoincompressibleflows gov-erned by the depth-integrated continuity and momentum equa-tions.

Thecontinuityequationis

∂

h

∂

t +

∇

· hu=0 (1)

whereh

(

x,y,t

)

=

ζ

+d is thewaterdepthandu

(

x,y,t

)

=

(

u,

v

)

is the velocityvector withthe depth-averaged componentsu andv

alongthexandycoordinates,respectively.Themassfluxisgiven byhu.Notethatinthetwo-dimensionalspace(x,y),thedivergence operator

∇

reducesto(

∂

/

∂

x,

∂

/

∂

y).Eq.(1)iswrittenindivergence formandexpressesconservationofvolumeintheincompressible flow.TheLagrangianformofEq.(1)reads

Dh

Dt =−h

∇

· u (2)

withD/Dt=

∂

/

∂

t+u·

∇

denoting thematerial derivative. Hence, changesinthewaterdepththatmovealongwithuareassociated withdivergence ofthe flow velocity u. Thisis thecasewith, for instance,spatialvariationsoftheseabed.

Themomentumequationis

∂

q

∂

t +

∇

·



uq+1 2gh 2_I



_=gh

_∇

_{d (3)}

where q=hu is the depth-integrated velocity (momentum per unit area normalised by the density), g is the gravitational ac-celeration, and I is the identity matrix. The momentum flux is

givenby u q+1

2gh2I andthe pressureisassumedto be

hydro-static. Eq. (3) merely maintains its conservative character in the absence of the source term, i.e. the bed is supposed to be uni-form. It is worth commenting that the depth-integrated velocity

qin Eq.(3) andthe massflux hu inEq.(1)are differentintheir role:theformerisa storagequantity,whereasthelatterone isa transportquantity.

The continuity equation (1) and the momentum equation (3) are the well-known inviscid shallow water equa-tionsinstrongconservationformforthesolutionofhandqbeing the conserved variables. Initially smooth h and q may become discontinuouswithinfinitetime.Theseshockspropagateatspeeds determined by the Rankine-Hugoniot relations arising from the conservation of mass andmomentum [4]. Various Godunov-type methods have been devised for the numerical solution of this system of nonlinear hyperbolic equations that is mathematically equivalent to the compressible Euler equations of gas dynamics [4,5]. Applications of approximate Riemann solvers to shallow water equations, such as Roe [2] and HLL [3], are found in, e.g. [9,12–14,28]. Furthermore, alternative methods have been proposed that avoid the necessity to solve a Riemann problem. Examples are the TVD-MacCormack scheme [29], the central-upwind method [6], the artificial viscosity approach [30], and energy stableschemes using numericaldiffusion operators based on entropy variables [8]. Typically,the abovementioned methods employ colocated variable arrangement:both hand qare stored atthesamesetofgridnodes.

For reasonsthat become clearlater on, the preferredformof theshallowwaterequationsforincompressibleflowis

∂ζ

∂

t +

∇

· q=0 (4)

∂

hu

∂

t +

∇

·

(

qu

)

+gh

∇ζ

=0 (5)

withtheprimitivevariables

ζ

andu.Theseequationsarewritten inweak conservationformwithouta splittingofthesurface gra-dientintoapressureflux anda sourcetermthat includesthe ef-fect ofthebedslope(cf.Eq.(3)). Wethusrefrainfromusingthe hyperbolicityproperty. Nonetheless,forlarge-scaleapplicationsin coastalandoceanengineering,thegoverningequationsareusually writteninthefollowingnon-conservationform

∂ζ

∂

t +

∇

· hu=0

∂

u

∂

t +

(

u·

∇

)

u=−g

∇ζ

whereasotherrelevantforcesshouldbeincluded,suchasthe fric-tionalforces(windshear andbottomroughness)andtheCoriolis forceduetotheEarth’srotation.Differentnumericalmethodshave beenproposed forthesolutionoftheseshallowwaterequations; see,e.g. [15–21,31].Thesemethodsall haveincommonthat they employstaggeredgrids:theprimitivevariablesarestoredat alter-nate locations;the waterlevel is storedata cell centre,andthe velocityisstoredatacellface.Usingthisvariablestorage arrange-ment is an effective wayto avoidodd-even decoupling between thesurfaceelevationandflowvelocity.

Although the mathematical structure of Eqs. (1) and (3) is equivalent to that of the Euler equationsof gas dynamics, being asystemofhyperbolic conservationlawsforwhichGodunov-type methodsare applicable,acrucialdifferenceisthephysical behav-ior of the flow field. The compressibleEuler equationsare classi-fiedashyperbolicandexamplesofflowsgovernedbythese equa-tions are unsteady sub- and supersonicflows and steady super-sonic flows [32]. On the other hand, the inviscid shallow water

equationshavebeenderivedfromdepthintegrationofthe incom-pressibleEulerequationsandtheassociatedflowbehaviorisof el-lipticnatureduetotheincompressibility constraint[32].Muchin thesame way,Eqs.(4) and(5)are hyperbolic forunsteadyflows andbutellipticforsteadyflows.Examplesofthelatterarerapidly variedsteadyflows(e.g.flowoveraweircrestandflowthroughan orifice).1_{Yet,forsteadystatesolutions,thecontinuityequation(4)}

degeneratestoaconstraint

∇

· q=0,i.e.thedepth-integrated ve-locityisasolenoidalfield.Thisconstraintcanbedealtwith effec-tivelybytheuseofstaggeredgrids.Thisalsoappliesforunsteady flows. Examples of (un)steadyflows withsudden changes in the bathymetrywillbepresentedinSection6.3.

Furthertothis, Eq.(3)isa transport equation forthevariable

q,whichistheadvected quantity,whereasuactsastransport ve-locity.ThisisalsoreflectedinthesecondtermofEq.(3),

∇

·

(

uq

)

=

(

u·

∇

)

q+h

(

∇

· u

)

u

wherethefirsttermontherighthandsideisrelatedtothe back-groundflow. The second term is closely linked to the change in the water depth along flow paths (viz. Eq. (2)), illustrating how the changein the depth-integrated velocity q andthe change in thewaterdepthhareintertwined.YetEqs.(1)and(3)areasetof hyperbolicequations,andthus attractiveasthey lend themselves todiagonalizationyieldingdistinct variablesalongtheir character-istics.

Eqs.(4)and(5),ontheotherhand,naturallydescribethe prop-agationofgravitywavesastheessentialtermsarethesurface gra-dient,i.e.thethirdtermofEq.(5),andthedivergenceofthemass flux, i.e.the second termof Eq.(4). Observethat the divergence termofEq.(5)includestwoparts,(q·

∇

)u and(

∇

· q)u,ofwhich thesecondpartiscommonlyrelatedtothewavepropagation(viz. Eq.(4)).Hence,thecombinationofterms(

∇

· q)uandgh

∇ζ

inthe momentum equation characterizes the wave-like propagation on topofthe backgroundflow. Thefirst partofthedivergenceterm representsadvectioninthebackgroundflow.Thisimpliesthatthe depth-averagedvelocityu=q/histransportedbythemassfluxq, likevolume (ormass). Thistransportedvelocity uis thusthe re-solved variable of the depth-integrated momentum equation (5). Thisfeatureturnsouttobebeneficialforaccuracyandstabilityas willbedemonstratedinthecurrentstudy.

InthispaperitisarguedthatthemodelbasedonEqs.(4)and (5)canpreserve shocksanddiscontinuities, iftheunderlying nu-mericalschemecomplies withsomeformofdiscreteconservation thatenablestocomputeshockswiththecorrectspeed.Yet,rather thanconservation itself,it wouldbe crucialforphysical accuracy to employ a discrete form ofthe Rankine-Hugoniot jump condi-tions whichrelate the fluid stateson both sidesof a discontinu-ity.Thereby,a finite differencescheme that reliesupon thistype ofconservationformisabletoproducethecorrectweaksolution totheaforementionedequations.Thisisthecentralthemeofthis paper.

3. Finitedifferenceschemesfora1Dgenericcontinuity equation

Before we outlinevarious numerical methods forthe solution of the shallow water equations, we first discuss the spatial dis-cretizationofagenericcontinuityequationtodemonstratethe ba-sicideasofthedifferentnumericalapproximationstofollowwith someusefulobservations.The numericalapproachto theshallow waterequationswillbedevelopedlaterinSection4.

1 The associated pressure distribution may not be hydrostatic in the region of curved flows.

Fig. 1. A part of a 1D computational grid.

Agenericcontinuityequationisaspecialcaseofaconservation law(see[4]),andisgiveninone-dimensionalformasfollows

∂

c

∂

t +

∂

uc

∂

x =0

wherec(x,t)isaconservedvariableandu(x,t)istheflowvelocity. Thisequation typically expresses conservation of quantity c (e.g. mass,volume, internal energy), implying that the rate of change of its amount within a closed surface is merely due to the flux

f

(

c

)

=uc across thesurface.Ifthenetfluxiszero,likeinan iso-latedsystem,thencdoesnotchangeovertime.

We define a computational grid inthe x_−space discretized in

M grid cells. A part of this computational grid is depicted in Fig.1. Within thisgrid we have cell centres located at xm, with

m=1,...,M,andcellinterfacesatxm+1/2,withm=0,...,M.We

assume an equidistant grid with a constant grid size



x. The gridcellscanberegardedascontrolvolumescontainingdensities

cm(t)≈ c(xm,t).Theirdiscrete valuesrepresentpiecewiseconstant

averagesofthedensitiesoverthecells.Consequently,thedensities arediscontinuousacrosscellinterfaces.

The desiredsemi-discretizedequation mustreflect the conser-vationpropertyinthesense thatcmcanonlybereducedby

mov-ing some of it to neighboring cells. Or its amount can be en-larged by moving something from neighboring cellsinto cell m. Hence,thetemporalvariationofcmcanbeseenasin-andoutflow

throughthecellinterfaces.Thisisexpressedintermsofafaceflux

fm+1/2=f

(

cm+1/2

)

≈ f

(

c

(

xm+1/2,t

))

andnaturally,wegetthenext

balanceprinciple



xdcm

dt = fm−1/2− fm+1/2 (6)

whichyieldsthefollowingsemi-discretizedformulation dcm

dt +

fm+1/2− fm−1/2



x =0

Notethatthebalanceequation(6)canalsobederivedfromthe ap-plicationofthefinitevolumemethod.Thismethodisdesignedto obtainacorrectweaksolutionthrough,amongothers,the conser-vation requirement [4,33].Although global conservation can eas-ilybeobtainedwithanyfinitevolumescheme,duetothenatural cancellationoffluxes inthe interiorof the domain,the question remainshowtodeterminethenumericalfluxatthecellinterface, asthedensitycisnotknownatthecellinterfaces.Thisisthe cen-tralissuebeingaddressedinthefinitevolumemethod:toexpress thecell-facevalueintermsofcell-centrevalues.Interpolationthus determinesthevariationoftheconservedquantity cbetweencell centroids.Below,we willseethat thechoiceofinterpolationwill affectbothaccuracyandstabilityofthenumericalmethod.Thisis particularlythecaseforhyperbolicPDEs withdiscontinuous solu-tions.Onthecontrary,globalconservationprovestobemore rele-vanttoaconvection-diffusiontypeequation forwhichtheoption offlux interpolation is lesscriticalanyway. This is especiallythe casewhen there are no sources and sinks. In such a case, tele-scopiccancellingofinternalfluxeswouldberatherbeneficial.

By employing a staggered variable arrangement the discrete scalarscareresidedincellcentresasdefinedbefore,butthe dis-cretevelocitiesuarepositioned atthecellfaces. Theirvaluesare designated by um+1/2, with m=0,...,M. To demonstrate proper

treatment of a discontinuity at the discrete level, we consider a

stationaryflow.This caseismotivatedby therequirementofflux conservation, as becomes clear later on. We have the following spatialdiscretization

fm+1/2− fm−1/2



x =

um+1/2cm+1/2− um−1/2cm−1/2



x =0 (7)

Thedensitycatthefacesofthecontrolvolumeneedstobe inter-polatedintermsofcell-centrevalues.Linearinterpolationor arith-meticaveragingisthecommonlyusedoneandisgivenby cm+1/2≈

1

2

(

cm+cm+1

)

Substitutionyieldsasecondorderschemeforauniformgrid,

0= fm+1/2− fm−1/2



x ≈

um+1/2

(

cm+cm+1

)

− um−1/2

(

cm−1+cm

)

2



x

(8)

This scheme can give rise to a spurious checkerboard mode for densityc.Indeed,iftheflowvelocityisdivergencefreeoruniform, i.e.um+1/2=um−1/2,weobtain

cm+1=cm−1

which involves alternate c values only, allowing unwanted oscil-lations witha wave length of 2



x,cm=

(

−1

)

mc˜with dc˜/dx=0.

Thisodd-evendecouplingisnotdetectedbytheschemeand, asa consequence,thecheckerboardmodecanbeunbounded.

Topreventspuriousoscillations,afirstorderupwindschemeis chosen,asfollows cm+1/2≈



_c m, ifum+1/2>0 cm+1, ifum+1/2<0 (9)

Letusconsiderapositiveflow.Oneobtains

0= fm+1/2− fm−1/2



x ≈

um+1/2cm− um−1/2cm−1



x (10)

whichclearlycuredtheodd-evendecouplingproblem.

We reexamine the schemes (8)and (10) butfrom a different pointofview,namelyconservationofaflux.Certainphysical quan-tities may undergo rapid changes leading to the formation of a shock in a finite time. For instance, discontinuous water depths and flow velocities may occur in rapidly varied flows like flows oversteepbottomslopes,hydraulicjumps,tidalbores,andbroken waves.Ontheother hand,continuityisrequiredregardingtothe flowortransportacrossan abruptchange.Thiscanreadilybe de-ducedfromtheRankine-Hugoniotrelationorjumpcondition[1].

Letusconsiderastationaryjumpatx1<xs<x2,wherethe

sub-scripts1and2denotejustupstreamandjustdownstreamofthe jump,respectively.Althoughthe quantitiescandumayexhibit a discontinuity at xs, flux uc must be continuousacross the jump,

sinceduc/dx=0.Hence,thefollowingfluxconservationmusthold

u1c1=u2c2 (11)

Commonly,condition(11)enforcesconservationoffluxacrossthe shock. However, this conditioncan also be interpreted as a pre-requisite foran accurate andstable approximation ofthe flux at

thecellinterface.Bynature,discreteunknownsarepiecewise con-stantfunctionsonacomputationalgrid,andconsequently,are dis-continuous across cell faces. Generally, their smoothness proper-tiesare influencedbythe gridresolution.Inturn,thismayaffect the numerical accuracy. Variations inu and c, particularlyacross shocks, tend to be far stronger than changes in flux uc. Further-more, since grid cells on either side of a discontinuity contain no informationabouteach other,any approximationof unknown

c usingarithmetic averaging, linear interpolations or central dis-cretizations can give rise to spurious oscillations. Yet it is desir-ablethatEq.(11)holdsexactlyatthediscretelevel;giventhe up-stream state ofthe flow, thedownstream state iscompletely de-terminedbytheRankine-Hugoniotrelationinquestion.Hence, nu-merical errors associated withinadequate resolution of gradients arereducedsignificantly,whereasstabilitymightenhanced asthe odd-evendecouplingofanydependentvariablemaybeprevented. However,aswewillseeinduecourse,whenanumericalscheme doesnot complywiththe jumpcondition(11),theodd-even de-couplingproblemcropsup,ifallowedbythescheme.

Letusgobacktoourschemes.FromEq.(10)wehave um+1/2cm=um−1/2cm−1

implying that flux uc is continuous when evaluated from either sideofthecell.Thisschemestrictlyconservesfluxcellbycelland remains accurate irrespective of thegrid size. This is called cell-wise (or pointwise) conservationand theassociated jump condi-tionisgivenby

[uc]m=um+1/2cm− um−1/2cm−1=0 (12)

withtheoperator[· ]mdenotingthejumpofafluxacrosscellm.

Bycontrast,fromEq.(8),wederivethefollowing

[uc]m=um+1/2cm− um−1/2cm−1=um−1/2cm− um+1/2cm+1 (13)

Obviously,continuityoffluxucacrossthecellcannotbeobtained by the achieved scheme asthe right hand side gives rise to in-ternal sources of density c. This is the case, for instance, if the flow velocity u isdivergence free anddensity c exhibitsan odd-even oscillation.Toappreciatethe natureoftheunderlying prob-lem,itshouldbepointedoutthatthemagnitudeoftheerrordue todisbalance(13)canberelevant,especiallywhenacoarsegridis employed. Providedthisscheme isconverging,one canobviously achieveahighaccuracybyrefiningthemesh,butthisapproachis untenableinlarge-scalesimulations.

4. Finitedifferenceschemesforthe1Dshallowwater equations

Thissectiondealswiththevariousspatialdiscretizationsofthe inviscid,incompressibleshallowwaterequations.Inordertomake thematterclearandunambiguous,butwithnolossofgenerality, we willpresenta detailedstudyofthedifferentmethodsapplied to theone-dimensionalcontinuity andmomentum equations.We will discusstwo versionsof theseequations inthe remainder of thispaper.Thefirstsystemofequationstobeinvestigatedisgiven by

∂ζ

∂

t +

∂

q

∂

x =0 (14)

∂

hu

∂

t +

∂

qu

∂

x +gh

∂ζ

∂

x =0 (15)

withtheprimitivevariables

ζ

(x,t)andu(x,t).Furthermore,h

(

x,t

)

=

ζ

(

x,t

)

+d

(

x

)

andd(x) representsthebedtopography.These one-dimensional equations follow from Eqs. (4) and (5), respec-tively. The second one isa set ofhyperbolic equationsbased on Eqs.(1)and(3),andisgivenin1Dform

∂

h

∂

t +

∂

q

∂

x =0 (16)

∂

q

∂

t +

∂



uq+1 2gh2



∂

x =gh

∂

d

∂

x (17)

withthefluxvariablesh(x,t)andq(x,t).

Eachoftheabove setsof1Dequationsserves asaproxy fora detailedanalysisofthedifferentnumericalapproaches tobe dis-cussed below. We will first consider the spatial discretization of continuityequation(14)andmomentumequation(15).Thereafter, thediscretizationofthehyperbolicequations(16)and(17)willbe outlined.TimeintegrationwillbepresentedinSection5.

With the staggered grid arrangement, the primitive variables are carriedat alternate grid points. Boththe waterlevel

ζ

m and

waterdepthhmarelocatedatthecellcentreoftheassociatedgrid

cell, xm, whereas thisgrid cell is bounded by the cell interfaces

atxm−1/2 andxm+1/2,wherethevelocitiesum± 1/2arelocated;see

alsoFig.1.Notethatthebedleveldmisresidedatthecellcentre.

A two-point stencil can be obtained for both the continuity equationandalltermsofthemomentumequation,exceptthe ad-vection term, approximated with second order accuracyon stag-gereduniform meshes. This provides theability to construct nu-mericalschemesthatareconsistentwiththeRankine-Hugoniot re-lationsforconservationofmassfluxandmomentumflux.

Letusconsiderthecontinuityequation(14).Central discretiza-tioninspaceyields

d

ζ

m

dt +

qm+1/2− qm−1/2



x =0 (18)

Forastationaryflow,wehave qm+1/2− qm−1/2=0⇒

[q]_m=0 (19)

Eq.(19)representstheRankine-Hugoniotrelationwithrespect to the cellwise conservation of mass flux applied at the staggered grid.

Since um+1/2 is the primitive variable, we compute the mass

flux,asfollows

qm+1/2=hm+1/2um+1/2 (20)

The waterdepth atthe cellface hm+1/2 mustbe defined further.

Thefollowingfirstorderapproximationallowsconditional preser-vationofnon-negativityofthewaterdepth(seeSection5) hm+1/2≈ ˆhm+1/2=



_h

m, ifum+1/2>0 hm+1, ifum+1/2<0

(21)

withthehatdenotingan upwindvalue. Alternatively,thesurface elevation

ζ

can be taken upwind instead of the water depth h, whilethe bedleveldis determined by meansofa min operator (see,for instance,[34]). This maylead toa moreaccurate result, particularlyatcoarsegrids.Inthisstudy,however,werestrict our-selvestotheuseofEq.(21).Thiscompletesthespatial discretiza-tionofthecontinuityequation.

Next, we consider the second term of the momentum equation (15). We have assumed an equidistant grid in the

x_−space,andwe donot considerthe domainboundaries. Due to theseassumptions,andthefactthat theconsidered termis writ-teninaconservativeform,bothfinitedifferenceandfinitevolume methodsyieldthesamediscretization,whichisgivenby

∂

qu

∂

x

|

m+1/2≈

(

qu

)

m+1−

(

qu

)

m



x (22)

withm+1/2theindexindicatingthecell-centrepointinacontrol volume ofum+1/2. The faces ofthis control volume areindicated

withtheindicesmandm+1.Inthiscontext,thevelocityof mo-mentumum+1/2 isa cell-volume quantity,whereas the massflux qm+1/2 inEq.(20)isacell-facequantity.Bothqanduatthefaces

ofthecontrolvolumeneedtobeinterpolated.Sincethemassflux

qiscontinuous,thefollowinginterpolation qm≈ qm= 1 2



qm−1/2+qm+1/2



(23)

isaproperone. Here,theoverbardenotes thearitmethicaverage ofaquantityacrossthecellorcellinterface.Thisaveragingis sec-ondorderaccurateonauniformgrid.Forthetransportedvelocity

um(itmaynotbesmooth!)weapplyafirstorderupwindscheme,

asfollows um≈ ˆum=



_u m−1/2, ifqm>0 um+1/2, ifqm<0 (24)

Secondorder accuracycan be obtainedby meansof ahigher or-derupwindscheme,withorwithoutafluxlimiter.Thisschemeis treatedasadeferredcorrectiontothefirstorderupwindscheme, whilethis correction isadded to the righthandside ofthe mo-mentum equation [35]. In this way, second order accuracy is achieved, andthe sizeof the stencilassociated withthe first or-derupwindschemeisunalteredandremainslow.

The approximationofthethirdtermofEq.(15)isobtainedby meansofsecondordercentraldifferences,asfollows

h

∂ζ

∂

x

|

m+1/2≈ hm+1/2

ζ

m+1−

ζ

m



x with hm+1/2= 1 2

(

hm+hm+1

)

(25)

To demonstrate cellwise conservation of the momentum flux, wethusrestrictourselvestoastationary,positiveflow.Inthiscase

qisconstant,andsowehaveqm=qm−1/2=qm+1/2,andwithqm>

0andqm+1>0,oneobtainsfromEq.(22)

∂

qu

∂

x

|

m+1/2≈ qm+1uˆm+1− qmuˆm



x = qm+1/2um+1/2− qm−1/2um−1/2



x (26)

Weprovethecellwiseconservationofmomentumfluxasfollows qm+1/2um+1/2− qm−1/2um−1/2+ghm+1/2

(

ζ

m+1−

ζ

m

)

=0⇒ qm+1/2um+1/2− qm−1/2um−1/2+ghm+1/2

(

hm+1− hm

)

= ghm+1/2

(

dm+1− dm

)

⇒ qm+1/2um+1/2− qm−1/2um−1/2+ 1 2g



h2 m+1− h2m



= ghm+1/2

(

dm+1− dm

)

⇒



qu+1₂gh2



m+1/2= ghm+1/2[d]m+1/2 (27) Eq.(27) representsthe Rankine-Hugoniotrelation of momentum flux conservation applied at the staggered grid. Note that this relation holds on the interface xm+1/2, while the jump in qu is

taken upstream of cell edge xm+1/2 by virtue of flux continuity

(cf. Eq. (12)). It must be emphasized that this property remains validwhensecond orderaccuracybymeansofadeferred correc-tion scheme is included. It should also be noted that the water depthatthecellface,Eq.(25)asoccurredinthesurfacegradient,

istakenasthearithmeticaverage.Thisisrequiredtoobtaina well-balancedschemerepresentingthehydrostaticbalancebetweenthe pressurefluxandthetopographyterm[7,10]

1 2g

h2

m+1/2=ghm+1/2[d]m+1/2 (28)

that enables to preserve the quiescent water over varying bathymetry,i.e.

ζ

=constant,andu=0.

Continuity ofthemass fluxandmomentum flux acrossa dis-continuity, asgivenbyEqs. (19)and(27),respectively,requires a lossinthe(kinetic)energy(implyingan increaseofentropy)2 _[1].

Thislossis duetotheviscous dissipation andisprovidedby the upwindscheme (24).This numericalviscosity is alsorequired to dampoutoscillationsarisingintheregionsofsteepgradientsand numericaldispersionerrors.Althoughnumericalviscosityisan im-plicit and uncontrollable feature of the applied upwind scheme, the corresponding value of viscosity is solely determined by the cellwise conservation of mass and momentum fluxes across the discontinuity. Thisimplies acorrect magnitudeneededfor a vis-cositytohandlethediscontinuitiesintheprimitivevariables.This willbedemonstratedinSection6.3.1.

Theaboveapproximationsareoutlinedandproposedearlierin [24] and[27].They are also implemented in NEOWAVE [36] and SWASH[34].Inaddition,similardiscretizationequationswere pre-sentedbyZhou [37],althoughhesuggestedtoapply Eq.(25), in-steadofEq.(21),toapproximatethemassfluxfollowingEq.(20). Itmustbestressedthat,contrarytotheapproachesofDoyenand Gunawan[27]andZhou[37],thespatial discretizationspresented sofar havebeenobtainedfromequationsthat arenot writtenin strongconservationform, viz.Eqs.(14)and(15).Infact,they can evenbederivedfromshallowwaterequationsinnon-conservative form(see,e.g.[24]).

Wenowconstructastaggeredschemeforthesolutionof con-tinuity and momentum equations written in strong conservation form,viz.Eqs.(16)and(17).Themomentumvariableqis consid-eredastheunknownofthemomentumequation(17),whichisthe localmomentumperunitlength,whereasu=q/histheflowrate per unit width. Obviously,q istreatedasthe transportedstorage quantity,whereas u actsastransport velocity.Thispoint of view hasbeenadoptedby,amongothers,Stansby[18]andZijlemaand Stelling[38].Thecontinuityequation(16)isdiscretizedas

dhm

dt +

qm+1/2− qm−1/2



x =0 (29)

whichsatisfiesthejumpcondition(19)insteady-statecase. ForspatialdiscretizationoftheadvectiontermofEq.(17),first orderupwinding withrespectto qis employed, whiletakingthe arithmeticmeanofu,asfollows

∂

uq

∂

x

|

m+1/2≈

(

uq

)

m+1−

(

uq

)

m



x (30) with um≈ um= 1 2



um−1/2+um+1/2



(31) and qm≈ ˆqm=



_q m−1/2, ifum>0 qm+1/2, ifum<0 (32)

TheothertermsofEq.(17)arediscretizedcentrally,asfollows

∂

h2

∂

x

|

m+1/2≈

h2

m+1− h2m



x

2 As such, a flow across the discontinuity experiences a deceleration. That is to say, the combined approximations (19) and (27) can not be applied in re- gions of contraction with flow acceleration. For instance, Stelling and Duinmeijer

[24] demonstrated an erroneous increase in energy head for the subcritical flow over sudden contractions due to strict conservation of momentum.

and h

∂

d

∂

x

|

m+1/2≈ hm+1/2

dm+1− dm



x

Velocity um+1/2 remains to be defined, since qm+1/2 is the

re-solved variable.Stansby[18]suggested thefollowing ‘central’ ap-proach

um+1/2= qm+1/2 hm+1/2

(33)

However, as an alternative, this quantity may also be obtained fromEqs.(20)and(21),asfollows

um+1/2= qm+1/2 ˆ hm+1/2

(34)

which is referred to as the ‘upwind’ approach. Note that um+1/2

maynotbesmoothforflowwithspace-varyingdepth.

Wenowdemonstratethatthejumpconditioncannotbe guar-anteed. Let us consider a steady-state, positive flow. We have

qm+3/2=qm+1/2=qm−1/2,and

∂

uq

∂

x

|

m+1/2≈ um+1qˆm+1− umqˆm



x = =



um+1/2+um+3/2



qm+1/2−



um−1/2+um+1/2



qm−1/2 2



x = 1 2



um+3/2+um−1/2



qm+1/2− um−1/2qm−1/2



x (35)

Substitutioninthediscretizedmomentumequationyields



_u m+3/2+um−1/2 2



qm+1/2− um−1/2qm−1/2+ 1 2g



h2 m+1− h2m



=ghm+1/2

(

dm+1− dm

)

(36)

Thus,thejumpcondition(27) cannotberecoveredexactlyatthe staggered grid dueto an arithmetic averaging of flow velocity at

xm+1/2. Infact, um+1/2 isnot involved in Eq.(36),indicating that

thisschemeallowsodd-evendecoupling.Thiswilllead to inaccu-ratemomentumflux,aswillbedemonstratedinSection6.

The causeisthe reverserolesofthe transportedquantity and thetransport velocity.Inthisapproach,they areq andu, respec-tively, which should be the other way around, i.e. u and q, be-cause ofthe constraint

∂

q/

∂

x=0. Indeed, q is constant across a grid cell that is assured by the continuityequation (19).Next, if

∂

(

qu

)

/

∂

x=0, then u must be constant. Using the discretization accordingtoEq.(26)yields

um+1/2=um−1/2

implyingthat uisconstantoverthesamecell.However,with ap-proximation(35),oneobtains

um+3/2=um−1/2

displayingtheodd-evendecoupling,whichinducesacheckerboard pattern.

5. Timeintegration

Thestaggeredschemes,asoutlinedinSection4,areintegrated intime bymeans ofan explicitleapfrogschemeusingstaggering in time [31]. The flow velocity um+1/2 and the depth-integrated

velocity qm+1/2 are evaluatedat eachhalf time step

(

n+1/2

)

t,

whereasthesurfaceelevation

ζ

mandthewaterdepthhmateach

wholetimestepn



t,with



tthetimestepandn=0,...,N

indi-catingthetimeleveltn_=n

_

_{t with,intotal,N+1timesteps.}

We considerthreevariantsofthestaggeredschemes tobe in-tegratedintime.

VariantI

The firstvariant hasbeen discussed earlierinDoyen and Gu-nawan[27].Theyfirstsolvethecontinuityequationtoobtainhnm+1,

followedbythemomentumequationtoobtainu_mn+3₊₁/_/2₂.Hence,the continuityequation(29)isdiscretizedintime,asfollows

hn+1 m − hnm



t + qn_m+1₊₁/_/2₂− qn+1/2 m−1/2



x =0 (37)

Thisequationrepresentsthevolumebalance



x

hn+1 m − hnm

+



t

q_mn+1₊₁/_/2₂− qn+1/2 m−1/2

=0

demonstrating that the volume (per unit width) is conserved over time for an isolated physical domain. Furthermore, with

qn_m+1₊₁/_/2₂=_hˆn m+1/2u

n+1/2

m+1/2,onecanverifybyinsertionofEq.(21)into Eq.(37)thatthewaterdepthhnm+1remainsnon-negative[24],if

max



un_m+1₊₁/_/2₂,0



− min



un_m+1₋₁/_/2₂,0





t



x ≤ 1 (38)

Condition(38)isasufficientbutnotnecessarycondition.Byvirtue ofEq.(21),the waterdepth atthe cell interface, ˆhn_m+1_+1/2, is non-negativeunderthesamecondition. Thus, flooding neverhappens fasterthan onegrid cellper time step.Onthe other hand,when a wet cell becomesdry, the water depth can be arbitrarysmall, makingthe solutionofthemomentum equation meaningless. In-stead,themomentumequation(39)isnotsolvedandthevelocity

un_m+3₊₁/_/2₂issettozeroifthewaterdepthhˆn_m+1_+1/2isbelowathreshold value.Throughoutthiswork,thisvalueequals10−8 m,whichdoes notadverselyaffectthemassbalance.Apartfromthisminor mea-sure,nofurther modificationswithrespecttowetting anddrying arerequired.

Thetime-integratedversionofthemomentumequation,as pro-posedbyDoyenandGunawan[27],isgivenby

h_mn+1_+1/2un_m+3₊₁/_/2₂− hnm+1/2u n+1/2 m+1/2



t + q_mn+1₊₁/2uˆ_mn+1₊₁/2− qn+1/2 m uˆ n+1/2 m



x +1 2g

(

hn+1 m+1

)

2−

(

hnm+1

)

2



x =gh n+1 m+1/2 dm+1− dm



x (39) or,alternatively, hmn+1+1/2umn+3+1//22− h n m+1/2unm+1+1//22



t + q_mn+1₊₁/2uˆ_mn+1₊₁/2− qn+1/2 m uˆnm+1/2



x +ghnm+1+1/2

ζ

n+1 m+1−

ζ

mn+1



x =0 (40) with

ζ

n+1

m =hnm+1− dm. Notethat the third termis evaluated

us-ingtheupdatedhnm+1 toobtaina stablesolution[27],andthat no

systemsofequationsareinvolved.The secondtermofEq.(39)is evaluatedbymeansofEqs.(23)and(24).Noticetheroleofmass fluxqinEq.(39)beingthetransportvelocity.Assumingauniform bed,the momentumbalance canbe derived fromEq.(39), yield-ing



x

q_mn+3₊₁/_/2₂− qn+1/2 m+1/2

+



t



qn_m+1₊₁/2uˆn_m+1₊₁/2− qn+1/2 m uˆ n+1/2 m + 1 2g



(

hn+1 m+1

)

2−

(

hnm+1

)

2



=0

providedthatqn_m+1₊₁/_/2₂=hn_m+1/2un_m+1₊₁/_/2₂ istheconservedmomentum perunitlength.

Thetimestepforthisexplicitschemeislimitedby the follow-ingCFLcriterion



|

un_m+1₊₁/_/2₂

|

+



ghˆn m+1/2



t



x ≤ 1

Note that this CFL criterion determines the maximum time step whichsatisfiesautomaticallycondition(38).

VariantII

Thesecond variantresemblesthefirstoneinthesolution pro-cedure, exceptthat the momentum per unit length q is updated insteadofflowvelocityu,bysolvingthefollowingequation qn_m+3₊₁/_/2₂− qn+1/2 m+1/2



t + u_mn+1₊₁/2qˆ_mn+1₊₁/2− un+1/2 m qˆnm+1/2



x + 1₂g

(

h n+1 m+1

)

2−

(

hnm+1

)

2



x =gh n+1 m+1/2 dm+1− dm



x (41)

whichconservesqn_m+3₊₁/_/2₂ inthecaseofa uniformbedtopography. ThesecondtermofEq.(41)isevaluated usingEqs.(31) and(32). Since q is the resolved variable, the requirement of strict non-negativewater depths can not be properly enforced by Eq.(37). Inaddition,qis thetransportedquantity,whereasuisthe trans-port velocity which needs to be reconstructed, either based on Eqs.(33)or(34).Like thefirstvariant,anoflowconditionatthe cellinterface betweena dryandwetcellis imposed,so that the massbalanceismaintained.

Asimilarapproachinwhichthelocalmomentumqisupdated canbefoundinMadsenetal.[25].

VariantIII

SWASH isa wave-flow model intended to simulateboth long andshort wavesandrapidlyvaried flows [34]. Thismodelsolves the shallow water equations including the non-hydrostatic pres-surethattakesintoaccountthefrequencydispersionoftheshort wavesandvertical accelerations.For thisreason, andcontraryto thefirstvariant,inthisvariantthemomentumequationissolved firstandthenthecontinuity equation.We willdealwiththe fol-lowingexplicitmomentumequation(cf.Eq.(40))

hnm+1/2unm+1+1//22− h n−1 m+1/2unm−1+1//22



t + q_mn−1₊₁/2uˆ_mn−1₊₁/2− qn−1/2 m uˆnm−1/2



x + ghnm+1/2

ζ

n m+1−

ζ

mn



x =0 (42) Toobtainthefinaldiscretizedmomentum equationforthe primi-tivevariableum+1/2,westartwithanexplicitschemeforthe

con-tinuityequation(29),whichisgivenby hn m− hnm−1



t + ˆ hn−1 m+1/2u n−1/2 m+1/2− ˆh n−1 m−1/2u n−1/2 m−1/2



x =0 (43)

Next,we derive an equation by takingthe arithmetic average of Eq.(43)inpointsmandm₊1_,yielding

hn_m+1/2− hnm−1+1/2



t + qn_m∗₊₁−1/2− qn∗_−1/2 m



x =0 (44) with qn_m∗−1/2=1 2



ˆ hn−1 m−1/2u n−1/2 m−1/2+hˆ n−1 m+1/2u n−1/2 m+1/2



Multiplying Eq.(44) withu_mn−1₊₁/_/2₂ andsubtracting the resultfrom Eq.(42),afterwhichitisdividedbyhn_m+1/2,gives

un_m+1₊₁/_/2₂− un−1/2 m+1/2



t + 1 hn_m+1/2



_qn−1/2 m+1 uˆmn−1+1/2− qnm−1/2uˆ n−1/2 m



x − un−1/2 m+1/2 qn_m∗₊₁−1/2− qn∗_−1/2 m



x



+g

ζ

n m+1−

ζ

mn



x =0 (45)

After the momentum equation (45) is solved, the continuity equationisthensolved,whichisgivenby(cf.Eq.(43))

ζ

n+1 m −

ζ

mn



t + ˆ hn m+1/2u n+1/2 m+1/2− ˆhnm−1/2u n+1/2 m−1/2



x =0 (46)

Like the first variant, this schemepreserves non-negativity of thewaterdepthundertheCFLcondition(38)andremains stable, if



|

un_m−1₊₁/_/2₂

|

+



ghˆn m+1/2



t



x ≤ 1

In addition, assuming a uniform bathymetry, Eq. (45) conserves momentum perunit length hnm−1+1/2umn−1+1//22 over time (cf.Eq.(42)).

Notice that variants I and III are mathematically equivalent, al-thoughtheydifferintheorderofthesolutionwithrespecttotime stepping.

Concerning the temporal order of accuracy, all three variants are second order accurate in time in the absence of the advec-tiontermby virtueofstaggering intime. Toobtainsecond order accuracy in both space andtime, a predictor-corrector technique is adopted.In the predictorstep, the momentum equations(39), (45)or(41)issolved todetermine theprovisional flowvelocities

u∗_m+1/2 ormomentumvaluesq∗_m+1/2,respectively. Inthecorrector step, thesevalues are correctedby means of a second order up-wind schemeina deferred correctionfashion, whilemarching in timeusingtheMacCormackapproachtoachievesecondorder ac-curacyintime.See[34]forfurtherdetails.Theaccuracytestofthis approachwillbediscussedinSection6.2.

6. Numericalexperiments

6.1. Methodstobeassessed

Anumberofnumericalmethodsforthesolutionofthe depth-integratedshallowwaterequationshasbeendiscussedinthe pre-vioussections.Inthissection,wepresentthreestaggeredschemes, whichinturnwillbeassessed inlightofthenumericalaccuracy, stability,robustnessand conservationproperties.In essence,they aredistinguishedbytheuseofdependentvariableswhichare ei-therprimitiveorfluxvariables.

SchemeI-

Stagg

(

ζ

,u)

Thefirststaggeredschemeisthethirdvariantoftime integra-tionschemesasoutlinedinSection5andisgivenbyEqs.(45)and (46), augmented withEqs. (21), (23) and (24). This method em-ploys the primitive variables

ζ

andu. A key feature of this ap-proach is that the depth-averaged velocity u is transported by the massflux q, whichis computedwithEqs.(20)and (21). An-otherkeyaspectisthepreservationofnon-negativityofthewater depthbasedoncondition(38).Themethodsatisfiesthe Rankine-Hugoniotjumprelations (19)and(27),implyingcontinuityofthe massfluxandmomentumflux acrossdiscontinuities,respectively. Furthermore, the scheme is well-balanced, i.e. it preserves ex-actly the balance between the pressure flux and the bed slope term, so that artificial flows are excluded when quiescent still-water conditions are considered. This first scheme is denoted as

Stagg

(

ζ

,u)

. Note that this approach is implemented in, e.g. SWASH[34].Alsonotethat thisschemeissimilar totheone de-scribedin[27].

SchemeII-

Stagg

U

(h,q)

The second staggered scheme consists of Eqs. (37) and (41) combined with Eqs. (31) and (32). This is also the second variant of time integration schemes presented in Section 5. This scheme uses the water depth h and the depth-integrated veloc-ityq(ormomentumperunitlength)astheresolvedvariables.The depth-averagedvelocityuisobtainedbymeansofthe‘upwind’ ap-proach,i.e.Eq.(34).Wedenotethisschemeas Stagg

U

(h,q)

. Like thefirstscheme,thisschemepreservessteady-state flowsat rest.However, themethoddoesnotcomplywiththejump condi-tion (27), so that it allows to produce odd-even oscillations (viz.

Fig. 2. Root-mean-square error in water level ζ(in m; black squares) and in velocity u (in m/s; red dots) at t = 5 T as function of grid size obtained with different schemes. The numbers in the plots indicate the rate of convergence. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Eq.(36)).Inaddition,theschemecannotensurenon-negativityof thewaterdepth,sinceaconditionlikeEq.(38)cannotbededuced fromthisscheme.

SchemeIII-

Stagg

C

(h,q)

The third staggered scheme is similar to the second one, ex-cept the depth-averaged velocity u is estimated via the ‘cen-tral’ approach, i.e. Eq. (33). This scheme is designated as

Stagg

C

(h,q)

.

6.2. Convergencetest

To test theactual order ofaccuracy of thediscussed schemes we will employ the method ofmanufactured solutions [39]. We consider a progressivewave ina one-dimensional domain ofthe followingform

ζ

(

x,t

)

=a0cos



2

π



x L− t T





with a0,L andT the wave amplitude, length andperiod,

respec-tively.Thebedisuniform.Tosatisfy continuityequation(14),the mass flux is then given by q

(

x_,t

)

₌c0

ζ

(

x,t

)

with c0=



gd the wave speed. Next, thismanufactured solution is substituted into momentumequations(15)or(17)yielding non-zeroterms,which can be added to the momentum equation. These terms act as forcing termsforpreventingtheprogressivewavefrombeing de-formed by nonlinearity. Also, thegeneration of higherharmonics isprecluded.Inordertotestalltermsofthediscretizedequations, including momentumadvection,the amplitudeofthe chosen so-lutionwasset toa0=0.5d.Theboundary conditionatx=0was q

(

t

)

=a0c0cos

(

2

π

t/T

)

andazerovelocitywasappliedatx=10L.

Thisboundarylocationissufficientlymovedawayfromthedomain ofinterest soasto not tointroduce reflectionerrors inthe solu-tion.Both

ζ

anduwereinitiallysettozero.

A grid refinement study wasconducted to assess the conver-gencepropertiesofthethreediscussedschemes.Forthispurpose second order central differences were implemented for the ap-proximation of momentum advection by means of the deferred correction approach [35]. Furthermore, for the computation of

hm+1/2 the aritmethic averagewasemployed,i.e.hm+1/2≈ hm+1/2,

viz. Eq. (25), instead of Eq. (21). As a consequence, schemes

Stagg

U

(h,q)

and Stagg

C

(h,q)

are identical. Computa-tional grids consisting of 10, 20, 40 and 80 grid cells per wave length were used.Inorder tokeep theCourantnumberc0



t/



x

constant, thetimestep



twasreducedforeach gridrefinement. ACourantnumberof0.2waschosen,sothat50,100,200and400 timestepsperwaveperiodwereemployed,respectively.

The root-mean-square errors for

ζ

(x,t) and u(x,t) were mon-itored after five wave periods in a part of the domain, namely 0≤ x≤ 5L. The results are displayed in Fig. 2. The observed er-ror is a combination of spatial and temporal contributions ob-tainedondifferentgridsandtimesteps.Giventheuseofrelatively coarse grids, the rate of convergence of schemes

Stagg

(

ζ

,

u)

and Stagg

C

(h,q)

isinclinedtosecondorderwithrespectto



x(andto



t)forboth

ζ

andu.

6.3.Idealizedtestcases

Inthis section we present three idealizedtest caseswith the objectiveofassessingthecapabilitiesandlimitationsofeachofthe methodsdescribedinSection6.1.Inparticular,wewilladdressthe following:

• the continuity of mass and momentum flux across a bottom step;

• theavoidanceofunwantedeffectsinducedbythecheckerboard modes;and

• the accuracy with which either the steady state or unsteady flowsolutionsareresolved.

Thethreetestcasesinclude(i)expansioninopenchannelflow; (ii)transitionfromsupercriticaltosubcriticalflowacrossabottom step;and(iii)unsteadydambreakflowondrybed.

6.3.1. Subcriticalflowoverbackwardfacingstep

Weconsiderafrictionlesschannelwithaconstantwidthanda jumpinthebottomprofileasdepictedinFig.3.Thelengthofthe channelis10km.Upstreamofthisjumpthebedlevelis5mbelow datumlevel(d1=5m),anddownstreamthebedlevelis10m

be-lowdatumlevel(d2=10m).Hence,theheightofthebottomstep

is5m. In addition, upstream ofthe jump a discharge of10m2_/s

isgiven,whereasdownstreamawaterlevelof0misimposed.The maximumFroudenumberforthistestis0.3,sotheflowis subcrit-icaleverywhere. Based on elementary balanceprinciples inopen channel flow [40,41], theresulted flow experiences an expansion downstreamofthejumpwithalossintheenergyhead.Bothmass andmomentumconservationmustholdeverywhere.Thisexample isanidealizedstudyofrapidlyvariedsteadyflowwithanabrupt changeinthewaterdepthwithwhichthedepthandmomentum oneithersideofthejumpisgiven.

Acrossthejumpwehavethefollowingequations u1h1=u2h2=10m2s−1 u2 1h1+ 1 2g

(

h1+5

)

2_=u2 2h2+ 1 2gh 2 2

withindices1and2indicatingthelocationsupstreamand down-streamofthejump,respectively.Notethatthelefthandsideofthe momentum balanceincludes the hydrostaticpressure against the wall (of thestep). Since the downstream waterdepth is fixed in thetest,wewill examinetheupstreammomentum fluxobtained withdifferentmethods.Bymeansofarootfindingmethodits the-oreticalvalueisobtainedas

u1q1+ 1 2gh

2

Fig. 3. Sketch of expansion flow over a bottom step.

Fig. 4. Computed net mass flux (upper panel) and momentum flux (lower panel) along the channel for the expansion flow with a jump in the bed. The location of this jump is x = 0.

Tocheck thecontinuityofmomentumfluxacross thejumpat thediscretelevel,weconsiderthefollowingadapteddownstream momentumflux u2q2+ 1 2gh 2 2− 5 2g

(

h1+h2

)

withthelasttermrepresentingthetopographyterm.Notethatfor eachmethodthedownstreammomentumfluxwascomputedwith theindices 1and2 corresponding tothe consecutivegrid points wheretheactualjumpoccurredinbetween.

The meshwidthwassetto10mforall methods,whereas the time step was0.5 s. This time step was fixed a priori, andthen leftunchangedthroughoutsimulations.ItfulfilledtherequiredCFL conditionassociatedwiththeusedschemes(maximumCFLis0.8). Initially, both surface elevation and velocity were zero ev-erywhere. The steady state results of the methods are de-picted in Fig. 4. The first panel illustrates mass conservation: the net mass flux must be zero. The schemes

Stagg

(

ζ

,

u)

and Stagg

U

(h,q)

appeared to be strictly mass conservative, whereasthe Stagg

C

(h,q)

methodismassconservativeupto machineprecision. Thesecond panel demonstratesthe continuity ofmomentum flux across thejump. All methods clearly showed perfectcontinuity of the momentum flux, implying that they all satisfyjumpcondition(27).Thiscanbeexplainedbythefactthat

thegradientofthesolutioniszeroeverywhereexceptatthe bot-tom step,so that the set offinite difference equations automati-cally reducesto the algebraicexpressions (19)and(27).Schemes

Stagg

U

(h,q)

and Stagg

C

(h,q)

underestimated the up-streammomentumflux,whereasthisfluxascomputedby Stagg

(

ζ

,u)

isprecisely137.909m3_s−2_.

Although all discussed schemes comply with the Rankine-Hugoniot relations in the jump region, still the overall numeri-cal solution may not be correct. As pointed out by, e.g. Stelling andDuinmeijer [24], Bernetti etal. [42] andMurillo and Garcia-Navarro [43], the numerical evaluation of the hydrostatic force exerted by the bottom step on the fluid is responsible for this. For instance, Bernetti et al. [42] derive the jump condi-tions for mass and momentum whereby the total energy across the discontinuity must be dissipated correctly (i.e. the entropy condition is fulfilled). In this way, non-physical solutions are excluded.

Letusnowreturntotheexpansionflowtest.Wecalculatethe downstreammomentumfluxusingthejumpconditionsofBernetti etal.[42],asfollows u2q2+ 1 2gh 2 2− 1 2g[h1+

(

ζ

1+d2

)

]

(

d2− d1

)

=137.909m 3_s−2