An implicit wetting and drying approach for non-hydrostatic baroclinic flows in high aspect ratio domains

An implicit wetting and drying approach for non-hydrostatic baroclinic flows in high aspect

ratio domains

Candy, A. S.

DOI

10.1016/j.advwatres.2017.02.004

Publication date

2017

Document Version

Final published version

Published in

Advances in Water Resources

Citation (APA)

Candy, A. S. (2017). An implicit wetting and drying approach for non-hydrostatic baroclinic flows in high

aspect ratio domains. Advances in Water Resources, 102, 188-205.

https://doi.org/10.1016/j.advwatres.2017.02.004

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

Advances

in

Water

Resources

journalhomepage:www.elsevier.com/locate/advwatres

An

implicit

wetting

and

drying

approach

for

non-hydrostatic

baroclinic

flows

in

high

aspect

ratio

domains

A.S.

Candy

∗

Department of Earth Science and Engineering, Imperial College London, UK

Environmental Fluid Mechanics Section, Faculty of Civil Engineering and Geosciences, Delft University of Technology, The Netherlands

a

r

t

i

c

l

e

i

n

f

o

Article history: Received 20 April 2016 Revised 1 February 2017 Accepted 5 February 2017 Available online 6 February 2017

Keywords:

Wetting and drying Non-hydrostatic Baroclinic

High aspect ratio domains Multi-scale simulation Vertical inertia Finite element method

a

b

s

t

r

a

c

t

Anewapproachtomodellingfreesurfaceflowsisdevelopedthatenables,forthefirsttime,3D consis-tentnon-hydrostaticbaroclinicphysicsthatwetsanddrysinthelargeaspectratiospatialdomainsthat characterisegeophysicalsystems. Thisis keyinthe integrationofphysicalmodels topermitseamless simulationinasingleconsistentarbitrarilyunstructuredmultiscaleandmulti-physicsdynamicalmodel. AhighordercontinuumrepresentationisachievedthroughageneralGalerkinfiniteelementformulation that guaranteeslocal and global massconservation, and consistenttracer advection.A flexiblespatial discretisationpermitsconformingdomainboundsandavariablespatialresolution,whilstatypicaluseof fullyimplicittimeintegrationensurescomputationalefficiency.Notablythisbringsthenaturalinclusion ofnon-hydrostaticbaroclinicphysicsandaconsiderationofverticalinertiatofloodmodellinginthefull 3Ddomain.Thishasapplicationinimprovingmodellingofinundationprocessesingeophysicaldomains, wheredynamicsproceedsoveralargerangeofhorizontalextentsrelativetoverticalresolution,suchas intheevolutionofatsunami,orinurbanenvironmentscontainingcomplexgeometricstructuresata rangeofscales.

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

1. Introduction

Flooding has huge impacts on the economy of a region and the livelihood of its people. Significant progress has been made tomodelandpredict theimpactoftheseinundationevents, cap-turingthecharacteroftheirsourceandresultantbehaviour.Many challengesstill existand in particular inconcurrently simulating thephysicalprocessesinvolvedfromthelargeplanet-scaleforcings downtothesmallhumanscalesofanurbanenvironment.Thisis highlightedinthereview(MedeirosandHagen,2013)asoneofthe keylimitationsofexistingwettinganddrying(WD)models.Inan urbanfloodingscenarioforexample,modelledwatercolumndepth couldbedownto1cmoverahorizontalrangeoftensorhundreds ofkilometres,leadingtoaveryhighaspectratioof∼ 10−7_.

Inundation flow models typically use simplified formulations oftheNavier–Stokesequations,commonlytheSaint–Venant shal-lowwater equations(SWEs).Thesesimplifications make assump-tions,suchasahydrostaticpressureandwell-mixedwatercolumn,

∗ _{Corresponding author at: Environmental Fluid Mechanics Section, Faculty of}

Civil Engineering and Geosciences, Delft University of Technology, Stevinweg 1, 2628 CN, Delft, The Netherlands.

E-mail address: a.s.candy@tudelft.nl

whichare notnecessarilyvalidinthewhole rangeofscales rele-vantto theinundation. Non-hydrostaticprocesses become impor-tant, forexample, inthe dispersive effectsof short waveswhere the ratioof vertical andhorizontal scales of motionare not suf-ficiently small. The study of Oishi et al. (2013) considering the 2011T¯ohokutsunamiinJapan,founditwascriticaltoinclude non-hydrostaticeffectstocorrectlymodelprocessesonthesmallscale, apointfurtherhighlightedbyCuietal.(2014).

Typically, multi-physics over a broad range of scales is ap-proachedusing multiplemodelrunsata hierarchyofscales such that domains are nested, with varying complexity and included physics. As an alternative, efforts to integrate the physics and scalesofseparatemodelsintosingleEarthsystemmodelsis grow-ing, where it is important individual components function in a generalcontext,andarenottoorestrictiveindiscretisationchoice. Althoughthiscanbeachievedweaklywithofflinecommunication betweenmodels,the‘holygrail’isaflexible singlemodelcapable of simulating a range of physics andscales, with inherit consis-tencyandconservation.

This work pushes the boundaries in two key regards: Firstly, addinganovelapproachtoWDinthe‘thin-film’ familysolving a full3DpressureratherthantheusualSWEapproximationinvery challengingacutelylargeaspect ratiodomains typicalof

geophys-http://dx.doi.org/10.1016/j.advwatres.2017.02.004

icalsystems— a firstforWD.Secondly,thisbringsthemodelling ofWDprocessestogetherwithnon-hydrostaticbaroclinicflow dy-namicsinasinglesimultaneousandseamlesssystemmodel.This is critical for tightly coupled processes, for example in tracking groundinglinemovement underan iceshelf oceancavity, thatis stronglyinfluencedbynon-hydrostaticandbaroclinicoceanflows. Accuratelytracking aninundationinterface istechnically chal-lenging. Of the Eulerian type WD approaches (reviewed in

MedeirosandHagen,2013),wheretheunderlyingspatial discreti-sationispredominantlyindependent ofspaceandtimesuch that matrix operators can be cached and there is no need for com-plex contourtracking,therearefourtypes:elementremoval, lim-itingthecomputationaldomaintothewetregion(seeCasulliand Stelling, 1998, UnTRIM Casulli and Walters, 2000; D’Alpaos and Defina, 2007; Defina, 2000 and the WASH123D code Lin et al., 2004); thin film approaches (see Bates and Anderson, 1993, the FVCOM modelChen etal., 2003,POM model Oey, 2005andalso

Begnudelli and Sanders, 2006); depth extrapolation from wet to dry cells (Bradford and Sanders, 2002; Lynett et al., 2002); and negativedepth(Henicheetal.,2000;JiangandWai,2005)applied inROMS(Warneretal.,2013),includingtheuseofporousmedia belowtheseabed(vantHofandVollebregt,2005;Ipetal.,1998) andbathymetrymovement(Kärnä etal.,2011).

Underlyingmodeldiscretisationslargelysteerthischoice, with thefirstby farthemostcommonforexplicittime stepping mod-els, applied in QUODDY,ADCIRC, MIKE21, Delft3Dand in one of the first Eulerian methods (Leendertse, 1970), subsequently re-viewed in Balzano (1998). Whilst robust, stabilityconstraints re-strict movement of the interface to one cell per time step (



t), since the Courant number must be maintained lessthan one in drying regions (Stelling and Duinmeijer, 2003) to ensure a non-negativeboundonwaterdepth,astrictlimitationon



