Semi‐implicit subgrid modelling of three‐dimensional free‐surface flows

In this paper, a semi‐implicit numerical model for two‐ and three‐dimensional free‐surface flows will be formulated in such a fashion as to intrinsically account for subgrid bathymetric details. It will be shown that with the proposed subgrid approach the model accuracy can be substantially improved without increasing the corresponding computational effort. Copyright © 2010 John Wiley & Sons, Ltd.


INTRODUCTION
Numerical modelling of free-surface hydrodynamics for environmental applications has been a challenging task for several decades. The ever increasing computer power has stimulated researchers to formulate stable and efficient numerical models for accurate two-and three-dimensional flow simulations.
To bypass the severe restriction imposed on the time step size by the Courant-Friedrichs-Lewy stability condition, an alternating direction implicit method was proposed by Leendertse in 1970 [1]. This successful method has been extensively studied by several authors (see, e.g. [2] and the numerous references contained therein).
A detailed characteristic analysis on the two-dimensional governing equations was carried out in [3] to determine the precise terms that have to be discretized implicitly in order to obtain a semi-implicit method that is stable at a minimal computational cost. The resulting algorithm, which is relatively simple and robust, was soon extended to three-dimensional problems [4] and was employed in several applications. This method has been further extended to apply on orthogonal unstructured grids, which provide higher flexibility for fitting boundaries and for local mesh refinements [5].
Nowadays data of land heights and water depths can be collected with a very high resolution. Airborne scanning laser altimetry, for instance, can produce maps of surface height over large areas with a height precision of about ±15 cm and a spatial resolution of 1 m (see, e.g. [6]). Consequently, data bases with bathymetric data are built with a remarkably fine resolution. 442 V. CASULLI AND G. S. STELLING The availability of superior computing power, solid and stable semi-implicit methods, structured or flexible unstructured grids are useful tools but, alone, are still insufficient to faithfully account for complex topographic features. The availability of detailed bathymetric data within a coarser grid model can and should be used to further improve a model accuracy.
Accounting for subgrid-scale variations in numerical models is a longstanding problem that has been the subject of inquiry for at least four decades. Some of the earlier work is related to averaging for LES. This approach has been extensively pursued in atmospheric boundary layer modelling. A recent review on volume averaging can be found in References [7,8].
Within this context, the method proposed here is a simplified form that takes into account celerity but no new terms arising from body forces. As recently shown in Reference [9], the use of subgrid bathymetric details permits efficient and very accurate simulation of the wetting and drying dynamic with a simple, semi-implicit two-dimensional shallow water model.
In the present investigation, the subgrid approach presented in [9] is intrinsically extended to three-dimensional flow modelling. Subgrid bathymetric resolution is incorporated into the semiimplicit numerical method described in [5] in order to account for detailed bottom location and local bottom shear stress, thus permitting a substantial improvement of model accuracy without increasing the computational effort.

Governing differential equations
For convenience, unless otherwise stated, the notation used here is taken from References [5,9]. Thus, the momentum equations and the incompressibility condition for hydrostatic flows are taken to be where u, v and w are the velocity components in the horizontal x-, y-and vertical z-directions, respectively; t is the time; is the water surface elevation; f is the Coriolis parameter; g is the gravitational acceleration; and h and v are prescribed non-negative coefficients of horizontal and vertical eddy viscosity, respectively. Integrating the continuity equation over depth, and using a kinematic condition at the freesurface, leads to the following free-surface equation: where h is the specified bathymetry and H = max(0, h + ) is the total water depth. The boundary condition at the free-surface and at the sediment-water interfaces are assumed to be where is a non-negative bottom friction coefficient and u * and v * are the velocity components near the bottom. When wetting or drying is expected, the differential equations (1)-(4) are defined on a timedependent horizontal domain (t) defined as (t) = {(x, y) : H (x, y, t)>0} [9].

Unstructured orthogonal grid
To solve Equations (1)-(4) numerically, a fixed domain is chosen in such a fashion that (t) ⊆ for all t 0.
is then covered by an unstructured orthogonal grid [5] consisting of a set of 443 non-overlapping convex polygons i , i = 1, 2, . . . , N p , separated by N s sides j , j = 1, 2, . . . , N s . Within each polygon, a center must be identified in such a way that the segment joining the centers of two adjacent polygons and the side shared by the two polygons have a non-empty intersection and are orthogonal to each other. Once has been covered with an unstructured orthogonal grid, each polygon i may have an arbitrary number of sides. Let S i denote the set of sides of the ith polygon. The left and the right polygons that share the jth internal side are identified by the indices ( j) and r ( j), respectively. Moreover, ℘(i, j) denotes the neighbor of polygon i that shares side j with the ith polygon. The non-zero distance between the centers of two adjacent polygons that share the jth internal side is denoted with j .
Along the vertical direction a simple finite difference discretization, not necessarily uniform, is adopted. By denoting with z k+1/2 a given level surface, the vertical discretization step is defined by The three-dimensional space discretization consists of bathymetry shaved volumes contained in prisms whose horizontal faces are the polygons at two consecutive level surfaces, and whose height is z k . The discrete water surface elevation n i , assumed to be constant within each polygon, is located at the center of the ith polygon.

