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

A new approach to modelling free surface ﬂows is developed that enables, for the ﬁrst time, 3D consis- tent non-hydrostatic baroclinic physics that wets and drys in the large aspect ratio spatial domains that characterise geophysical systems. This is key in the integration of physical models to permit seamless simulation in a single consistent arbitrarily unstructured multiscale and multi-physics dynamical model. A high order continuum representation is achieved through a general Galerkin ﬁnite element formulation that guarantees local and global mass conservation, and consistent tracer advection. A ﬂexible spatial discretisation permits conforming domain bounds and a variable spatial resolution, whilst atypical use of fully implicit time integration ensures computational eﬃciency. Notably this brings the natural inclusion of non-hydrostatic baroclinic physics and a consideration of vertical inertia to ﬂood modelling in the full 3D domain. This has application in improving modelling of inundation processes in geophysical domains, where dynamics proceeds over a large range of horizontal extents relative to vertical resolution, such as in the evolution of a tsunami, or in urban environments containing complex geometric structures at a range of scales.


Introduction
Flooding has huge impacts on the economy of a region and the livelihood of its people. Significant progress has been made to model and predict the impact of these inundation events, capturing the character of their source and resultant behaviour. Many challenges still exist and in particular in concurrently simulating the physical processes involved from the large planet-scale forcings down to the small human scales of an urban environment. This is highlighted in the review ( Medeiros and Hagen, 2013 ) as one of the key limitations of existing wetting and drying (WD) models. In an urban flooding scenario for example, modelled water column depth could be down to 1cm over a horizontal range of tens or hundreds of kilometres, leading to a very high aspect ratio of ∼ 10 −7 .
Inundation flow models typically use simplified formulations of the Navier-Stokes equations, commonly the Saint-Venant shallow water equations (SWEs). These simplifications make assumptions, such as a hydrostatic pressure and well-mixed water column, which are not necessarily valid in the whole range of scales relevant to the inundation. Non-hydrostatic processes become important, for example, in the dispersive effects of short waves where the ratio of vertical and horizontal scales of motion are not sufficiently small. The study of Oishi et al. (2013) considering the 2011 T ōhoku tsunami in Japan, found it was critical to include nonhydrostatic effects to correctly model processes on the small scale, a point further highlighted by Cui et al. (2014) .
Typically, multi-physics over a broad range of scales is approached using multiple model runs at a hierarchy of scales such that domains are nested, with varying complexity and included physics. As an alternative, effort s to integrate the physics and scales of separate models into single Earth system models is growing, where it is important individual components function in a general context, and are not too restrictive in discretisation choice. Although this can be achieved weakly with offline communication between models, the 'holy grail' is a flexible single model capable of simulating a range of physics and scales, with inherit consistency and conservation.
This work pushes the boundaries in two key regards: Firstly, adding a novel approach to WD in the 'thin-film' family solving a full 3D pressure rather than the usual SWE approximation in very challenging acutely large aspect ratio domains typical of geophys- ical systems -a first for WD. Secondly, this brings the modelling of WD processes together with non-hydrostatic baroclinic flow dynamics in a single simultaneous and seamless system model. This is critical for tightly coupled processes, for example in tracking grounding line movement under an ice shelf ocean cavity, that is strongly influenced by non-hydrostatic and baroclinic ocean flows.
Accurately tracking an inundation interface is technically challenging. Of the Eulerian type WD approaches (reviewed in Medeiros and Hagen, 2013 ), where the underlying spatial discretisation is predominantly independent of space and time such that matrix operators can be cached and there is no need for complex contour tracking, there are four types: element removal, limiting the computational domain to the wet region (see Casulli and Stelling, 1998 , UnTRIM Casulli and Walters, 20 0 0;D'Alpaos and Defina, 20 07;Defina, 20 0 0 and the WASH123D code Lin et al., 2004 ); thin film approaches (see Bates and Anderson, 1993, the FVCOM model Chen et al., 2003, POM model Oey, 2005 and also Begnudelli and Sanders, 2006 ); depth extrapolation from wet to dry cells ( Bradford and Sanders, 2002;Lynett et al., 2002 ); and negative depth ( Heniche et al.,20 0 0;Jiang and Wai,20 05 ) applied in ROMS ( Warner et al., 2013 ), including the use of porous media below the sea bed ( vant Hof and Vollebregt, 2005;Ip et al., 1998 ) and bathymetry movement ( Kärnä et al., 2011 ).
Underlying model discretisations largely steer this choice, with the first by far the most common for explicit time stepping models, applied in QUODDY, ADCIRC, MIKE21, Delft3D and in one of the first Eulerian methods ( Leendertse, 1970 ), subsequently reviewed in Balzano (1998) . Whilst robust, stability constraints restrict movement of the interface to one cell per time step ( t ), since the Courant number must be maintained less than one in drying regions ( Stelling and Duinmeijer, 2003 ) to ensure a nonnegative bound on water depth, a strict limitation on t ( Walters, 2005 ). Depth extrapolation also suffers this restriction with elements switching states ( Medeiros and Hagen, 2013 ), whereas thin film and negative depth options can be time-integrated implicitly.
Current approaches to unstructured mesh geophysical fluid modelling is considered in detail in Danilov (2013) , with its potential importance best highlighted in Danilov et al. (2013) . Indeed, this review states that whilst unstructured mesh models may not replace structured modelling approaches completely, there are cases where this type of approach could be optimal. In particular, allowing a flexible approach to the vertical discretisation could improve accuracy and model efficiency in domains where there are sharp changes in bathymetry relative to horizontal spatial resolution, strong non-hydrostatic gradients in pressure, strong vertical inertial flows, or when it would be more optimal to reduce or increase the number of layers in shallow and deep regions, respectively. Moreover, these could be critical in the fringes of the ocean boundary, along geometrically complex coastlines and in interactions with other types of physical systems, such as an urban environment, or the complex shallowing in ice shelf ocean cavities. Within this discretisation type, WD models can more accurately model a wider range of scales in larger domains.
One of the early finite volume (FV) approaches UnTRIM ( Casulli and Walters, 20 0 0 ) permits unstructured meshes with the constraint that, like structured models, the domain elements are orthogonal where circumcentres are inside their respective elements. Its non-hydrostatic advance ( Casulli and Zanolli, 2002 ) is applied in the SUNTANS model, with the same orthogonality restriction. It contains a WD algorithm ( Wang et al., 2009 ) stabilised with a technique from Ip et al. (1998) that applies an increased drag to satisfy an additional constraint on volume flues in dry regions. Similar approaches are also made in Cui et al. (2010 , FVCOM, MIKE21) with non-hydrostatic corrections added (e.g. Cui et al., 2012 ). FV is low order only and models generally explicit.
Unstructured finite element (FE) methods offer high order continuum approximations which are more accurate and naturally include less diffusive and dispersive advection schemes. WD has been built into 2D barotropic flow models such as QUODDY with dry element removal in Greenberg et al. (2005) ; ADCIRC, a SWE method for explicit hydrostatic modelling of storm surges with dry removal ( Dietrich et al., 2006 ); TELEMAC, initially using element removal ( Bates and Hervouet, 1999 ) and now negative depth; and SLIM ( Kärnä et al., 2011 ) with a repositioned sea bed SWE method and adoption of implicit t advance.
WD is combined with solvers capable of modelling baroclinic processes in Funke et al. (2011) ;Warner et al. (2013) , with the former using thin film high order FE and the latter explicit finite difference with negative depth WD and mode splitting. The former performs well in relatively modest aspect ratio domains, but performance is strictly limited by use of direct solvers (LU decomposition), restrictions on dry element aspect ratios and erroneous unphysical flows that develop in dry regions.
Here a general approach for WD with FEs is considered in full 3D, building on established methods for modelling fluid flow on fully 3D unstructured meshes ( Piggott et al., 2008 ) which vary in resolution and support a multiscale of physical processes, including non-hydrostatic and baroclinic dynamics in the large aspect ratio domains found in geophysical domains. Under the constraints of a global number of degrees of freedom, this allows the focus of computational resources on small scale regions and areas of interest, whilst capturing the large scale flows elsewhere in the domain. An additional advantage is no constraint on the internal mesh structure, nor that it is fixed in time. It is not constrained to layers, and can be completely, or partially in select regions, fully anisotropically unstructured.
To allow efficient time integration over a range of element sizes, an implicit treatment necessitates a continuum approach to interface tracking. A thin film is applied, which as ( Medeiros and Hagen, 2013 ) notes, generally satisfies mass and momentum conservation without significant special treatment, and produces a realistic and smooth wetting front. WD is included in a natural manner, through additional terms in the momentum equation and modified boundary conditions. Indeed, the numerical treatment is careful to ensure the solution remains in the Sobolev solution space of the original physically-based weakly formulated Galerkin problem. Prognostic variables, including tracers, are self-consistent through the FE formulation and notably, through use of a combined pressure variable, consistency with the free surface is naturally inherent.
In the following Sections 2 -4 , the new consistent approach for WD in large aspect ratio geophysical domains is developed, with details of mesh movement in Section 5 and additional strategies noted in Section 6 . This is validated in Section 7 and evaluated in Section 8 .

3D Boussinesq with piezometric pressure
The non-hydrostatic Boussinesq equations for a rotating stratified fluid, are solved in a time-dependent domain ⊂ R 3 (see Fig. 1 ), bounded by the surface . This is split into the free surface boundary η , and the remaining bound b = \ η . These are defined for the prognostic variables of velocity u :  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.
where μ is the tensorial dynamic viscosity, −n g and g the gravitational acceleration direction and magnitude respectively, and ρ : × [0 , T ) → R the density. The latter is split into a background ρ 0 , and perturbation density ρ , such that ρ = ρ 0 + ρ . Since the hydrostatic component of pressure of the equilibrium state does not have an important contribution dynamically, it is subtracted from the momentum equation and the full pressure p , is replaced by a piezometric pressure, commonly applied in coastal engineering applications (e.g. Labeur and Pietrzak, 2005 ), defined as p := p + ρ 0 g n g · r + p a , for a position vector r , relative to a position where hydrostatic pressure is zero. Atmospheric pressure at the interface is denoted p a .
Redefining the prognostic pressure with this particular choice of piezometric pressure forms a combined free surface -pressure prognostic p( p , η) eliminating the need to solve a separate, commonly used, wave equation for the free surface, denoted by the injective function η: × [0, T ) → η . The prognostic pressure p now contains non-hydrostatic components and the hydrostatic pressure due to perturbations in the free surface elevation. This remaining hydrostatic pressure ρ 0 g η, is the boundary condition for p at η , and through (3) , we find p η = ρ 0 gη. (4)

Boundary conditions
With the inclusion of the free surface height in the prognostic pressure, the kinematic free surface boundary condition of Appendix A is expressed n · n g ∂ p ∂t η = ρ 0 g n · u , on η .
(5) This is the boundary condition for η and now a required constraint for the combined p( p , η) prognostic variable. This is joined by the u and p boundary constraints More general conditions, for open ocean boundaries or flux inputs, can be applied without fundamental changes to the approach.

Coordinate system and frame of reference
Note additionally that the direction of gravity, describing the normal n g , is not restricted to a Cartesian z -component, such that the development is relatively independent of the coordinate reference frame. It is free to vary arbitrarily within R 3 , aligned with the local direction of gravitational acceleration, and it is possible for example, to apply this method to the spheroid shell of the Earth in a Cartesian coordinate reference frame.

Spatial and temporal discretisation
The non-linear system of Eqs. (1) and (2) , combined with boundary conditions (5) and (6) , are solved for p( p , η) , and velocity u , using a Chorin projection method ( Chorin, 1967 ) to enforce incompressibility. This is a modified predictor -corrector scheme based on Gresho et al. (1984) in which a predictor u n +1 * is obtained from momentum conservation that is not divergence free, such that a correction u n +1 = u n +1 * − ∇ φ is then calculated subject to the divergence-free constraint ∇ · u n +1 = 0 . For each time step, this proceeds for a number of Picard iterations until sufficiently converged.

Temporal discretisation
Discretisation in time is achieved by the θ -method ( Iserles, 2012 ) in all cases, for a time step t , such that the explicit forward Euler, Crank-Nicolson and backward Euler time-stepping schemes can be obtained with choices of θ = 0 , 1 2 and 1, respectively. The modified Navier-Stokes with implicit free surface system (1) and (2) at a time n is therefore where R n + θ = θ R n +1 + (1 − θ ) R n contains the advective mass flux term, together with viscosity and other source terms. A choice θ ∈ [ 1 2 , 1] leads to an implicit time stepping scheme that allows simulations to use large time steps, which are not restricted by the Courant-Friedrichs-Lewy (CFL) condition ( Courant et al., 1928 ) with respect to the velocity and wave speed. In practice for the simulations presented here, for the required level of accuracy and stability, Courant numbers up to 10 are applied.

Combined free surface -pressure Chorin corrector
Under a Galerkin FE spatial discretisation the temporally discretised momentum (7) and continuity (8) equations are tested with the velocity φ i and pressure ψ i basis functions, respectively. The trial functions u and p are defined in terms of their respective basis functions also, and Appendix B describes their form and the nomenclature used in more detail. This leads to the space-time discrete momentum equation for a Picard iteration step and end of time step, respectively. The starred variable u n +1 * is the current best approximation to u n +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 ˜ u n +1 , and used in the calculation of updated non-linear operators, such as mass flux ˜ A n +1 . The velocity space mass matrix M u additionally contains the diagonal or block-diagonal (depending on the chosen discretisation) component of viscosity from R , which is to be treated implicitly in pressure. The discrete cross-space gradient operat or C i j := − φ i ∇ ψ j d , contains an inner product over velocity and pressure spaces, leaving sources in S u .
Subtracting (9) from (10) and multiplying by θ t C T M −1 u yields a discrete Poisson equation for the correction
For G θ = G (1 −θ ) = 0 , the system of equations enforces incompressibility with weakly applied no normal flow boundary conditions.

Discrete modified kinematic boundary condition
Discretisation of the free surface boundary condition (5) is now required to provide the boundary integral terms in (12) , with a θ time discretisation described by n n +1 · n g p n +1 − n n · n g p n = ρ 0 g t θ n n +1 · u n +1 + (1 − θ ) n n · u n .
Discretisation of (13) in space using the test and trial functions φ i and ψ i , introduced in Section 3.2 gives M s p n +1 − p n with the surface integr al M s,i j := η n g ψ i ψ j d . Applying this discrete combined p( p , η) kinematic condition (14) to the discrete continuity (12) , we find Substituting (15) into the momentum Eq. (11) , with the correction defined p := p n +1 − p n +1 * , yields 3.5. Combined free surface -pressure system solution During a single Picard iteration, the first velocity predictor step solves the discrete linearised momentum Eq. (9) , to establish an updated intermediate velocity u n +1 * , from the current best approximation to the velocity and pressure, and their value at the previous time step.
The predictor u n +1 * obtained is not divergence free in general, and in order to enforce the incompressibility condition, a pressure correction p is calculated to project this velocity into the divergence free subspace by solving (16) above. The velocity correction is made to update the intermediate velocity, consistent with the new intermediate pressure, and projected to the divergence free subspace using the difference of (9) and (10) , where Finally, the interface tracking step adjusts the free surface position following (4) in light of the new pressure field, in a direction −n g , parallel to the gravitational vector.

Wetting and drying of the simulation domain
The free surface boundary is split into distinct wet and dry regions (illustrated in Fig. 1 ), defined by the combined p( p , η) at the The conditions in dry regions differ from those in wet in two defining ways. Firstly, the water column depth is maintained at a threshold minimum d 0 above the bottom bathymetry defined by h : → R , and secondly, the surface boundary condition on the combined p( p , η) prognostic variable is modified to enforce this constraint in the solver. With the depth constraint the free surface evolution described by (4) provides the first constraint The second is found by modifying the combined kinematic condition (5) , which in light of the depth restriction gives In wet regions w , the constraints (18) and (19) reduce back to the free surface conditions (4) and (5) , respectively. In dry regions w , (19) is a no normal flow condition, and effectively im poses a rigid lid approximation.

Conditioning of the pressure calculation
The spatial domains of geophysical processes are typically large aspect ratio, due to the gravitational influence and disparity in dynamics parallel and perpendicular to geoid surfaces. The solution of a non-linear fluid flow system in these types of domains including non-hydrostatic dynamics, with implicit time evolution and a predictor -corrector approach such as Section 3.5 is shown in Kramer et al. (2010) to lead to an ill-conditioned pressure system. In the limit of large domain aspect ratio and long time steps, the system behaves approximately as though it has a rigid lid, where the free surface is fixed with η = 0 and u · n = 0 .
The dry regions introduced by the WD process significant exacerbate ill-conditioning, since a rigid lid condition is applied directly and the region contains elements with acutely large aspect ratios due to their defining shallow water column depth.
The correction u n +1 = u * − ∇ φ (17) is calculated subject to the divergence-free constraint ∇ · u n +1 = 0 . This leads to the following pressure Poisson equation for φ which corresponds to the discrete Poisson operator C T M −1 u C in the formulations (16) and (42) above. For no normal flow boundary conditions where u · n = 0 , at interfaces with bedrock or in the case of the rigid lid approximation for the ocean-air interface, the coupling between velocity and pressure results in the corresponding boundary condition on (20) as the Neumann expression ∂φ ∂n ensuring the velocity constraint is consistently preserved. Applying a kinematic condition instead leads to the homogeneous Dirichlet condition φ = 0 on η . The redefinition of pressure in (3) to form the piezometric pressure here allows standard pressure splitting approaches to treat baroclinic and barotropic dynamics (e.g. Shchepetkin and McWilliams, 2005 ) in the general case of domains discretised with fully-unstructured meshes. These schemes themselves aid the conditioning of pressure solves in geophysical models ( Maddison et al., 2011 ), where there is a large disparity of scales and resolution of the dominant physical processes. This piezometric variable satisfies the same Eq. (20) , with a modified right-hand side source term and boundary condition. Discretisation of the kinematic condition (5) defined in terms of the piezometric pressure using implicit backward Euler in time gives a Robin condition for φ, such that With the barotropic wave speed c = gH , for a distance H , and noting that the ratio of the terms in (22) scale as the square of the time it takes for a barotropic wave to travel a distance H relative to the length of a time step, and we see that the condition for free surface flows (22) tends to that of the rigid lid (21) in the large time step limit. So although adjusting a system to apply a free surface kinematic boundary condition on the top surface as opposed to a rigid lid does improve conditioning for modest aspect ratios, as the disparity in scales increases and the aspect ratio becomes smaller, or equivalently larger time steps are taken, the ill-conditioning of a rigid lid system is soon recovered due to the quadratic dependence.
The multigrid preconditioner of Kramer et al. (2010) for unstructured meshes on high aspect ratio domains helps better condition the Poisson problem in general, without consideration of WD, using a combination of algebraic multigrid and a geometric vertical prolongation operator. This solver method itself makes it feasible to run non-hydrostatic unstructured mesh simulations of fluids in geophysical domains.
Whilst the relatively moderate aspect ratio wet areas can be treated by the multigrid preconditioner approach, specific methods to handle the acute aspect ratio and direct rigid lid condition applied in dry regions are required, if this general fully 3D and non-hydrostatic WD approach is to be applied to real geophysical systems.

Quantification of the ill-conditioning
The discrete form of the Laplacian operator that appears on the left-hand side of (20) , seen in (16) , has eigenvalues λ i ∼ k 2 i , for wavenumbers k i . The conditioning of the matrix is determined by the ratio of the maximum λ ∞ and minimum λ min eigenvalues. This is, equivalently, the ratio of the minimum and maximum wavenumbers, k min and k max , squared For high aspect ratio problems H / L 1, for H and L characteristic length scales of the solution domain in the vertical and horizontal, respectively (see Fig. 1 ), we find On a spheroid, such as the Earth, the characteristic 'horizontal' length scale L is the extent of the encompassing surface geoid, with H the height in a direction parallel to gravitational acceleration. Conditioning of linear system that results from the discretisation of the Poisson equation is approximately proportional to the square of the aspect ratio of the global domain. Equivalently, the element edge-lengths, which are constrained to resolve processes important to the simulation, can also be used to characterise the scaling, such that condition number is proportional to ( x / z ) 2 , with x and z characteristic edge-lengths in local horizontal and vertical directions, respectively. Ideally, entries into the matrix of the linear system that arise from dry cells would be removed, in a process similar to lifted Dirichlet boundary conditions (e.g. Karniadakis and Sherwin, 1999 ) and the solver limited to variables on the wet sub-system, much like an element removal approach. For an implicit approach it is not clear how this would be accomplished without adversely affecting the natural evolution of the interface. Instead, under implicit integration, treatment of the ill-conditioning highlighted by (24) needs to be addressed.

Vertical velocity relaxation in dry areas
To close the system (1) -(2) , an equation of state is required. Details of the form of this function do not influence the development that follows, and a general treatment of the evolution of density is considered, such that with its temporal discretisation following Section 3.1 as Development of the approach proceeds with a discretisation of the density transport Eq. (25) , in a slightly different linearisation to that of (26) , of the form with vertical velocity w , starred variables representing the best current guess, and the source term s n +1 ρ containing details of spatial gradients of density locally aligned to the geoid. Subtracting (28) from (27) gives a transport equation, that mirrors (11) , describing the variation over the Picard iteration process Substitution of this temporally discrete density transport Eq. (29) into the momentum Eq. (7) leads to The FE weak form of (30) is developed by testing with velocity basis functions φ i and applying integration by parts twice at the free surface to obtain The density of air just above the free surface interface ρ a , can in most cases be neglected as a small effect, in the same way as the atmospheric pressure. Assuming that, for shallow waters, the vertical velocity w is linearly related to the distance from the bottom of the ocean or in a depth-averaged sense, and ignoring the density variations ρ in the surface integral above, the terms in (31) above containing explicit reference to the vertical velocity w can be grouped into an absorption term The inverse time scale for the vertical velocity relaxation is defined by (32) . As the Picard iterations proceed, and u n +1 * → u n +1 , the magnitude of this stabilising term, marked by † in (33) , relaxes to zero. Although the absorption coefficient σ will be relatively small in wet regions, and the contribution from † small overall, it is important to include these terms throughout in order to maintain consistency and as a result, accuracy.
The following conditions on vertical density gradient, the free surface, vertical viscosity, and vertical absorption provide a wellconditioned pressure Poisson equation 3. Vertical viscosity where a is a tolerable aspect ratio of element length scales (e.g. unity in the isotropic case), x and z characterise local resolution scales, ν zz is a kinematic viscosity, and σ zz an absorption. Note that the viscosity of (36) must be treated implicitly or semi-implicitly in pressure (e.g. diagonal or block diagonal) in order to control the condition number of the pressure Laplacian. Implementation of the viscosity in stress form is appropriate here since tensor forms directly smooth horizontal velocities in the vertical.
In the case of WD, where (35) does not hold, we must ensure (37) is satisfied by a suitable choice of the absorption σ zz . From (37) , in order to make the resulting pressure matrix feel like an O ( a ) aspect ratio domain, we need The form of σ zz in (38) defines the inverse time scale for the vertical velocity relaxation in (32) . Note that this form of σ zz is spatially varying, and in particular, the characteristic local length scales x and z are non-homogeneous across the geoid surface. In WD simulations, these fields contain large deviations, indicative of the regions affecting conditioning of the pressure Poisson equation.

Discretisation for high aspect ratio domains
The momentum Eq. (33) discretised in space-time at any given Picard iteration step is This is equivalent to (11) previously, noting the term marked † in (33) is zero at the end of a time step.
The new form of the combined kinematic boundary condition (19) leads to a time discretised form modified from (13) to include the no normal flow component applied in dry regions, described by Moreover, the surface integral M s is modified such that the discrete modified kinematic conditon (14) becomes This kinematic condition change modifies the pressure correction, and the discrete continuity (15) becomes The predictor -corrector method of Section 3.5 solves the nonlinear system with the updated p( p , η) Poisson corrector (42) combined with discrete linearised momentum (39) , and a velocity correction determined from (40) .

Self-consistency and physical basis of the solution
Mass, momentum and tracer quantities are self-consistent and conserved, properties inherited from their underlying Galerkin FE weak formulations ( Piggott et al., 2008 ) and use of a thin-film ( Funke et al., 2011;Medeiros and Hagen, 2013 ). The constrained discrete Soblev solution space of the weak form modified with the additional terms marked † in (33) converges on the solution space of the original form without these, as the Picard process proceeds. In a similar manner to Petrov-Galerkin and variational multiscale ( Hughes et al., 1998 ) residual-based stabilisation methods ( Candy, 2008 ), this ensures consistency, that the solution found is a valid solution of the original weak Galerkin formulation, a true discrete solution to the governing continuum equations and is therefore physically-based.
Moreover, just like streamline-upwind Petrov-Galerkin (SUPG) stabilisation, the additional terms themselves are defined from physical properties of the flow. For example, (32) includes contributions from g , the local vertical density gradient and water column depth. This is supported along with local discretisation parameters such as time step and element size used to quantify unresolved scales, in a similar way to multiscale turbulence closures ( Candy, 2008 ).

Determination of characteristic length scales
Accurate calculation of the characteristic length scales is critical to the success of the approach, particularly due to the quadratic dependence in (38) .
The calculation of the characteristic horizontal length scale could be simply the minimum or maximum edge length of the element projected to a 2D horizontal geoid. A more accurate approximation can be determined from the smallest and largest circumscribing circular bounds of this projection. The length scale x is then a function of these minimum and maximum extents. This is a natural approach for models employing anisotropic mesh elements.
The vertical length scale is less ambiguous to determine, since unique intersections with η and b exist ∀ x ∈ , due to the construction 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 further increases accuracy, since FE assembly integrations are performed this way, with options to develop mean or area-weighted means. This is trivially extended to superparametric elements which are typically used in the top layer for accurate representation of geoid curvature.
Arguably the best characterisation of tetrahedral element size is determined from the Jacobian transformation matrix which projects a FE from global to local parameterised space. The determinant of the transformation Jacobian intersected with the local (to quadrature point) surface geoid plane and gravitational acceleration vector will give characteristic length scales for the element in the required horizontal and vertical directions, respectively. This approach also naturally handles element anisotrophy and meshes which are fully unstructured in 3D. The merits of this choice are examined in Section 7.3 .

Correction to velocity relaxation in shallow regions
Under no forcing the momentum Eq. (33) tends to relax the implicit velocity u n +1 * to the state in the previous time step u n , but this can be too strong in very shallow areas. This is corrected by reducing the magnitude of the explicit part of the velocity that we relax to, by adding −γ u n to the right of the momentum equation, with For a water column depth d , this relaxation scales away the velocity in the vicinity of dry regions where d < 2 d 0 , and relaxes to zero in dry regions, where d = d 0 .

Discrete function space updates
The free surface evolution results in many quantities varying in time, such as the free surface normal vector n in (13) . Moreover, this includes the mesh, and hence spatial discretisation, which leads to a change of the discrete function spaces S h , and their spanning basis sets resulting in new forms of mass and other matrices in discrete forms such as (14) . For conservation and accuracy it is necessary to update the discrete non-linear system during the Picard iteration to reflect these changes. There are various techniques to handle this conservatively, through the definition of a grid velocity, for example. In this formulation, the domain discretisation is updated at the end of a Picard iteration to reflect the new free surface height predicted, with the normal n , mass matrix and other matrices representing advection and surface integrals recalculated under the new domain discretisation. It is therefore the case that the discrete matrices M u , M s , C, AG , and Q ; free surface normal n , basis functions φ and ψ, and domain and free surface η , are always the best known approximation, i.e. the starred n + 1 case. A subtle exception is at the end of the final Picard iteration, where the update is not made, to ensure the domain and derivative parameters are those the prognostic variables were calculated on.
The generalised approach that includes the non-linear advection term in the governing Eq. (1) precludes the discretised spatial operator from being self-adjoint. Evaluation of this non-linear term requires sub-cycling, and for under-resolved high Froude number or rapidly-varying flows this could require a large number of iterations to converge, unless the continuum system is linearised, or local resolution increased.

Surface representation and interface tracking
At the end of each Picard iteration, as outlined in Section 3 , the free surface position is updated using (18) to reflect the new pressure at the interfac e p n +1 | ηs . Due to the minimum threshold d 0 , the perturbation of the interface in the direction of the gravitational acceleration is limited. If the pressure p n +1 at the interface implies it should move below this level, it is fixed at the threshold level above the bottom bathymetry (i.e. η = h + d 0 ). The pressure remains unaffected, and is allowed to deviate from the interface position η. Conversely, as soon as p n +1 produces a water column depth greater than d 0 the free surface interface moves upwards. Correspondingly, the domain discretisation is updated with the mesh stretched in the direction −n g , parallel to the gravitational vector, to meet the new free surface bound.
Spatial representation of η is inherited from the function space used to approximate the combined p( p , η) . Irrespective of the order of variation, such as quadratic for the P DG 1 − P 2 element pair, the interface is approximated by a piecewise linear function as far as the domain representation is concerned. This satisfies the minmax property, such that the extent of the surface is bounded by the nodal positions that define its representation. This, together with the minimum threshold level prevents elements from becoming inverted or excessively small.

Remeshing
It is not a requirement that the domain is remeshed anew to these adjusted bounds. Since only one of the domain boundaries is perturbed through the above process and in a direction aligned to the gravitational vector field, locally orthogonal to other bounds of the domain, it is possible to apply a relatively simple r -adaptive transform. The domain mesh is stretched linearly in this direction to fit the new boundary. It is also possible to limit the perturbation to the nodes on the free surface, or to apply more complicated r -or h -adaptive strategies to achieve a hybridised coordinate system (see Kleptsova et al., 2010;Bleck, 2002;Burchard and Petersen, 1997 ) for more accurate solutions or better-represented features. The implementation of the approach described here in the model code (Fluidity, Piggott et al., 2008 ) functions with and supports these methods.

Additional stabilising approaches for dry areas in high aspect ratio domains
Two supplementary approaches to control conditioning are presented, acting directly to prevent strong erroneous flows developing in the thin film and modifying behaviour in neighbouring wet regions. This is exacerbated by the fact the physical system is solved in a weak sense, which whilst better for conditioning can permit large fluxes across the interface. Unlike the above, these approaches transform the solution space and have the potential to affect the solution in unphysical ways. They are presented as additional techniques which can be employed to enable a solution to be reached, but require careful application.

Manning-Strickler drag and dry region stability
In the case of inundation flows where WD is applied, a parameterisation of drag that is commonly employed is the Manning-Strickler formulation, defining the bottom stress n · μ∇ u = n 2 g | u | u 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 formulation itself has a stabilising effect, and more so in the very shallow dry regions, with a drag applied along the bottom boundary proportional to d −1 / 3 . In practice, the Manning-Strickler bottom stress is sufficient to prevent significant erroneous flow developing in dry areas. In the cases of acute high aspect ratio, long time steps or particularly steep bathymetric gradients, the Manning coefficient can be increased in dry regions and their proximity to increase the stabilising effect, with where ˆ n replaces n in (44) , and for n dry a new Manning coefficient (with usual standard units of sm −1 / 3 ) large in size, relative to the standard coefficient n .

Horizontal bulk eddy viscosity in dry regions
A second solution to increase stability, is to damp flow directly in dry regions with a bulk volume viscosity or a source-absorption sponge, both allowing the approach to remain implicit. This stabilisation is applied throughout the domain, or selectively in dry regions and their immediate proximity, with the large horizontal viscosity introduced to control spurious horizontal fluxes, with ν dry a constant eddy viscosity coefficient and d ≥ d 0 ∀ x ∈ . This horizontal viscosity is continuous in space without discontinuous jumps in intensity across the WD interface, acting in the proximity of dry regions where d < 2 d 0 .

Validation and application: Numerical tests
Performance of the implicit WD formulation described in Sections 2 -5 , and additional strategies of Section 6 are examined in four test scenarios in acutely high aspect ratio domains, to a degree found in geophysical systems.

Implementation and verification
The approach has been implemented and validated in the FE fluid dynamics code Fluidity ( Piggott et al., 2008 ). This simulation framework contains many tools for geophysical modelling, is parallelised with sophisticated load balancing and supports adaptive mesh methods allowing computational effort to be focused on regions of dynamic interest. It functions for a spatially variable gravitational acceleration vector, and hence can be used for large-scale simulations on the Earth's spheroid. The implementation includes a suite of test cases to routinely verify the new algorithm in a formal sense, in an automated continuous verification build engine ( Farrell et al., 2011 ) to ensure robustness of the code and resiliency in light of further development. The unstructured meshes used in the following cases were built by means of the open source software Gmsh 1 .
The balance and LBB stability properties of the P DG 1 − P 2 velocity-pressure pairing (see Cotter et al., 2009 ) aid conditioning and are used in all applications considered here. All four cases have been run on the purely continuous pairing P 1 − P 1 also, but due to the pressure filtering required, did not perform as well, and in all but modest aspect ratio cases were too ill-conditioned to reach convergence. The behaviour of P DG 1 − P 2 and P 1 − P 1 and their relative performance in regular aspect ratio problems is presented in Cotter et al. (2009) .
Due to the aspect ratios considered, all cases use the multigrid preconditioner described in Kramer et al. (2010) for iterative solution of the conditioned symmetric pressure Poisson linear system in combination with Conjugate Gradient (CG, Hestenes and Stiefel, 1952 ). The momentum system is solved in a more standard approach with Symmetric Successive Over-Relaxation (SSOR, Young, 1971 ) preconditioning and the iterative Restarted Generalised Minimal Residual (GMRES, Saad and Schultz, 1986 ) algorithm, where the calculation is restarted after k = 30 iterations. The iterative SSOR-GMRES process is performed using algorithms built into the established and well-verified PETSc library ( Balay et al., 1997 ). In contrast to the study ( Funke et al., 2011 ), it was found that two Picard iterations provide sufficient convergence of the coupled system in the cases studied. In all cases, both linear systems are solved to a convergence criteria specified by a relative error tolerance of 10 −7 , which is considered sufficiently accurate. The quadrature based subgrid resolution described in Funke et al. (2011) is also used, with a quadrature degree of eight.

First Balzano sloped channel benchmark
The first two sets of numerical tests are from the suite of problems in Balzano (1998) , selected since they exhibit the problematic ill-conditioning in as simple a setup as possible. No analytical solution is available, so the problem configuration is chosen consistently with Balzano (1998) to be able to draw comparisons. The base benchmark case is developed from the originally 2D domain consisting of a 13.8km long slope with a depth of 5m at one end which tends to zero at the other. Recently developed schemes, such as the flux-limiting WD method for FE SWE models presented in Gourgue et al. (2009) and the non-hydrostatic algorithm proposed in Funke et al. (2011) , have been benchmarked on these cases. These model in 3D, but force dynamics to occur predominantly in the directions where the extremes in extent occur, with 10 elements introduced in the third direction in the former and 1-2 in the latter, which is followed here to a width of 1km. With the assumption solutions are laminar, this extrusion into 3D space will not change the physical behaviour. The sloped bottom bathymetry is defined h (x ) = x/ 2760 , for the x -coordinate direction indicated alongside the surface geoid computational mesh in Fig. 2 (a). The base case single-layer mesh contains vertically-aligned nodes and a horizontal element size of 500m.
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 . 02 sm −1 / 3 is applied at the bottom boundary. The gravitational acceleration is set to 9 . 81ms −2 and the fluid is initially at rest. Time discretisation is performed with Crank-Nicholson integration (i.e. θ = 1 2 ) and a time step of 600s. In this case the WD threshold is set at d 0 = 0 . 5 mm . The free surface is forced at the deep open boundary with a sinusoidal variation of amplitude 2m, such that water column thickness oscillates between 3-7m, with a period of 12h.
In the series of tests considered here, the horizontal extent is varied from 1.38 × 10 2 m to 1.38 × 10 6 m, centred about the defined benchmark length of 1.38 × 10 4 m. This provides a range of element aspect ratios from 10 −4 to 10 −8 , a domain aspect ratio up to 3 . 62 × 10 −6 and spatial scales spanning over 10 orders of magnitude 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 direction also scaled to preserve element shape. The time step is also scaled to ensure the wave Courant number is constant. The WD threshold d 0 , and vertical extent are kept constant across all cases. The free surface evolution of the intermediate case with a horizontal extent of 1.38 × 10 4 m is shown in Fig. 2 at 10min intervals, matching ( Balzano, 1998 ) and ( Gourgue et al., 2009 ), for the initial drying and then wetting phase, respectively. The results are physically reasonable, and comparable to other formulations ( Funke et al., 2011 andGourgue et al., 2009 for example). In particular, the free surface interface does not suffer from either underestimation with negative water column thickness, nor produce oscillations during the wetting process observed in Balzano (1998) for some of the 10 methods examined. This behaviour is characteristic of the solutions across the range of aspect ratios.
Through a modification of the optimum aspect ratio parameter a in (38) , there is a corresponding change in the aspect ratio felt in the discrete pressure matrix of elements in dry regions. The parameter a is varied over the range a ∈ [ 10 −4 , 10 4 ] in a suite of 1001 simulations of the base Balzano case. Solver iteration number is used as an indicator of conditioning, and plotted in Fig. 3 for both the pressure and velocity calculations as mean and maximum values over the course of a WD phase. The parameter range has been spaced equally in log-space in order to give a good representation of the behaviour over the large range of domain aspect ratios. This is achieved with a discrete parameter space defined for Whilst the conditioning of the velocity solver is largely unaffected, the number of iterations required for pressure convergence increases dramatically as the magnitude of the parameter a increases. As the aspect ratio parameter becomes acutely large with | a | → ∞ , behaviour tends to that of the system without the scheme applied. It is clear that the vertical relaxation scheme has a positive impact on conditioning, reducing the number of required iterations in the pressure solution in this test by a factor of 20. With Picard iteration numbers also reduced as a consequence, this effect is multiplied for overall performance.
Changes in the parameter demonstrate the scheme significantly improves conditioning in the base case. Now an optimal aspect ratio a = 1 is specified and actual changes to the domain extents considered. Again a suite of simulations are run to span the parameter space and determine conditioning, and the number of iterations required for convergence of the pressure is shown in Fig. 4 , with and without the relaxation consitioning. In the range considered, the improvement is reduced by a factor of up to 20 and results highlight the approach eliminates a dependence of conditioning on aspect ratio.

Second Balzano shelf channel benchmark
This case also originates in Balzano (1998) and differs from the first by the inclusion of a shelf break in the sloped bathymetry, defined in Appendix C . The horizontal domain is discretised in a way to ensure accurate bottom boundary representation, ensuring element faces align with the discontinuous changes in surface gradient ( Fig. 5 ). Except for the change in bathymetry, discretisation proceeds in the same manner as the first Balzano case of Section 7.2 , and is again run over a range of aspect ratios.
The free surface evolution in the case with minimum element aspect rati o 10 −6 is shown in Fig. 5 , again characteristic of the formulation over the range of aspect ratios. In addition to the oscillatory and retention problems already mentioned, Balzano noticed a runoff problem with some methods in this test case, where water remains on the shelf during the dry period instead of flowing Table 1 Pressure solver iterations to convergence in the second Balzano benchmark with a domain aspect ratio of 3 . 62 × 10 −6 for five approaches to calculating characteristic height.
With the irregular bathymetry of this case, we consider the effect of how the length scales that are passed to the relaxation scheme are calculated, as discussed in Section 4.7 . The characteristic height z varies both by element and over elements, and can be calculated at quadrature points for increased accuracy. Noting the role of these length scales in the vertical velocity relaxation inverse time scale (38) , we see that errors in how they are determined influence conditioning in the same manner as that of perturbations of a from the optimum value of 1, except to a greater degree due to the quadratic dependence which, following Fig. 3 , reduces the effectiveness of the conditioning.
Five approaches are considered ( Table 1 and Section 4.7 ). The methods 'minimum', 'maximum' and 'mean' each refer to the minimum, maximum and mean of the set of 6 vertical lengths z calculated from the 4 tetrahedral element vertices. The minimum of these performs poorly in all phases, so its value was limited by a lower bound in the 'minimum capped' approach, which prevents the applied absorption becoming too large. This produced better conditioning than the maximum in the drying phase, and whilst improved in the wetting phase, the maximum here still produced better conditioning. The mean behaves very well in the drying phase, but only satisfactorily during wetting. This implies all three of these norms are not capturing all of the important parameters to determine an optimum z . The Jacobian approach using the determinant of a contracted transformation matrix at quadrature points provides the best conditioning during the wetting phase. The conditioning in the drying phase is not consistently the best, but the lowest mean number implies it is best overall. In the Balzano shelf case examined here, the number of iterations required for convergence is approximately halved by a careful consideration of the calculation of z .

Thacker parabolic basin benchmark
The Thacker parabolic bowl ( Thacker, 1981 ) is an idealised ocean basin that thins at its edges, with bathymetry defined in Appendix D . It is a challenging free surface flow problem with WD that has previously been used in intercomparison studies ( Balzano, 1998;Funke et al., 2011;Gourgue et al., 2009;Kärnä et al., 2011 ). An analytical solution for the evolution of the free surface is known (also in Appendix D ) when both dissipation and Coriolis are absent, and the case suitable for the evaluation of spatial and temporal accuracy, and volume conservation.
The base case domain size matches that of Balzano (1998)  Conditioning is examined for domain aspect ratios ranging over four orders of magnitude, from the base 5 . 68 × 10 −5 down to 5 . 68 × 10 −8 . This is achieved through vertical scaling the domain and d 0 , with the maximum equilibrium water column depth varying between 50 m − 5 cm . With the characteristic horizontal edge length ∼ 10 4 m close to the edges where the domain drys, element aspect ratios vary similarly ∼ 5 × 10 −5 − 5 × 10 −8 . A cross section of the resulting single-layer basin domain for the h 0 = 5 cm case is shown in Fig. 6 (b), with the initial perturbation η 0 ensuring a minimum thickness of d 0 is applied.
Edge element length scales are defined isotropically by which for the case x = 10 4 m in Fig. 6 , result in a range from 100km in the middle down to 10km at a distance 3.8 × 10 5 m from the centre, in an approach following ( Funke et al., 2011 ). Numerical evolution of the free surface for the highest aspect ratio case, shown in Fig. 7 (a)-(b), is observed to fit the analytical solution very well, even with elements of a very high aspect ratio ( 5 × 10 −8 ). Like the results from the more modest domain size a phase shift is observed, which also seen in Funke et al. (2011) , is a feature produced by the thin layer in the dry areas. The contribution of numerical dissipation inherent in the scheme is eliminated, We can eliminate numerical dissipation inherent in the scheme as a contributor, as we find that with solves iterated to convergence, volume is conserved up to a relative factor of 1 . 0 × 10 −11 , which can be attributed to numerical round off error. This phase shift is reduced with an increase in mesh resolution (see Funke et al., 2011 ), which contributes to the increase in accuracy observed in Fig. 9 (b). In the time series taken at the edge of the domain, it is clear when the location becomes dry in both the analytical and numerical solution, and where the factor of d 0 is maintained in the latter (here 0.5mm).
The radial velocity at the free surface at two locations is presented in Fig. 7 (c)-(d) at approximately the same relative locations as those considered in Casulli and Zanolli (2007) , and is compared to the analytical solution provided in Appendix D . In the main body of fluid the solution is a very good match, with the same shift observed in η as in Fig. 8 above and ( Funke et al., 2011 ). Close to the edge of the basin, u r is not as well predicted as η. This is partly due to the continuous nature of the thin-film ap-  proach, which solves for u in both wet and dry regions. The spatial discretisation local to this point is relatively coarse, and additionally is not aligned to a radial direction, which makes u r particularly challenging to calculate. This and the phase error, can be mitigated by increasing spatial resolution and constraining mesh structure to align with flow direction in inundation regions. Importantly, accuracy of Funke et al. (2011) is maintained, whilst the difficulty in solving the linear systems is much reduced. The position of the free surface in a vertical slice of the domain along the line indicated in Fig. 6 and after a period of thirty days, to include two each of the WD phases, is shown in Fig. 8 . Spatially, the numerical solution is a good fit to the analytical solution and its resolution of the WD front comparable to studies in more modest aspect ratio domains ( Funke et al., 2011;Gourgue et al., 2009 ). The use of the vertical velocity relaxation approach and iterative solvers for the linear systems does not have a significant impact on accuracy of the solution, and provides a formulation for high aspect ratio domains that performs as well as those of modest size.
An evaluation of error e ( t ), at a time t , is made under an L 2 norm of the absolute difference, such that for h and η a defined in Appendix D . The minimum water depth is included in the analytical solution, since this is the free surface height the formulation converges to, and the domain encompasses both wet and dry regions.
Solution convergence with respect to the smallest horizontal characteristic edge length x is considered in Fig. 9 (a) for the base domain, where the time step is linearly scaled to maintain a constant CFL number. Meshed domains are generated by scaling the metric (46) . With this WD formulation, we obtain the linear convergence in error to characteristic edge length observed in Funke et al. (2011) .
The impact of domain aspect ratio on the accuracy of the calculation of free surface height after the initial wetting phase is considered in Fig. 9 (b). Notably the error does not increase significantly with an increase in the magnitude of the aspect ratio, and is far from linear. The increase could be accounted for, to some extent, by the fixed relative tolerance on the iterative solvers of the linear systems. Adjusting this tolerance to increase convergence in Fig. 9. Thacker benchmark convergence properties. Accuracy relative to (a) edge length in base domain of Fig. 6 and with respect to aspect ratio in (b). Error is evaluated at the point of time that the initial wetting period is complete. Linear and quadratic gradients are indicated by dashed lines. The diagonal cross, marked by * , points to a case with a relative tolerance reduced to 10 −10 .
cases with very small edge lengths could help to increase accuracy at this level. A small improvement in accuracy is seen in the highest aspect ratio case considered in Fig. 9 (b) where the relative error tolerance of 10 −7 described in 7 is reduced to 10 −10 . It is a significant result that a solution can be found for these cases with very high aspect ratios and additionally, that the approach does not have an appreciable impact on accuracy.

Basin inundation
This case considers the inundation of water into an initially dry basin, with the effect of bathymetric features on WD front propagation also examined. The base domain is shown in Fig. 10 and consists of a basin with horizontal extent 100m × 100m, and an inlet of width 10m, its centre positioned 15m in from one of the corners. The domain is discretised with elements of a characteristic edge length of 5m. The problem is forced with a normal inlet velocity of 0 . 5 ms −1 to model a levee breach into a flood plain on an urban scale.
To provide a more natural forcing, instead of applying a flux directly on the boundary, the inlet is extended back 10m and is maintained wet throughout by developing a sloped bathymetry back, down to a depth of b = 20 m , as seen Fig. 10 (b). The normal inlet velocity is then applied to the face that has been extended back, with velocity slip conditions on the adjacent sides. This was found to avoid problems with the inflow at the edges of the breach. At the outflow on the far boundary at y = 100 m , a natural Neumann condition is applied perpendicular to the boundary, such that ∂v / ∂y = 0 , for velocity v in the y −direction. All other boundaries are closed, with no normal flow conditions. Other velocity components are free and left unconstrained.
To ensure accuracy of the calculation of prognostic variables is not affected by the use of relatively large time steps with potential impact on the conditioning analysis, t is set conservatively at 10s to give a maximum Courant number of 1.
In a similar approach taken for the Balzano slope case, we consider the influence of the optimum aspect ratio parameter a on conditioning in the base domain with aspect ratio 10 4 , over a parameter space spanned by 1001 simulations shown in Fig. 11 . Conditioning of the pressure solver is significantly increased, by over a factor of six in this modest aspect ratio case. Again velocity is only slightly affected, and felt through the coupling, a consequence of better pressure conditioning. When varying simulation domain extent, with an optimal choice of a = 1 , similar behaviour is observed and shown in Fig. 12 . In practice, large gradients in bathymetry have an impact on conditioning. This is studied with the introduction of a depression in the domain to form a hollow, and conversely, a raised hill. Both interact differently with the incoming wetting front. These features Fig. 11. Impact of the optimal aspect ratio parameter on solver conditioning in the flood plain basin benchmark over a WD phase. 1001 individual simulations shown. are introduced to the domain with a Gaussian perturbation, which is defined at all points r on the horizontal surface on the do- , for h 0 the maximum deviation in height, which occurs at the centre where r = r 0 . The inverse variance vector σ defines width, and consequently the gradient, of the obstacle. In the scales of the base case, the magnitude of the perturbation | h 0 | is 5m, with a width of 10m, defined by σ = 1 10 , 1 10 . The perturbation is positioned at 30m in from each of the bounding edges at the corner closest to the inlet. At the start, the minimum water thickness of d 0 = 1 cm is applied above the bathymetry, to provide the initial thin dry flood basin.
The above is used to generate an inundation into a domain containing a large hollow with h 0 = −5 m . Conditioning is further decreased with the presence of the hollow, with a mean number of 105 iterations required in pressure for the modest aspect ratio case. Compared to the flat case, the number of iterations required increases at a greater rate, and the positive effect on conditioning of the vertical relaxation scheme is further pronounced. Additionally, the large gradients in bathymetry adversely affect conditioning of the velocity solver early in the simulation where large velocities develop around the steep slopes to fill the hollow. This can be seen in the example snapshot results shown in Fig. 14 . Initially flow is strong from the breach, and predominantly flows into the hollow, whose surface oscillates in a similar manner to that seen in the Thacker parabolic bowl benchmark of Section 7.4 . Once the hollow is filled, the free surface peaks and a hydraulic jump develops between the fast-flowing inlet from the breach and the formed lake. The fluid then gains momentum in the direction of the inlet flow and is seen to build up on the opposite boundary. A clear front has developed by this stage, and begins to propagates across the plain towards the open boundary. It is also possible to see the larger velocities that develop at the front and ahead in the thin dry regions. This is motivation for the application of velocity conditioning discussed in the following.
In the case of the hill, with h 0 = 5 m in the base domain, the effect on velocity is more significant, particularly as the WD front meets the bathymetric intrusion. In this case it is necessary to apply a regularisation to the momentum equation to improve conditioning, which is achieved through an application of a bulk volume viscosity, as introduced in Section 6.2 . The domain-wide horizontal viscosity ν L , is varied over the range ν L ∈ [ 10 −4 , 10 4 ] m 2 s −1 through (45) in the mid aspect ratio 10 −6 case, with conditioning shown in Fig. 13 . As a general trend, the number of iterations required increases with strengthening of the horizontal viscosity. There is however a point at which there is a noticeable dip, where the increase in intensity improves conditioning. This decrease in the mean number of iterations is due to improvement of conditioning made when the front approaches and traverses the hill protrusion. Limiting application of this conditioning to dry regions and its proximity, as described in Section 6.2 , allows WD fronts to encounter steep changes in bathymetry without the corresponding impact on conditioning of the velocity solver, in this implicit and continuous WD formulation.
Increasing the bottom drag through the Manning-Strickler parameterisation in this region as outlined in Section 6.1 also acts to improve conditioning. In particularly high aspect ratio cases with steep bathymetry, the velocity solver is too badly conditioned for efficient solution with a SSOR-GMRES iterative process without approaches such as the horizontal viscosity and drag discussed.

Conclusions
In this paper a novel approach to efficiently modelling WD inundation processes in 3D, capturing non-hydrostatic and baroclinic 0 Fig. 14. Inundation into a flood basin of side length 100m and threshold value 1cm containing a hollow bathymetric feature. Three successive visualisations with (a) the hollow filling, (b) a hydraulic jump and (c) propagation further into the plain, are shown at 7,350, 20,400 and 37,650s into the simulation, respectively. The left panels contain contour plots of free surface perturbation, overlaid with magnitude-scaled vectors of depth-integrated velocity. The right panels present the 3D domain stretched in the vertical by a factor of 40, to better show the change in free surface height, with the inlet breach and connecting reservoir seen on the right side. Contours of the magnitude of surface velocity are plotted together with vectors indicating the surface flow direction. Note the velocity fields presented are those in a continuous P 1 space, calculated through a Galerkin projection from the discontinuous P DG 1 prognostic velocity field (see Candy, 2008 ). physics, in the high aspect ratio domains that characterise geophysical systems has been proposed.
This has identified the ill-conditioning present in implicit continuum WD methods applied in fully 3D fluid flow models. Following a quantification of the highly spatial and temporally variable contributing factors, regularisation of the governing weak form leads to a linear system that appears as a unit aspect ratio problem. The result is that the approach can be used to model WD in multiscale geophysical domains, seamlessly alongside other challenging physics, such as baroclinic and non-hydrostatic flow, without a severe and limiting impact on the iterative solvers typically required for efficient simulation of multi-physics 3D dynamics.
The approach has been demonstrated effective over a wide range of spatial scales and correspondingly, aspect ratios. The predicted behaviour on convergence is verified in numerical tests in both domain and element aspect ratios representing up to 8 orders of magnitude difference, with discrete domains containing spatial scales spanning 10 orders of magnitude.
The approach imposes no restrictions on space and time discretisation, permitting an arbitrarily flexible mesh choice (including generalised vertical coordinates), order of representation and implicit time integration. All are important for system models simulating over a range of scales and physics. Discretisation can be chosen largely independent of WD considerations, with for example, spatial resolution focused on local physics modelling demands.
Use of a combined p( p , η) variable strictly enforces consistency between the full 3D pressure and free surface perturbation. Notably there is no need to interpolate η and its derivatives from η to the internal domain for inclusion in the momentum calculation. Consistency with other fields and conservation are achieved by the overall FE approach, which can provide a high order continuum representation. P 1 − P 1 and the heterogeneous element pairing P DG 1 − P 2 have been applied in the numerical tests. The implicit treatment of p( p , η) is inherited by p and η, and as a result, t may be based solely on accuracy considerations and not stability when considering free surface wave propagation. As discussed in D'Alpaos and Defina (2007) this may need careful consideration when a system is under-resolved with a relatively irregular bottom topography containing sharp gradients, or in high Froude number rapidly varying flow. A limitation to note is that the free surface interface cannot become unduly complicated, including folds, since the function η is by definition injective with only a single point permitted to lie on the surface for any point within the domain. As such it is not possible to model breaking waves, a common limitation to all of the Eulerian approaches cited. Unlike schemes applying additional viscosity or bed friction based on empirical numerical measures that potentially lead to stabilisation through unphysical means, the approach ensures physical consistency such that resultant solutions are enforced to exist in the space of solutions available to the original physically based weak form of the continuum governing Eqs. (1) -(2) . Physical consistency is verified in the numerical tests. Lastly, since the terms introduced specifically to improve conditioning are formulated in the continuum primitive form, this part of the approach could equally be applied in other WD implementations for an arbitrary underlying discretisation.
This approach will not be optimum for some WD problems, particularly due to the computational cost even with the aspect ratio problem solved, where a single layer SWE approximation is sufficient, or computational efficiency may demand lower order methods for real-time tsunami prediction, for example. However, this approach now enables the modelling of physical phenomena not possible previously, particularly those at the interfaces of traditionally separate fields, and with rapid ongoing develop-ment of computational resources, this approach and similar will grow in use and become more common practice -a way to bring WD to seamless massive multiscale multi-physics Earth system models.