t(Walters, 2005). Depth extrapolation also suffers this restriction with ele-ments switchingstates(MedeirosandHagen,2013),whereas thin filmandnegativedepthoptionscanbetime-integratedimplicitly.

Forspatialdiscretisations,WDprocedureswerefirstappliedto structured meshes(Casulliand Stelling, 1998; Stelling and Duin-meijer, 2003), with updates to include non-hydrostatic correc-tions(StellingandZijlema,2003),baroclinicsolvers(Warneretal., 2013)andrecentlysubgridinformation(CasulliandStelling,2010; D’Alpaos andDefina,2007;Defina,2000;Volp etal.,2016)to in-cludehigherresolutionbathymetryandfluxcalculations.

Current approaches to unstructured mesh geophysical fluid modelling is considered in detail in Danilov (2013), with its po-tential importance best highlighted in Danilov et al. (2013). In-deed,thisreviewstatesthatwhilstunstructuredmeshmodelsmay notreplacestructuredmodellingapproachescompletely,thereare cases where this type of approach could be optimal. In particu-lar,allowingaflexibleapproachtotheverticaldiscretisationcould improveaccuracyandmodelefficiencyindomainswherethereare sharpchanges inbathymetry relativeto horizontalspatial resolu-tion, strong non-hydrostaticgradients inpressure, strong vertical inertialflows,orwhenitwouldbemoreoptimaltoreduceor in-creasethenumberoflayers inshallowanddeepregions, respec-tively.Moreover,thesecouldbecriticalinthefringesoftheocean boundary, along geometrically complex coastlines andin interac-tions withother types ofphysicalsystems, such asan urban en-vironment, orthecomplex shallowing iniceshelf ocean cavities. Within this discretisation type, WD models can more accurately modelawiderrangeofscalesinlargerdomains.

Oneoftheearlyfinitevolume(FV)approachesUnTRIM(Casulli and Walters, 2000) permits unstructured meshes with the con-straint that, like structured models,the domain elementsare or-thogonalwherecircumcentresareinsidetheirrespectiveelements. Its non-hydrostatic advance (Casulli andZanolli,2002) is applied intheSUNTANSmodel,withthesameorthogonalityrestriction.It

containsaWDalgorithm(Wangetal.,2009)stabilisedwitha tech-niquefromIpetal.(1998)thatappliesanincreaseddragtosatisfy anadditionalconstraintonvolumefluesindryregions.Similar ap-proachesarealsomadeinCuietal.(2010,FVCOM,MIKE21)with non-hydrostaticcorrectionsadded(e.g.Cuietal., 2012).FVislow orderonlyandmodelsgenerallyexplicit.

Unstructuredfiniteelement(FE)methodsofferhighorder con-tinuumapproximationswhicharemoreaccurateandnaturally in-clude less diffusive and dispersive advection schemes. WD has beenbuiltinto2D barotropicflow modelssuch asQUODDYwith dryelement removal in Greenberget al. (2005); ADCIRC, a SWE methodforexplicithydrostaticmodellingofstormsurgeswithdry removal(Dietrich et al., 2006);TELEMAC, initially usingelement removal(BatesandHervouet,1999) andnownegativedepth;and SLIM(Kärnä etal.,2011)witharepositionedseabedSWEmethod andadoptionofimplicit



tadvance.

WD is combinedwith solvers capable of modelling baroclinic processes in Funke et al. (2011); Warner et al. (2013), with the formerusingthinfilm highorderFEandthe latterexplicitfinite difference with negative depth WD and mode splitting.The for-merperforms well inrelatively modestaspectratio domains,but performanceisstrictlylimitedbyuseofdirectsolvers(LU decom-position),restrictions ondryelementaspect ratiosanderroneous unphysicalflowsthatdevelopindryregions.

HereageneralapproachforWDwithFEsisconsideredinfull 3D, building on established methods for modellingfluid flow on fully3D unstructuredmeshes(Piggottetal., 2008)which varyin resolutionandsupportamultiscaleofphysicalprocesses,including non-hydrostaticandbaroclinic dynamicsin thelarge aspect ratio domainsfoundingeophysicaldomains.Undertheconstraintsofa globalnumberofdegreesoffreedom,thisallowsthefocusof com-putationalresources on small scaleregions and areasof interest, whilstcapturingthelargescaleflowselsewhereinthedomain.An additionaladvantage isnoconstrainton theinternal mesh struc-ture,northatitisfixedintime.Itisnotconstrainedtolayers,and can be completely, orpartially in selectregions, fully anisotropi-callyunstructured.

Toallowefficienttimeintegrationoverarangeofelementsizes, animplicittreatmentnecessitatesa continuumapproachto inter-facetracking.Athinfilmisapplied,whichas(MedeirosandHagen, 2013) notes,generallysatisfiesmassandmomentumconservation withoutsignificantspecialtreatment,andproducesarealisticand smoothwettingfront.WDisincludedinanaturalmanner,through additionaltermsinthemomentumequationandmodified bound-aryconditions. Indeed,the numerical treatment is carefulto en-surethesolutionremainsintheSobolevsolutionspaceofthe orig-inal physically-based weakly formulated Galerkin problem. Prog-nosticvariables, including tracers, are self-consistent through the FEformulation andnotably, through useof a combined pressure variable,consistencywiththefreesurfaceisnaturallyinherent.

InthefollowingSections2–4,thenewconsistentapproachfor WDin large aspect ratiogeophysical domains is developed,with details of mesh movement inSection 5 andadditional strategies notedinSection6.ThisisvalidatedinSection7andevaluatedin

Section8.

2. Governingcontinuumequations

2.1. 3DBoussinesqwithpiezometricpressure

The non-hydrostaticBoussinesq equations fora rotating strat-ified fluid, are solved in a time-dependent domain



⊂ R3 _(see Fig.1),boundedbythesurface



.Thisissplitintothefreesurface boundary



_η,andtheremainingbound



b=



\



η.Theseare

de-fined forthe prognostic variables of velocity u:



× [0,T

)

→R3_, andpressurep:



_{× [0,}T

)

_→R,overthetimeinterval[0,T),such

Fig. 1. A schematic of an example high aspect ratio geophysical inundation domain considered here, with a ‘horizontal’ length scale L spanning its extent on a geoid surface, and the ‘vertical’ length scale H . In reality, these length scales differ by many orders of magnitude.

that

ρ0



∂

_∂

u_t +u·

∇

u



−

∇

·

μ∇

u+

∇

p=−g

ρ

ng, (1)

∇

· u=0, (2)

where

μ

is thetensorial dynamicviscosity, −ng andg the

gravi-tationalaccelerationdirectionandmagnituderespectively,and

ρ

:



× [0,T

)

→Rthedensity.Thelatterissplitintoabackground

ρ

0, andperturbationdensity

ρ

,such that

ρ

=

ρ

0+

ρ

.Since the hy-drostaticcomponentofpressureoftheequilibriumstate doesnot havean importantcontributiondynamically, itissubtracted from themomentum equationandthefullpressure p,isreplacedbya