Semi-implicit finite volume discretization
In the present derivation, rather than assuming constant depth, constant velocity, and constant friction coefficient along each edge of the grid, a spatial variability is allowed. Hence, in a properly oriented co-ordinate system where the jth edge is aligned with the y axis, the local horizontal velocity components perpendicular to the jth edge and at the kth vertical layer will be denoted with u n j,k (y). Thus, assuming a constant pressure gradient along each edge of the computational grid, a detailed momentum balance, previously given by Equation (7) in [5], is now modified to account for spatial variability: where z n j,k (y) is, in general, the distance between two consecutive level surfaces, unless the bottom or the free-surface comes first; G n j,k (y) is a finite difference operator that accounts for the explicit contributions from the discretized Coriolis, advection, horizontal viscosity, and hydrostatic pressure; t is the time step; is an implicitness factor to be taken in the range ∈ [ 1 2 , 1] [10]; and m j (y) and M n j indicate the bottom and the top layers, respectively. A formulation for G n j,k (y) can be specified in several ways, such as by using a stable Eulerian-Lagrangian scheme [5,9] or a conservative formulation that allows for accurate simulation of rapidly varied flows at every Froude number [11,12]. When an Eulerian-Lagrangian discretization is used, a form for G n j,k (y) is taken to be G n j,k (y) = z n j,k (y) where u L j,k (y) denotes the velocity component normal to the jth side of the grid and v L j,k (y) is the tangential velocity component in a right-hand coordinate system [13]. Both components are 444 V. CASULLI AND G. S. STELLING interpolated at time t n at the end of the Lagrangian trajectory based on the values at adjacent grid points. The Lagrangian trajectory is calculated by integrating the velocity backwards in time from the y position along the jth edge at t n+1 to its location at time t n . h is the explicit discretization of the horizontal viscosity terms.
The values of u n+1 j,k (y) above the free-surface and below the bottom in Equation (6) are eliminated by means of the boundary conditions (5), which yield the following difference formulae: After including the boundary conditions (8), Equations (6) can be written in compact matrix form as where A n j (y) is a symmetric, positive definite, tridiagonal matrix that includes bottom friction and vertical viscosity terms; U n+1 j (y) is a vector containing the unknown velocities u n+1 j,k (y); DZ n j (y) is a vector whose entries are the vertical increments z n j,k (y); and G n j (y) is a vector containing the known explicit terms G n j,k (y). Of course, Equations (9) are defined on wet points only, i.e. where DZ n j (y)>0. On dry points U n+1 j (y) = 0 is assumed. On wet points, Equations (9) constitute a set of M n j −m j (y)+1 linear equations that involves the unknowns U n+1 j (y), n+1 ( j) , and n+1 r ( j) . Hence, the new free-surface elevation needs to be computed in order to determine the new velocity field.
Having assumed a constant free-surface elevation within each polygon, upon integration of the free-surface equation (4) over i , one obtains the following semi-implicit finite volume approximation: where is a nonlinear function representing the water volume within the ith water column for a specified water level n i and for any given bathymeric detail; and s i, j is a sign function associated with the orientation of the normal velocity defined on the jth edge [9].
It is worth noting that, other than assuming a constant pressure gradient along each edge of the computational grid, Equations (6)-(8) correctly represent the interaction of bottom friction with the vertical viscosity for any bathymetric distribution. Likewise, other than assuming a constant water level within each computational polygon, Equation (10) represents a precise mass balance over each computational cell and for any subgrid specification (see [9], Section 4.3, for further details).