piezometricpressure,commonlyappliedincoastalengineering ap-plications(e.g.LabeurandPietrzak,2005),definedas

p:=p+

ρ

0 gng· r+pa, (3)

for a position vector r, relative to a position where hydrostatic pressureiszero.Atmospheric pressureattheinterfaceisdenoted

pa.

Redefining the prognostic pressure with thisparticular choice ofpiezometricpressureformsa combinedfree surface– pressure prognostic p

(

p,

η

)

eliminatingthe needtosolve aseparate, com-monlyused,waveequationforthefreesurface,denotedbythe in-jectivefunction

η

:



× [0,T)→



η.Theprognosticpressurepnow containsnon-hydrostaticcomponentsandthehydrostaticpressure dueto perturbationsin thefreesurfaceelevation. Thisremaining hydrostaticpressure

ρ

0g

η

,is theboundaryconditionforpat



η, andthrough(3),wefind

p



_η=

ρ0

g

η

. (4)

2.2.Boundaryconditions

With the inclusion of the free surface height in the prognos-tic pressure, the kinematic free surface boundary condition of

AppendixAisexpressed n· ng

∂

p

∂

t





_η=

ρ

0gn· u, on



η. (5)

Thisistheboundaryconditionfor

η

andnowarequiredconstraint forthecombined p

(

p_,

η

)

prognosticvariable.Thisisjoinedbythe uandpboundaryconstraints

u· n=0,

∀

x∈



b, and

p=pa,

∀

x∈



η. (6)

Moregeneralconditions,foropenoceanboundariesorfluxinputs, canbeappliedwithoutfundamentalchangestotheapproach.

2.3. Coordinatesystemandframeofreference

Note additionally that the directionof gravity, describing the normalng,isnotrestrictedtoaCartesian z-component,such that

thedevelopmentisrelativelyindependentofthecoordinate refer-enceframe.ItisfreetovaryarbitrarilywithinR3_{,alignedwiththe} local directionof gravitational acceleration,and it ispossible for example,to apply thismethodto thespheroid shell oftheEarth inaCartesiancoordinatereferenceframe.

3. Spatialandtemporaldiscretisation

The non-linear system of Eqs. (1) and (2), combined with boundaryconditions(5)and(6),aresolvedforp

(

p,

η

)

,and veloc-ityu,usinga Chorinprojection method(Chorin,1967) toenforce incompressibility.Thisisa modifiedpredictor– corrector scheme based on Gresho et al. (1984) in which a predictor un+1

∗ is

ob-tained frommomentum conservationthat isnot divergence free, suchthatacorrectionun+1_=un+1

∗ −

∇

φ

isthencalculatedsubject tothedivergence-freeconstraint

∇

· un+1_{=0.Foreach timestep,} this proceeds for a number of Picard iterations until sufficiently converged.

3.1. Temporaldiscretisation

Discretisation in time is achieved by the

θ

-method (Iserles, 2012)inallcases,foratimestep



t,suchthattheexplicitforward Euler, Crank-Nicolsonand backward Eulertime-stepping schemes canbe obtainedwithchoicesof

θ

=0, 1

2 and1,respectively. The modified Navier-Stokes withimplicit free surface system (1)and

(2)atatimenistherefore

ρ0

un+1

_

− u_t n =Rn+θ₋

_∇

_pn+θ₋

_ρ

n+θ_{gng, (7)}

∇

· un+1_{=0, (8)}

whereRn+θ₌

_θ

_Rn+1₊

₍

₁₋

_θ

₎

_Rn _{containsthe advectivemass flux}

term, together with viscosity and other source terms. A choice

θ

∈[1

2,1] leadsto an implicittime stepping scheme that allows simulations to use large time steps, which are not restricted by theCourant-Friedrichs-Lewy(CFL)condition(Courantetal., 1928) with respect to the velocity and wave speed. In practice for the simulationspresentedhere,fortherequiredlevelofaccuracyand stability,Courantnumbersupto10areapplied.

3.2. Combinedfreesurface– pressureChorincorrector

UnderaGalerkinFEspatialdiscretisationthetemporally discre-tisedmomentum(7)andcontinuity(8)equationsare testedwith the velocity

φ

i and pressure

ψ

i basis functions, respectively. The

trialfunctionsuandparedefinedintermsoftheirrespective ba-sis functionsalso, and Appendix B describes their formand the nomenclatureusedinmoredetail.Thisleadstothespace-time dis-cretemomentumequation

ρ0

Mu



t



un+1 ∗ − un



+

θ

A˜n+1_un+1 ∗ +

(

1−

θ

)

Anun =

θ

C pn+1 ∗ +

(

1−

θ

)

C pn+Su, (9)

ρ0

Mu



t



un+1− un



₊

_θ

_An+1_un+1₊

₍

₁₋

_θ

₎

_An_un =

θ

C pn+1₊

₍

₁₋

_θ

₎

_{C p}n_+S u, (10)

fora Picarditeration stepandendoftime step, respectively.The starred variable un+1

∗ is the current best approximation to un+1, calculated from pressure at the previous time level n. The best guess of the solenoidal velocity at a time level n+1 is denoted

˜

un+1,andusedinthecalculationofupdatednon-linearoperators, such asmassfluxA˜n+1_{.Thevelocity spacemassmatrixMu} addi-tionallycontainsthediagonalorblock-diagonal(dependingonthe chosen discretisation)componentofviscosityfromR, whichisto betreatedimplicitlyinpressure.Thediscretecross-spacegradient operatorC_{i j}:₌₋_

φ

_i

∇

ψ

_jd



_,containsaninnerproductover ve-locityandpressurespaces,leavingsourcesinSu.

Subtracting(9)from(10)andmultiplying by

θ



tCT_M−1 u yields

adiscretePoissonequationforthecorrection

ρ

0

θ

CT

(

un+1− un_∗+1

)

=

θ

2



tCTM−1u C

(

pn+1− pn∗+1

)

. (11) 3.3. Discretecontinuity

DiscretisationofthecontinuityEq.(8)iswritten

θ

CTun+1+

(

1−

θ

)

CTun+GT_θun+1+G₁T_−θun=0, (12) where G_{θ,i j}:₌_

η

θ

n

φ

i

ψ

jd



andG(1−θ ),i j:=η

(

1−

θ

)

n

φ

i

ψ

jd



.

ForG_θ=G_{(1−θ )=}0_,thesystemofequationsenforces incompress-ibilitywithweaklyappliednonormalflowboundaryconditions.

3.4. Discretemodifiedkinematicboundarycondition

Discretisationofthefreesurfaceboundarycondition(5)isnow requiredtoprovidethe boundaryintegraltermsin(12),witha

θ

timediscretisationdescribedby

nn+1· ngpn+1− nn· ngpn

=

ρ0

g



t



θ

nn+1_{· u}n+1₊

₍

₁₋

_θ

₎

_nn_{· u}n



_{. (13)}

Discretisationof(13)inspaceusingthetestandtrialfunctions

φ

i

and

ψ

i,introducedinSection3.2gives Ms

pn+1_{− p}n

ρ

0 g



t =

GT

θun+1+GT1−θun, (14) withthesurfaceintegralMs,i j:=ηng

ψ

i

ψ

jd



.

Applying this discrete combined p

(

p_,

η

)

kinematic condition

(14)tothediscretecontinuity(12),wefind

θ

CT_un+1₊

₍

₁₋

_θ

₎

_CT_un_+M s

pn+1_{− p}n

ρ

0 g



t =

0. (15)

Substituting(15)intothemomentumEq.(11),withthecorrection defined



p:₌pn+1_{− p}n+1 ∗ ,yields



θ

2_CT_M−1 u C+ Ms g

(

t

)

2



p =−

θ

CTun∗+1+

(

1−

θ

)

CTun



t − Ms g

(

t

)

2

(

p n+1 ∗ − pn

)

. (16)

3.5.Combinedfreesurface– pressuresystemsolution

DuringasinglePicarditeration,thefirstvelocitypredictorstep solves the discrete linearised momentum Eq. (9),to establish an updatedintermediatevelocityun+1

∗ ,fromthecurrentbest approx-imationtothevelocityandpressure,andtheirvalueatthe previ-oustimestep.

The predictorun+1

∗ obtained isnot divergence free ingeneral, andinordertoenforcetheincompressibilitycondition, apressure correction



p is calculated to project this velocity into the diver-gencefree subspace by solving (16)above. Thevelocity correction

is madeto update the intermediate velocity, consistent with the new intermediate pressure, and projected to the divergence free subspaceusingthedifferenceof(9)and(10),where

un+1_=un+1

∗ +

θ



t

ρ0

M−1u C



p. (17)

Finally,theinterfacetrackingstepadjuststhefreesurface posi-tionfollowing(4)inlightofthenewpressurefield,inadirection −ng,paralleltothegravitationalvector.

4. Highaspectratiowettinganddryingdomains

4.1. Wettinganddryingofthesimulationdomain

Thefreesurfaceboundaryissplitintodistinctwetanddry re-gions(illustratedinFig.1),definedbythecombinedp

(

p,

η

)

atthe surfacesuchthat



_η=



w∪



d,with



w:r∈



η,

∀

p

(

r

)

≥

ρ0

g

(

h

(

r

)

+d0

)

, and



d:r∈



η,

∀

p

(

r

)

<

ρ

0g

(

h

(

r

)

+d0

)

.

The conditions in dry regions differ from those in wet in two definingways.Firstly, thewatercolumndepthismaintained ata thresholdminimum d0 above the bottom bathymetry definedby h:



→R, and secondly,the surface boundary condition on the combined p

(

p,

η

)

prognostic variable is modified to enforce this constraintinthesolver.Withthedepthconstraintthefreesurface evolutiondescribedby(4)providesthefirstconstraint

η

(

r

)

=max



1

ρ0

gp

(

r

)

,h

(

r

)

+d0

, forron



_η. (18) Thesecondisfoundby modifyingthecombinedkinematic condi-tion(5),whichinlightofthedepthrestrictiongives

n· ng

∂

∂

tmax

(

p,

ρ0

g

(

h+d0

)

)

=

ρ0

gn· u, on



η. (19) Inwetregions



w,theconstraints(18)and(19)reducebackto

thefreesurfaceconditions(4)and(5),respectively.Indryregions



w,(19)isa nonormalflow condition,andeffectivelyimposes a

rigidlidapproximation.

4.2.Conditioningofthepressurecalculation

Thespatialdomainsofgeophysicalprocessesaretypicallylarge aspectratio,duetothegravitationalinfluenceanddisparityin dy-namicsparallel andperpendicular to geoid surfaces.The solution ofanon-linearfluidflowsysteminthesetypesofdomains includ-ingnon-hydrostatic dynamics,withimplicittime evolutionanda predictor – corrector approach such as Section 3.5 is shown in

Krameretal.(2010) toleadtoan ill-conditionedpressuresystem. Inthelimitoflarge domainaspectratioandlongtime steps,the systembehavesapproximatelyasthough ithasa rigidlid,where thefreesurfaceisfixedwith

η

=0andu· n=0.

ThedryregionsintroducedbytheWDprocesssignificant exac-erbateill-conditioning,sincearigidlidconditionisapplieddirectly andtheregioncontainselementswithacutelylargeaspect ratios duetotheirdefiningshallowwatercolumndepth.

Thecorrectionun+1_=u∗−

_∇

_φ

_{(17)iscalculatedsubjecttothe} divergence-freeconstraint

∇

· un+1_{=0.Thisleadstothefollowing} pressurePoissonequationfor

φ

∇

2

φ

=

∇

· u∗, (20)

whichcorrespondstothediscretePoissonoperatorCT_M−1 u C inthe

formulations (16) and (42) above. For no normal flow boundary conditions where u· n=0, at interfaces with bedrock or in the caseoftherigidlidapproximationfortheocean-airinterface,the couplingbetweenvelocityandpressureresultsinthe correspond-ingboundaryconditionon(20)astheNeumannexpression

∂φ

∂

n =0, on



η, (21)

ensuringthevelocityconstraintisconsistentlypreserved.

Applying a kinematic condition instead leads to the homoge-neousDirichlet condition

φ

=0 on



_η. The redefinition of pres-surein(3)toformthepiezometricpressurehereallowsstandard pressuresplittingapproachestotreatbaroclinicandbarotropic dy-namics (e.g. Shchepetkin and McWilliams, 2005) in the general caseofdomainsdiscretisedwithfully-unstructuredmeshes. These schemesthemselvesaidtheconditioningofpressuresolvesin geo-physicalmodels(Maddisonetal.,2011),wherethereisalarge dis-parityofscalesandresolutionofthedominantphysicalprocesses. ThispiezometricvariablesatisfiesthesameEq.(20),witha modi-fiedright-handsidesourcetermandboundarycondition.

Discretisation of the kinematiccondition (5)defined interms ofthepiezometricpressureusingimplicitbackward Eulerintime givesaRobinconditionfor

φ

,suchthat

n· ng

φ



t2 =g

∂φ

∂

n. (22)

Withthe barotropic wave speed c=

gH, for a distanceH, and notingthat

φ



t2



g

∂φ

∂

n ≈ H g



t2 =



_H c



t

2 ,

theratio of the terms in (22)scale as the square of the time it takesfora barotropic wave totravel a distance Hrelative to the lengthofatime step,andweseethat theconditionforfree sur-faceflows (22)tendstothatoftherigidlid(21)inthelargetime steplimit.So althoughadjustinga systemtoapply afree surface kinematicboundary conditiononthe topsurfaceasopposedtoa rigidliddoesimproveconditioningformodestaspectratios,asthe disparityinscalesincreasesandtheaspectratiobecomessmaller, orequivalentlylargertime steps aretaken, theill-conditioningof arigid lid systemis soonrecovered due tothe quadratic depen-dence.

The multigrid preconditioner of Kramer et al. (2010) for un-structuredmeshesonhighaspect ratiodomainshelpsbetter con-dition the Poisson problem in general, without consideration of WD,using a combination ofalgebraic multigrid anda geometric verticalprolongation operator. This solver methoditself makes it feasibleto run non-hydrostaticunstructured mesh simulations of fluidsingeophysicaldomains.

Whilst the relatively moderate aspect ratio wet areas can be treatedby the multigrid preconditioner approach, specific meth-odstohandletheacuteaspect ratioanddirectrigidlidcondition applied in dry regions are required, if this general fully 3D and non-hydrostaticWD approachistobe appliedtorealgeophysical systems.