Solution algorithm
Substitution of U n+1 j (y) from (9) into (10) yields the following discrete wave equation: where n j and b n i are scalar quantities given by Equations (12) can be assembled into a sparse, mildly nonlinear system of N p equations for n+1 i , i = 1, 2, . . ., N p . Details of an efficient iterative solver for this system have been given in Reference [9].
Once the new free-surface elevations n+1 i have been computed, the local velocities along each edge of the computational grid are readily derived from Equation (9).
Finally, the new vertical component of the velocity is diagnostically derived from a finite volume discretization of the incompressibility condition (3).
In practical applications where the bathymetric data are specified as a piecewise constant function over sub-polygons and over sub-edges, respectively, the integral in (11) might be replaced by a summation over sub-polygons, whereas the integrals in (13)-(14) might be replaced by a summation over sub-edges. Figure 1 shows a simple grid where each edge has 4 sub-edges and each polygon has 16 sub-polygons.
It should be noted that if only one vertical layer is specified, then m j (y) = M n j = 1, z n j,1 (y) = H n j (y) and the present algorithm reduces to a semi-implicit scheme that is consistent with the two-dimensional, vertically integrated shallow water equations [9].
Another interesting aspect of the present formulation is given by the fact that the size and the structure of the mildly nonlinear system (12) is independent of the vertical resolution and of the prescribed subgrid details. Subgrid and vertical resolution affect the assembly of Equation (12) and the number of equations for horizontal velocities. Moreover, since a major part of the computational effort is required to determine the free-surface elevation from Equations (12), a fine subgrid and vertical resolution can be adopted with an acceptable increase of the corresponding computing time.
In essence, a constant pressure gradient is assumed on each edge but velocity, friction, and depth are variable. This assumption not only allows to reduce the number of surface elevation points but also removes all the subgrid motions. Subgrid dynamic and near critical or supercritical flows where the pressure gradient variability cannot be neglected can be resolved by refining the computational grid.
Finally, to be noted is that the possible presence of dry sub-edges or dry sub-polygons has no influence on the computed velocities and water levels. Thus, when the subgrid approach is adopted, boundary fitting is naturally obtained, also in the presence of wetting and drying, without resorting to complex grid adaptation strategies [9].
In summary, the resulting algorithm described by Equations (9) and (12) is almost as simple as (13) and (15) from Reference [5]. The only, but important difference being that subgrid resolution greatly enhances the model accuracy.

NUMERICAL TEST
The above method has been tested on a number of realistic field-scale problems where bottom shear stress and wetting and drying are extremely important. The Venice Lagoon, for example, is a very complex system whose area covers about 420 km 2 and consists of several interconnected narrow channels with a maximum width of 1 km and up to 50 m deep encircling large tidal flats. The average water depth is less than 1.5 m. The Lagoon is connected to the Adriatic Sea through three narrow inlets, namely Lido, Malamocco, and Chioggia. Tides propagate from the Adriatic Sea into the Lagoon through the inlets. The historical city of Venice is located in the Lagoon near the Lido inlet and is constantly threatened by storm surges, which often inundate the city.
The bathymetric data, available on a fine regular grid with 25 m resolution, appears to represent the tree-like structure of the main channels with sufficient accuracy (see Figure 2). With such a fine grid, the entire Lagoon can be covered by N p = 671 030 square polygons separated by N s = 1 361 331 edges. The flow is driven by an idealized M 2 tide of 0.5 m amplitude and 12 h period that is specified as a boundary condition at the three inlets. By using a time step t = 300 s and a vertical resolution of 1 m, the three-dimensional flow is simulated for 10 tidal cycles and the corresponding results were assumed as a reference solution.
In a first simulation, four runs have been completed by using a regular horizontal grid of 25, 50, 100 and 300 m resolutions without subgrid refinement. As the grid is coarsened, the corresponding bathymetric data are taken as the average depth computed from the fine grid. Figure 3 shows the computed water volume for three consecutive tidal cycles. As expected, with coarser meshes the tide decays as it propagates to remote shallow areas. This is clearly seen by the computed surface elevation (Figure 4) observed at a remote station located in the Northern part of the Lagoon indicated by a mark in Figure 2.  Surface elevation (m) The previous simulations were repeated by using the available bathymetric data with 25 m resolution to provide subgrid details to coarser grid simulations. The resulting volumes and surface elevations shown in Figures 5 and 6 indicate that the proposed algorithm is much less sensitive to the grid size.
The above calculation were performed on a PC with a 2.16 GHz Core 2 Duo processor using an OpenMP implementation. Table I summarizes the problem size and the model performance. As indicated, a considerable saving in computational effort can be achieved by coarsening the grid while retaining subgrid details.
It should be noted that, without subgrid specifications, accurate bathymetry fitting can only be accomplished by using very fine meshes and/or highly distorted unstructured grids. Extreme mesh refinement and excessive grid distortion, however, lead to large and often ill-conditioned nonlinear systems (12)    discretization also requires an adequately small time step size (see. e.g. [11]) that contributes to an increase in the computational cost.
Other than keeping good accuracy at a drastically reduced computational cost, the present subgrid approach appears to be a very promising tool for applications to three-dimensional morphological models [14] and to multiscale numerical modelling. Additionally, the scalar transport algorithm described in [15] can be readily adapted to the present formulation so that mass conservation and maximum principle of any transported constituent can be assured.

CONCLUSIONS
A semi-implicit algorithm for three-dimensional free-surface flows has been generalized to intrinsically account for detailed bathymetric data at subgrid level.
Accurate mass and momentum balance in shallow flows over complex geometries is assured also in the presence of wetting and drying. The resulting algorithm is numerically stable, relatively simple, and extremely efficient.
The given example demonstrates that the present formulation can generate very accurate results, even with a coarse and structured finite difference mesh, without resorting to costly and unnecessary grid refinements.
The interested reader should readily see how to apply the present subgrid concept to his/her own advantage. Thus, for example, the development of a three-dimensional subgrid model in -coordinates should be rather straightforward.