4.3.Quantificationoftheill-conditioning

ThediscreteformoftheLaplacianoperatorthatappearsonthe left-hand side of (20), seen in (16), has eigenvalues

λ

i∼ k2i, for

wavenumberski.Theconditioningofthematrixisdetermined by

the ratioofthe maximum

λ

_∞ andminimum

λ

min eigenval-ues.Thisis,equivalently,theratiooftheminimumandmaximum wavenumbers,k_min andkmax,squared

κ



CT_M−1 u C



=





λ

∞

λ

min





∼









k2



∞

k2

min







=





k

_∞

k

_min





2. (23)

ForhighaspectratioproblemsH/L 1,forHandLcharacteristic lengthscalesofthesolutiondomainintheverticalandhorizontal, respectively(seeFig.1),wefind

k

_min∼ 1 H, and

k

∞∼ 1 L, andhence

κ



CT_M−1 u C



∼



H L

2 . (24)

On a spheroid, such as the Earth, the characteristic ‘horizontal’ length scale L is the extent of the encompassing surface geoid, withHthe heightinadirectionparallelto gravitational accelera-tion.Conditioningoflinearsystemthatresultsfromthe discretisa-tionof thePoissonequation isapproximately proportionalto the square ofthe aspect ratioofthe globaldomain.Equivalently, the element edge-lengths, whichare constrainedtoresolve processes importantto thesimulation,can alsobe used tocharacterise the scaling, such that condition number isproportional to (



x/



z)2_, with



xand



zcharacteristicedge-lengthsinlocalhorizontaland verticaldirections,respectively.

Ideally, entries intothe matrixof the linearsystem that arise from dry cells would be removed, in a process similar to lifted Dirichletboundaryconditions(e.g.KarniadakisandSherwin,1999) andthe solverlimited tovariables on thewetsub-system, much like an element removal approach.For an implicitapproach it is not clear howthis wouldbe accomplished without adversely af-fecting the naturalevolution ofthe interface. Instead, under im-plicit integration,treatment ofthe ill-conditioninghighlighted by

(24)needstobeaddressed.

4.4. Verticalvelocityrelaxationindryareas

Toclosethesystem(1)–(2),anequationofstateisrequired. De-tailsoftheformofthisfunctiondonotinfluencethedevelopment thatfollows,andageneraltreatmentoftheevolutionofdensityis considered,suchthat

∂ρ

∂

t +u·

∇

ρ

=0, (25)

withitstemporaldiscretisationfollowingSection3.1as

ρ

n+1₋

_ρ

n



t +u˜

n+1_·

_∇

_ρ

n+1_{=0. (26)}

Developmentoftheapproachproceedswithadiscretisationofthe density transport Eq. (25), in a slightly different linearisation to thatof(26),oftheform

ρ

n+1₋

_ρ

n



t +w n+1

∂ρ

∗n+1

∂

z +s n+1 ρ∗ =0, (27)

ρ

n+1 ∗ −

ρ

n



t +w n+1 ∗

∂ρ

n+1 ∗

∂

z +s n+1 ρ∗ =0, (28)

with vertical velocity w, starred variables representing the best currentguess,andthesourcetermsn_ρ+1 containingdetails of spa-tial gradients of density locally alignedto the geoid. Subtracting

(28) from (27) gives a transport equation, that mirrors (11), de-scribingthevariationoverthePicarditerationprocess

ρ

n+1₋

_ρ

n+1 ∗



t +

(

w n+1_{− w}n+1 ∗

)

∂ρ

n+1 ∗

∂

z =0. (29)

Substitution of this temporally discrete density transport Eq. (29)intothemomentumEq.(7)leadsto

ρ

0 un+1_{− u}n



t =R n+θ₋

_∇

_pn+θ +gng



t

∂ρ

n+1 ∗

∂

z

(

w n+1_{− w}n+1 ∗

)

−

ρ

∗n+1gng. (30)

TheFEweakformof(30)isdevelopedbytestingwithvelocity basis functions

φ

_i andapplying integrationby parts twice atthe freesurfacetoobtain



φ

i

ρ

0 un+1_{− u}n



t d



= 

φ

i



Rn+θ₋

_∇

_pn+θ +gng



t

∂ρ

n+1 ∗

∂

z

(

w n+1_{− w}n+1 ∗

)

−

ρ

∗n+1gng

d



− η

φ

in· ngg



t

(

ρ

n+1₋

_ρ

_a

₎₍

_wn+1_{− w}n+1 ∗

)

d



. (31)

The densityofairjustabovethefree surfaceinterface

ρ

a,canin

mostcasesbe neglectedasasmalleffect,inthesamewayasthe atmosphericpressure.

Assumingthat,forshallowwaters,theverticalvelocitywis lin-earlyrelatedtothedistancefromthebottomoftheoceanorina depth-averagedsense,andignoringthedensityvariations

ρ

inthe surfaceintegralabove,thetermsin(31) abovecontainingexplicit referencetotheverticalvelocitywcanbegroupedintoan absorp-tionterm

σ

zz=g



tmax



0,−

∂ρ

∗n+1

∂

z



+g



t

ρ0

d n· ng, (32)

whered=h+

η

isthewaterdepth,suchthat 

φ

i

ρ

0 un+1_{− u}n



t d



= 

φ

i



Rn+θ₋

_∇

_pn+θ −

σ

(

un+1_{− u}n+1 ∗

)







† −

ρ

n+1 ∗ gng

d



, with

σ

=



_{0 0 0} 0 0 0 0 0

σ

zz



. (33)

The inverse time scale for the vertical velocity relaxation is de-fined by (32).As thePicard iterationsproceed,andun+1

∗ →un+1, themagnitudeofthisstabilisingterm,markedby†in(33),relaxes to zero. Although the absorption coefficient

σ

will be relatively smallinwetregions,andthecontributionfrom†smalloverall,it isimportanttoincludethesetermsthroughoutinorderto main-tainconsistencyandasaresult,accuracy.

The following conditionson verticaldensity gradient,the free surface, vertical viscosity, andvertical absorption provide a well-conditionedpressurePoissonequation

1. Verticaldensitygradient

(

x

)

2

(

z

)

2 ≤ a 2_g

_

_t2_max



∂ρ

n+1 ∗

∂

z ,0



, (34)

2. Freesurfacevariation

(

x

)

2

(

z

)

2 ≤ a2

_(

_t

₎

2 _g d , (35) 3. Verticalviscosity

(

x

)

2

(

z

)

2 ≤ a2

_

_t

ν

zz



z2 , (36) 4. Verticalabsorption

(

x

)

2

(

z

)

2 ≤ a 2

_

_t

_σ

zz, (37)

wherea isa tolerableaspect ratio ofelement length scales (e.g. unityintheisotropiccase),



xand



zcharacteriselocalresolution scales,

ν

zzisakinematicviscosity,and

σ

zzanabsorption.Notethat

theviscosityof (36)mustbe treatedimplicitlyorsemi-implicitly inpressure(e.g.diagonalorblockdiagonal)inordertocontrolthe conditionnumberofthepressureLaplacian.Implementationofthe viscosityinstressformisappropriate heresincetensorforms di-rectlysmoothhorizontalvelocitiesinthevertical.

InthecaseofWD,where(35)doesnothold,we mustensure

(37) is satisfiedby a suitable choice ofthe absorption

σ

zz.From (37), in orderto make the resultingpressure matrix feel like an

O(a)aspectratiodomain,weneed

σ

zz=

(

x

)

2

a2



_t

(

_z

)

2. (38)

Theformof

σ

zzin(38)definestheinversetimescaleforthe

verti-calvelocityrelaxationin(32).Notethatthisformof

σ

zzisspatially

varying,andinparticular,thecharacteristiclocallengthscales



x

and



z are non-homogeneous across the geoid surface. In WD simulations,thesefieldscontainlargedeviations,indicativeofthe regionsaffectingconditioningofthepressurePoissonequation.

4.5.Discretisationforhighaspectratiodomains

ThemomentumEq.(33)discretisedinspace-timeatanygiven Picarditerationstepis

ρ

0 Mu



t



un+1 ∗ − un



+

θ

_A˜n+1_un+1 ∗ +

(

1−

θ

)

Anun =

θ

C pn_∗+1+

(

1−

θ

)

C pn− Q



un+1− un+1 ∗



+Su, (39)

withQi j:=_

φ

i

σφ

jd



.Thisbalancecomparedtoitsendoftime

stepstateismultipliedby

θ



tCT_M−1 u ,togive



ρ

0

θ

CT−

θ



tCTMu−1Q



un+1_{− u}n+1 ∗



=

θ

2

_

_tCT_M−1 u C

(

pn+1− pn∗+1

)

. (40)

Thisisequivalentto(11)previously,notingthetermmarked†in

(33)iszeroattheendofatimestep.

Thenew formofthe combinedkinematic boundarycondition

(19)leadstoatimediscretisedformmodifiedfrom(13)toinclude theno normalflow component appliedin dryregions, described by nn+1_{· ngmax}



_pn+1_,

_ρ

_0g

₍

_h+d0

₎



− nn_{· n} gmax

(

pn,

ρ

0g

(

h+d0

)

)

=

ρ

0g



t



θ

nn+1_{· u}n+1₊

₍

₁₋

_θ

₎

_nn_{· u}n



_.

Moreover,thesurfaceintegralMsismodifiedsuchthatthediscrete

modifiedkinematicconditon(14)becomes

Mw pn+1_{− p}n

ρ

0 g



t + Md h+d0



t =G T θun+1+GT1−θun,

whereMw,i j=wng

ψ

i

ψ

jd



andMd,i j=dng

ψ

i

ψ

jd



.This kine-maticconditionchange modifiesthepressure correction, andthe discretecontinuity(15)becomes

θ

CT_un+1₊

₍

₁₋

_θ

₎

_CT_un_+M wp n+1_{− p}n

ρ0

g



t +Md h+d0



t =0. (41) Substituting (41) into momentum (40), yields the discrete com-binedp

(

p,

η

)

Poissoncorrector,with(16)evolvingto



θ

2_CT_M−1 u C+ Mw g

(

t

)

2



p =−

θ

CTun∗+1+

(

1−

θ

)

CTun



t − Q



un+1 ∗ − un+1





t − Mw g

(

t

)

2

(

p n+1 ∗ − pn

)

−

ρ0

Md h+d0

(

t

)

2. (42)

The predictor – corrector method of Section 3.5 solves the non-linearsystemwiththeupdatedp

(

p,

η

)

Poissoncorrector(42) com-binedwithdiscretelinearisedmomentum(39),andavelocity cor-rectiondeterminedfrom(40).

4.6.Self-consistencyandphysicalbasisofthesolution

Mass, momentumandtracerquantities areself-consistent and conserved,propertiesinheritedfromtheir underlyingGalerkin FE weak formulations (Piggott et al., 2008) and use of a thin-film (Funke et al., 2011; Medeiros and Hagen, 2013). The constrained discreteSoblevsolutionspaceoftheweakformmodifiedwiththe additionaltermsmarked†in(33)convergesonthesolutionspace oftheoriginalformwithoutthese,asthePicardprocessproceeds. Inasimilar mannertoPetrov–Galerkinandvariationalmultiscale (Hughes etal., 1998) residual-basedstabilisationmethods (Candy, 2008),thisensures consistency,that thesolution foundisa valid solutionoftheoriginalweakGalerkinformulation,atruediscrete solution to the governing continuum equations and is therefore physically-based.

Moreover, justlike streamline-upwind Petrov–Galerkin(SUPG) stabilisation, the additional terms themselves are defined from physicalpropertiesoftheflow. Forexample,(32)includes contri-butionsfromg,the localvertical densitygradientandwater col-umndepth. This is supported along with local discretisation pa-rameterssuch astimestepandelement sizeusedtoquantify un-resolvedscales,inasimilarwaytomultiscaleturbulenceclosures (Candy,2008).

4.7.Determinationofcharacteristiclengthscales

Accuratecalculationofthecharacteristiclengthscalesiscritical tothe success of the approach, particularlydue to the quadratic dependencein(38).

The calculation of the characteristic horizontal length scale couldbesimplytheminimumormaximumedgelengthofthe el-ementprojectedtoa2Dhorizontalgeoid.Amoreaccurate approx-imationcanbedetermined fromthesmallestandlargest circum-scribingcircularboundsofthisprojection.The lengthscale



x is thenafunctionoftheseminimumandmaximumextents.Thisisa naturalapproachformodelsemployinganisotropicmeshelements. The verticallengthscale islessambiguoustodetermine,since uniqueintersectionswith



_ηand



bexist

∀

x∈



,duetothe

con-struction of geophysical domains Candy (in prep.); Candy et al. (2014), and similarly for internal layers. Evaluating length scale functions at Gaussian quadrature points rather than by element furtherincreasesaccuracy,sinceFEassemblyintegrationsare per-formedthisway,withoptionsto develop meanorarea-weighted means. This is trivially extended to superparametric elements whicharetypically usedinthe toplayerforaccurate representa-tionofgeoidcurvature.

Arguably the best characterisation of tetrahedral element size is determined from the Jacobian transformation matrix which projectsaFEfromglobaltolocalparameterised space.The deter-minantof the transformation Jacobian intersected with the local (to quadraturepoint)surfacegeoid plane andgravitational accel-erationvectorwillgivecharacteristiclengthscalesfortheelement intherequiredhorizontalandverticaldirections,respectively.This approachalsonaturally handleselement anisotrophyandmeshes whichare fullyunstructuredin 3D.The meritsofthis choiceare examinedinSection7.3.

4.8.Correctiontovelocityrelaxationinshallowregions

UndernoforcingthemomentumEq.(33)tendstorelaxthe im-plicitvelocity un+1

∗ tothe statein theprevious time stepun,but

thiscan betoo stronginvery shallowareas.Thisiscorrected by reducingthemagnitudeoftheexplicitpartofthevelocitythatwe relaxto,byadding−

γ

un _{totherightofthemomentumequation,}

with

γ

=max



2



1−₂d_d 0



, 0



. (43)

Forawatercolumndepthd,thisrelaxationscalesawaythe veloc-ityinthevicinityofdryregionswhered_<2d0,andrelaxestozero indryregions,whered=d0.

5. Meshmovementwithwettinganddrying

5.1. Discretefunctionspaceupdates

Thefreesurfaceevolutionresultsinmanyquantitiesvaryingin time, such asthefree surfacenormalvector n in (13).Moreover, this includes the mesh, and hence spatial discretisation, which leadstoachangeofthediscretefunctionspacesSh_,andtheir

span-ningbasissetsresultinginnewformsofmassandothermatrices indiscreteformssuch as(14).Forconservationandaccuracyitis necessarytoupdate thediscretenon-linearsystemduringthe Pi-carditerationtoreflectthesechanges.Therearevarioustechniques tohandle thisconservatively,through thedefinitionofa grid ve-locity,for example.Inthisformulation,the domain discretisation isupdated attheendofaPicard iterationtoreflectthenewfree surfaceheightpredicted,withthenormaln,massmatrixandother matricesrepresentingadvectionandsurfaceintegralsrecalculated underthenewdomaindiscretisation.It isthereforethecasethat thediscretematricesMu,Ms,C,AG,andQ;free surfacenormaln,

basis functions

φ

and

ψ

,anddomain



andfree surface



_η,are always thebest knownapproximation, i.e.thestarred n+1 case. AsubtleexceptionisattheendofthefinalPicarditeration,where theupdate isnot made, toensurethe domainandderivative pa-rametersarethosetheprognosticvariableswerecalculatedon.

The generalised approach that includes the non-linear advec-tionterminthegoverningEq.(1)precludesthediscretisedspatial operatorfrombeingself-adjoint.Evaluationofthisnon-linearterm requires sub-cycling,and forunder-resolved highFroude number orrapidly-varyingflows thiscould requirealargenumberof iter-ationsto converge, unlessthe continuum systemis linearised,or localresolutionincreased.

5.2. Surfacerepresentationandinterfacetracking

At the end of each Picard iteration, as outlined in Section 3, thefree surfacepositionis updated using(18)to reflectthenew pressure atthe interface pn+1

_|

_{ηs.Due to the minimumthreshold} d0,theperturbationoftheinterfaceinthedirectionofthe gravita-tionalaccelerationislimited. Ifthepressure pn+1 _{atthe interface} impliesitshould move belowthislevel,itis fixedatthe thresh-old levelabovethebottombathymetry(i.e.

η

=h+d0).The pres-sure remains unaffected, and is allowed to deviate from the in-terfaceposition

η

. Conversely, assoon as pn+1 _{produces a water} columndepthgreaterthand0 thefreesurfaceinterfacemoves up-wards.Correspondingly,thedomaindiscretisationisupdatedwith the mesh stretched in the direction _−ng, parallel to the

gravita-tionalvector,tomeetthenewfreesurfacebound.

Spatialrepresentationof

η

isinheritedfromthefunctionspace usedtoapproximatethecombined p

(

p_,

η

)

.Irrespectiveofthe or-der of variation,such asquadratic forthe PDG

1 − P2 element pair, theinterfaceisapproximatedbyapiecewiselinearfunctionasfar asthedomainrepresentationisconcerned.Thissatisfiesthe min-max property, such that theextent of thesurface isbounded by the nodal positions that define its representation. This, together

withtheminimumthresholdlevelpreventselementsfrom becom-inginvertedorexcessivelysmall.

5.3. Remeshing

It is not arequirement that thedomain is remeshedanewto theseadjusted bounds. Since only one ofthe domain boundaries isperturbedthroughtheabove processandinadirectionaligned tothegravitationalvectorfield,locallyorthogonaltootherbounds ofthedomain,itispossibletoapplyarelativelysimpler-adaptive transform. Thedomainmeshisstretched linearlyinthisdirection to fitthenew boundary.It isalsopossible tolimit the perturba-tion to the nodeson the free surface, or to apply more compli-catedr-orh-adaptivestrategiestoachieveahybridisedcoordinate system(see Kleptsovaetal., 2010;Bleck, 2002;Burchard and Pe-tersen,1997)formoreaccuratesolutionsorbetter-represented fea-tures. Theimplementation ofthe approachdescribed hereinthe modelcode(Fluidity,Piggottetal.,2008) functionswithand sup-portsthesemethods.

6. Additionalstabilisingapproachesfordryareasinhigh aspectratiodomains

Twosupplementaryapproachestocontrolconditioningare pre-sented, acting directly to prevent strong erroneous flows devel-oping in the thinfilm andmodifying behaviour in neighbouring wet regions. This is exacerbatedby the fact the physical system issolvedinaweaksense,whichwhilstbetterforconditioningcan permitlargefluxesacrosstheinterface.Unliketheabove,these ap-proaches transformthe solution space and have the potential to affectthe solutionin unphysicalways.Theyare presentedas ad-ditionaltechniqueswhichcanbeemployedtoenableasolutionto bereached,butrequirecarefulapplication.

6.1. Manning–Stricklerdraganddryregionstability

InthecaseofinundationflowswhereWDisapplied,a param-eterisation of drag that is commonly employed is the Manning– Stricklerformulation,definingthebottomstress

n·

μ∇

u=n2 _g

|

u

|

u

d1/3, on



b, (44)

where n is the Manning coefficient, d is the water depth and n here is the unit surface normal on the bottom surface



_b. This formulationitselfhasastabilisingeffect,andmoresointhevery shallowdryregions,withadragappliedalongthebottom bound-ary proportional to d−1/3_{. Inpractice, theManning–Strickler} bot-tom stress is sufficient to prevent significant erroneous flow de-velopingindryareas.Inthecasesofacutehighaspectratio,long time steps or particularly steep bathymetric gradients, the Man-ningcoefficientcanbeincreasedindryregionsandtheirproximity toincreasethestabilisingeffect,with

ˆ n=n+max



0,ndry 2d0− d d0



,

wherenˆreplacesnin(44),andforndryanewManningcoefficient (withusual standardunitsofsm−1/3_{) largeinsize,relativetothe} standardcoefficientn.

6.2. Horizontalbulkeddyviscosityindryregions

Asecondsolutiontoincreasestability,istodampflowdirectly indryregionswithabulkvolumeviscosityorasource-absorption sponge,bothallowingtheapproachtoremainimplicit.

This stabilisation is applied throughoutthe domain, or selec-tivelyindryregionsandtheirimmediateproximity,withthelarge horizontalviscosity

ν

L=max



0,

νdry

2d0− d d0



, (45)

introducedto controlspurious horizontalfluxes, with

ν

dry a con-stanteddyviscositycoefficientandd_{≥ d}0

∀

x∈



.Thishorizontal viscosity is continuous in space without discontinuous jumps in intensityacross the WD interface, actingin theproximity of dry regionswhered_<2d0.

7. Validationandapplication:Numericaltests

Performance of the implicit WD formulation described in

Sections2–5,andadditionalstrategiesof Section6 areexamined infourtestscenariosinacutelyhighaspectratiodomains,toa de-greefoundingeophysicalsystems.

7.1. Implementationandverification

The approach has been implemented and validated in the FE fluiddynamicscodeFluidity(Piggottetal., 2008).Thissimulation frameworkcontainsmanytoolsforgeophysical modelling,is par-allelisedwithsophisticated load balancingand supports adaptive meshmethodsallowingcomputationalefforttobe focusedon re-gionsofdynamicinterest.Itfunctionsforaspatiallyvariable grav-itationalaccelerationvector,andhencecanbeusedforlarge-scale simulationson theEarth’sspheroid.Theimplementation includes asuiteoftestcasestoroutinelyverifythenewalgorithmina for-mal sense, in an automated continuous verification build engine (Farrelletal.,2011)toensurerobustnessofthecodeandresiliency inlightoffurtherdevelopment.Theunstructuredmeshesusedin thefollowingcaseswere builtbymeans oftheopen source soft-wareGmsh1_.

The balance and LBB stability properties of the PDG 1 − P2 velocity-pressurepairing(see Cotteretal.,2009) aidconditioning andareusedinallapplicationsconsideredhere.Allfourcaseshave been run on the purely continuous pairing P1− P1 also,but due tothepressurefilteringrequired,didnot performaswell,andin allbutmodestaspectratiocasesweretooill-conditionedtoreach convergence.ThebehaviourofPDG

1 − P2andP1− P1 andtheir rela-tive performancein regularaspect ratioproblemsis presentedin

Cotteretal.(2009).

Duetotheaspectratiosconsidered,allcasesusethemultigrid preconditionerdescribed in Krameret al. (2010) foriterative so-lution ofthe conditioned symmetric pressure Poisson linear sys-tem in combination with Conjugate Gradient (CG, Hestenes and Stiefel, 1952). The momentum system is solved in a more stan-dardapproach withSymmetric SuccessiveOver-Relaxation (SSOR,

Young, 1971) preconditioning and the iterative Restarted Gener-alised Minimal Residual (GMRES, Saad and Schultz, 1986) algo-rithm, where the calculation is restarted after k=30 iterations. The iterativeSSOR-GMRES process isperformed usingalgorithms built into the established and well-verified PETSc library (Balay etal., 1997). In contrastto the study (Funke etal., 2011), itwas foundthat twoPicard iterationsprovide sufficientconvergenceof thecoupled systemin thecases studied.In all cases,both linear systemsaresolvedtoaconvergencecriteriaspecifiedbyarelative errortoleranceof 10−7,which isconsidered sufficientlyaccurate. ThequadraturebasedsubgridresolutiondescribedinFunkeetal. (2011)isalsoused,withaquadraturedegreeofeight.

1.0 0.5 0.0 -1.0 -1.5 -0.5 -2.0 1.4 0.4 0.0 0.2 0.6 0.8 1.0 1.2 1.0 0.5 0.0 -1.0 -1.5 -0.5 -2.0 1.4 0.4 0.0 0.2 0.6 0.8 1.0 1.2

Fig. 2. The first Balzano channel flow benchmark with, in this presented case, a horizontal extent of 1.38 × 10 4 m, corresponding to a minimum element aspect ratio of

∼ 10 −6 . The discretised horizontal surface (a) contains 108 triangular elements with a characteristic length scale of ∼ 500m. Vertical sections show the free surface position

at 10min intervals for the initial drying phase (b) and during wetting (c).

7.2.FirstBalzanoslopedchannelbenchmark

Thefirsttwosetsofnumericaltestsarefromthesuiteof prob-lems inBalzano (1998), selected since they exhibit the problem-atic ill-conditioningin as simple a setup as possible. No analyt-ical solution is available, so the problem configuration is chosen consistentlywithBalzano(1998)tobe abletodraw comparisons. Thebasebenchmarkcaseisdevelopedfromtheoriginally2D do-mainconsistingofa13.8kmlongslopewithadepthof5matone endwhichtendstozeroattheother.Recentlydevelopedschemes, suchastheflux-limitingWDmethodforFESWEmodelspresented inGourgue etal. (2009) and the non-hydrostatic algorithm pro-posed in Funke et al. (2011), have been benchmarked on these cases.These modelin 3D, butforce dynamics to occur predomi-nantlyinthedirectionswherethe extremesinextentoccur, with 10elementsintroducedinthethirddirectionintheformerand1– 2inthelatter,whichisfollowedheretoawidthof1km.Withthe assumptionsolutionsarelaminar,thisextrusioninto3Dspacewill notchangethephysicalbehaviour.Theslopedbottombathymetry isdefined h

(

x

)

₌x_/2760_, forthe x-coordinatedirectionindicated alongsidethe surfacegeoid computational meshin Fig. 2(a).The basecasesingle-layermeshcontainsvertically-aligned nodesand ahorizontalelementsizeof500m.

Following the benchmark description in Balzano (1998) (also in Gourgue et al. (2009)), no normal flow boundary conditions are applied at the bottom and shallow end of the domain, and additionally applied to the sides. A Manning–Strickler drag with

n=0.02sm−1/3 _{is applied at the bottom boundary. The} gravita-tionalacceleration isset to 9.81ms−2 and thefluid is initially at rest.Timediscretisation is performed withCrank-Nicholson inte-gration(i.e.

θ

=1

2) anda time stepof 600s. Inthiscasethe WD thresholdis set atd0=0.5mm. The free surface is forcedat the deepopen boundarywithasinusoidalvariationofamplitude2m, suchthatwatercolumnthicknessoscillatesbetween3–7m,witha periodof12h.

In the seriesoftestsconsidered here, the horizontalextent is varied from1.38× 102_{m to 1.38× 10}6_{m, centredabout the} de-finedbenchmarklength of1.38× 104_{m. Thisprovidesa rangeof} elementaspectratiosfrom10−4 to10−8_,adomainaspectratioup to3.62× 10−6 _{andspatialscales spanningover10ordersof} mag-nitude in a single domain. Element lengths are scaled with the domain length, such that element aspect ratio relative to global aspect ratio is maintained, with the extrusion in the third di-rection also scaled to preserve element shape. The time step is alsoscaled to ensure the wave Courant numberis constant. The WDthresholdd0,andvertical extentare kept constant acrossall cases.

Fig. 3. Impact of the optimal aspect ratio parameter on solver conditioning in the first Balzano benchmark over a WD phase. 101 individual simulations shown.

Thefreesurfaceevolutionoftheintermediatecasewitha hor-izontal extent of 1.38 × 104_{m is shown in Fig. 2 at 10min} in-tervals, matching (Balzano, 1998) and(Gourgue et al., 2009), for theinitialdryingandthenwetting phase,respectively.Theresults are physically reasonable, and comparable to other formulations (Funkeetal.,2011andGourgueetal.,2009forexample).In partic-ular,thefreesurfaceinterfacedoesnotsufferfromeither underes-timationwithnegativewatercolumnthickness,norproduce oscil-lationsduringthewetting processobservedinBalzano(1998)for someofthe10methodsexamined.Thisbehaviourischaracteristic ofthesolutionsacrosstherangeofaspectratios.

Through a modification of the optimum aspect ratio parame-ter a in(38),there is acorresponding change inthe aspectratio feltinthediscretepressurematrixofelementsindryregions.The parameter a isvaried overthe range a∈[10−4,104_{] ina suiteof} 1001 simulations ofthe baseBalzano case. Solver iteration num-ber is used asan indicator of conditioning,andplotted in Fig.3

forboththepressureandvelocitycalculationsasmeanand maxi-mumvaluesoverthecourseofaWDphase.Theparameterrange hasbeenspacedequallyinlog-space inordertogive agood rep-resentationofthebehaviouroverthelargerangeofdomainaspect ratios.Thisisachievedwithadiscreteparameterspacedefinedfor