A review on advanced numerical methods for free-surface hydrodynamics

Environmental hydrodynamics is typically characterized by free-surface hydrostatic flows. Within such a framework a set of simplified models are derived from the Reynolds averaged Navier–Stokes equations under various simplifying assumptions. Numerically, a simple method for the 1D open channel flow is derived first. This method is then extended to the 2Dxz\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2D_{xz}$$\end{document} laterally averaged model and to the 2Dxy\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2D_{xy}$$\end{document} vertically averaged model. Next, a semi-implicit method for the 3D model is described and discussed.


Introduction
The governing three-dimensional equations describing free-surface flows of an incompressible fluid are the well known Navier-Stokes equations which express the conservation of momentum and mass. Such equations in Cartesian coordinates have the following form (see, e.g., [1][2][3]) where u(x, y, z, t), v(x, y, z, t) and w(x, y, z, t) are the velocity components in the horizontal x, y and vertical z-directions; t is the time; p is the normalized pressure, that is the pressure divided by the constant density; g is the gravitational acceleration and ν is an eddy viscosity coefficient. Hereafter it is assumed that fluid flow is always confined between a prescribed solid bottom and a free-surface (see Fig. 1).
Assuming that the free-surface can be expressed as a single valued function z = η(x, y, t), the free-surface equation, also known as the kinematic condition of the free-surface, is given by where u s = u(x, y, η, t), v s = v(x, y, η, t) and w s = w(x, y, η, t) are the velocity components at the free-surface. Under the assumption that the bottom profile can be expressed as a single valued function z = −h(x, y), the total water depth is H (x, y, η) = h(x, y) + η(x, y, t) and a kinematic condition at the bottom boundary is where u b = u(x, y, −h, t), v b = v(x, y, −h, t) and w b = w(x, y, −h, t) are the velocity components at the bottom. Condition (3) states that the velocity component perpendicular to the solid boundaries must vanish. The tangential stress boundary conditions at the free-surface are specified by the prescribed wind stresses as follows where u a and v a are the horizontal wind velocity components and γ T is a nonnegative wind stress coefficient. Typically, γ T is taken to be γ T = C D ρ a ρ u 2 a + v 2 a where ρ a is the density of the air, ρ is the water density, ρ a ρ = 1.25 × 10 −3 and C D = 1.4 × 10 −3 is a drag coefficient.
The tangential boundary conditions at the sediment-water interface are given by specifying the bottom stress as follows ν(u z + u x h x + u y h y ) z=−h = γ B u b (6) where γ B is a nonnegative bottom friction coefficient. Typically, γ B is taken to be and C z is the Chezy friction coefficient which can be a constant (e.g., C z = 50), or related to a Manning's constant n by C z = H 1 6 n . Integration of the mass conservation equation over the depth yields thus, by using the kinematic boundary conditions (2)-(3) the following conservative form of the free-surface equation is obtained Equation (8) will replace (2) in the model formulations to be presented in the next Section (see, e.g., [4][5][6][7]).

Simplified mathematical models
In most environmental flows the vertical acceleration, as well as the vertical viscosity forces, are small when compared to the gravity acceleration and to the pressure gradient. Consequently, by neglecting the advective and the viscous terms in the third momentum equation, a simplified equation for pressure results This equation readily yields the following expression for the hydrostatic pressure where p a (x, y, t) is the atmospheric pressure at the free-surface which, without loss of generality, will be assumed to be a constant. Substitution of (9) into the Navier-Stokes equations yields the following simplified three-dimensional hydrostatic model (see, e.g., [4][5][6][7]) With properly specified initial and boundary conditions, Eq. (10) describe threedimensional, hydrostatic flows. Semi-implicit methods for the numerical treatment of this hydrostatic 3D model will be described in Sect. 6.
In lieu of solutions of the complete three-dimensional equations, the flow and circulation in a class of well-mixed estuaries, lakes and coastal embayments can be satisfactorily represented by the solution of a set of vertically averaged shallow water equations. The system of vertically averaged momentum equations can be derived by integrating vertically the momentum equations in (10) from the bottom z = −h to the free-surface z = η. Specifically, since H (x, y, η) = h(x, y) + η(x, y, t) is the total water depth, the vertically averaged velocities U and V are defined to be U = Similarly, by using the tangential stress boundary conditions (4) and (6), the vertical integration of the viscous terms yields  (10). Thus, the two-dimensional, vertically averaged shallow water equations are obtained from the free-surface equation (8) and from the vertical integration of the momentum equations in (10) which, after standard approximations of local velocities with their vertically average, yields the following 2D xy model (see, e.g., [8][9][10] where the overall friction coefficient γ = γ B + γ T must be nonnegative and can be reformulated to compensate for the specified approximations. Equation (11) constitute a system of three partial differential equations with three unknown functions U (x, y, t), V (x, y, t) and η(x, y, t). A semi-implicit numerical method for solving these equations will be discussed in Sect. 5.
Another alternative to the fully three-dimensional equations, is a simplified 2D xz model for narrow and deep estuaries which can be derived from the 3D equations under the assumption that the circulation of interest takes place in the vertical x-z plane. To this purpose let y = (x, z) and y = −r (x, z) be single value functions representing the left and the right walls, respectively, so that B(x, z) = (x, z) + r (x, z) ≥ 0 represents the width of the estuary (see Fig. 2).
The laterally averaged u-momentum equation can be derived by integrating umomentum equation in (10) from the right y = −r (x, z) to the left wall y = (x, z). To this purpose the laterally averaged velocities U and W , and the laterally averaged free-surfaceη are defined to be U = 1 B −r udy, W = 1 B −r wdy andη = 1 B −r ηdy, respectively. Thus, by using appropriate boundary conditions, and after standard approximations of the local variables with their laterally averaged quantities, the lateral integration of the u-momentum equation in (10) yields B −r νdy is the laterally averaged viscosity coefficient and the wall friction coefficient γ must be nonnegative and can be reformulated to partially compensate for the specified approximations.
Similarly, the lateral integral of the incompressibility condition in (10) yields Finally, upon integration of the free-surface equation in (10), one gets where A(x,η) = −r Hdy is the cross section area. In summary, then, the two-dimensional laterally averaged model is given by eqs. (12), (13) and (14), that is, With properly specified initial and boundary conditions, Eq. (15) form a 2D xz model that can be used to simulate water circulation in estuarine environment. A semi-implicit numerical method for solving these equations will be described in Sect.

4.
Finally, the one-dimensional shallow water equations for unsteady flows in open channel can be derived from either the two-dimensional models, or directly from the three-dimensional equations (10) by performing an average over the cross section. These equations, also known as the de Saint-Venant equations [10] are given by where, U (x, t) is the fluid velocity averaged over the cross section;η(x, t) is the laterally averaged free-surface elevation; A(x,η(x, t)) is the cross section area; and the overall friction coefficient γ = γ B +γ T must be nonnegative and can be reformulated to compensate for the specified approximations. Typically, γ B is taken to be γ B = g n 2 |U | R 1/3 , where R is the hydraulic radius defined as the ratio of the cross section area and the wetted perimeter.
Equations (16) constitute a system of two partial differential equations with two unknown functions U (x, t) andη(x, t). A semi-implicit numerical method for solving these equations will be discussed in the next Section. Figure 3 shows a schematic representation of the four differential models described above.
In solving numerically the differential equations (16) with simple explicit methods the occurrence of instability is a relatively common and disturbing event. Of course, a fully implicit discretization of these equations may lead to methods which are unconditionally stable. Implicit methods, however, may turn out to be unnecessarily complex. As a compromise, this Section will focus on the development of simple semi-implicit numerical methods which easily extend to multidimensional problems.
For the time being, the channel under consideration is assumed to have a rectangular cross section of constant width B = 1 so that the cross sectional area A simplifies to A(x, η) = H (x, η). Moreover, by using the lower case symbol u to denote the cross sectional average velocity, equations (16) in non conservative form are A severe limitation of standard explicit numerical methods for equations (17) is due to the stability restriction imposed by the Courant-Friedrichs-Lewy (CFL) condition [8][9][10][11]. Our first objective is to identify those terms in equations (17) whose implicit discretization enables to circumvent the CFL stability restriction of explicit methods.
Equations (17) in matrix notation can be written as In open channel flows the horizontal viscosity term 1 H (ν Hu x ) x is usually neglected. When ν = 0 system (18) is strictly hyperbolic and the corresponding characteristic Note that when |u| < √ g H the flow is subcritical and the characteristic speeds λ 1,2 have opposite directions (see Fig. 4). More importantly, the dominant term √ g H arises from the off diagonal terms g and H in matrix A. These are the coefficient of η x in the momentum equation and the coefficient of u x in the continuity equation. Therefore these derivatives must be discretized implicitly in order for the stability of the method to be independent of the celerity √ g H. The viscosity term and the non-linear advective term in the momentum equation will be discretized explicitly. However, for stability, the friction term in the momentum equation will be discretized implicitly, but the friction coefficient γ H will be evaluated explicitly so that the resulting approximation will be linear. The continuity equation will be considered in its original conservative form so that the differential equations being solved are Next, a staggered grid of N x equally spaced grid points is introduced. Let x and t be the space and the time steps, respectively. The discrete variables u and η are defined at alternate locations as shown in Fig. 5. Additionally, assuming that the bottom profile is prescribed everywhere, the discrete total water depths at the n-th time level are enforced to be nonnegative and are taken to be where η n i+ 1 2 is determined from the nearest grid values by taking, e.g., the average, the upwind or the maximum between η n i and η n i+1 . Then, a general semi-implicit discretization of equations (19) takes the form u n+1 In Eq. (22) F is any stable non-linear difference operator that includes a spatial discretization of the advective and viscous terms. A particular form for F can be chosen in a variety of ways thus, for example, by using a simple explicit upwind formula, or a more accurate and stable Eulerian-Lagrangian approach (for details see [4][5][6][7][8][9][10][11] and the references therein).
For any structure given to F, Eqs.  (20). This system has to be solved at each time step in order to determine the new values of the field variables from given initial data.
From a computational point of view, the system formed by Eqs.  (24) where, by setting G n Equation (24) can be assembled into a sparse, piecewise linear system of N x equations in which η n+1 i , i = 1, 2, . . . , N x are the only unknowns. This system takes the following form and T is the sparse, symmetric and tridiagonal matrix that arises from the linear terms in Eq. (24). At every time step, system (25) can be solved iteratively by a Newton type method in order to determine the new free-surface elevations η n+1 i . To this purpose note that the Jacobian matrix associated to system (25) is D(ζ ) + T where D(ζ ) is a diagonal matrix whose diagonal entries are Thus, a converging Newton type method for solving system (25) can be formulated as follows (see [12] for details): guess ζ 0 i > −h i for all i = 1, 2, . . . , N x . Then, until convergence, compute where the superscript denotes the iteration step (not to be confused with the time level The governing equations of a 2D xz model are the two-dimensional laterally averaged equations (15) in which, for simplicity, a constant width B = 1 is assumed. With this assumption equations (15) can also be written as follows The simplified free-surface boundary condition is given by prescribing the wind stress as and the simplified boundary condition at the sediment-water interface is specified by prescribing the bottom stress as where γ , γ B and γ T are nonnegative friction coefficients. Based on the characteristic analysis carried on the one-dimensional model, for stability, the gradient of surface elevation in the momentum equation and the velocity in the free-surface equation will be discretized implicitly. Moreover, in order to allow the use of arbitrarily small vertical resolution without compromising the stability of the resulting method, the vertical viscosity will also be discretized implicitly. An explicit discretization will be used for advection and for the horizontal viscosity.
The spatial mesh consists of rectangular cells of length x and height z k = z k+ 1 2 − z k− 1 2 . Each cell is numbered at its center with indices i and k. The discrete u velocity is then defined at half integer i and integer k; w is defined at integer i and half integer k (see Fig. 6). Finally, η is defined at integer i.
As in the one-dimensional model, nonnegativity of the total water depth at discrete level is enforced by specifying  In Eq. (32) F is an explicit, nonlinear finite difference operator which includes the contributions arising from the spatial discretization of the advective and the horizontal viscosity terms. A particular form for F can be chosen in a variety of ways thus, for example, by using a simple explicit upwind formula, or a more accurate and stable Eulerian-Lagrangian approach.
The values of u above free-surface and below bottom in (32) are eliminated by means of the boundary conditions (28)-(29) which are written in difference form as ν n+1 and ν n+1

semi-implicit finite difference discretization for the free-surface equation in (27) in conservative form is taken to be
Assuming that the computational domain is covered by a set of N x × N z computational cells, equations (32) and (35) constitute a piecewise linear system of at most N x (N z + 1) equations. The nonlinearity resides in the definition of the total water depth H i (η n+1 i ). This system has to be solved at each time step in order to calculate the new field variables u n+1 i+ 1 2 ,k and η n+1 i throughout the flow domain.
The system of N x (N z + 1) equations formed by (32) and (35) can be conveniently decomposed into a set of N x independent linear tridiagonal systems of N z equations, and one piecewise linear system of N x equations. Specifically, Eqs. (32) and (35) are first written in matrix form as where, along each water column the vectors U, Z , G and the tridiagonal matrix A are defined as follows: Since A is positive definite, A −1 is also positive definite and therefore ( Z) A −1 1Z is a non-negative scalar. Hence Eq. (38) constitute a well posed piecewise linear system of N x equations for η n+1 i which has the same size and structure as the corresponding system that was derived in Sect. 3 for the one-dimensional model. Thus, its unique solution can be determined by a converging Newton type method as described in Sect. 3.
Once the new free-surface location has been determined, Eq. (36) constitute a set of N x linear, tridiagonal systems with N z equations each. All these systems are independent from each other, symmetric and positive definite. Thus, they can be conveniently solved by a direct method.
Finally, by discretizing the continuity equation in (27) the new vertical component of the velocity in each water column can be computed recursively from bottom to top as Thus every time step can be summarized as a simple algorithm which is as follows

Numerical methods for the 2D xy model
When ν = 0 the two-dimensional shallow water equations (11) form a quasilinear hyperbolic system of partial differential equations with three unknown functions, namely u, v and η, and three independent variables, x, y and t. When ν > 0 this system is no longer hyperbolic but the stability of any explicit numerical scheme will still depend on its hyperbolic part, and hence on the celerity √ g H. In order to determine an efficient semi-implicit method whose stability is independent on the celerity, a characteristic analysis of the hyperbolic part of (11) is preliminarily investigated (see [8,9]). Thus Eq. (11) in non-conservative form are written as or, in matrix notation, The characteristic equation of system (41) is given by The triples (r , s, q) satisfying Eq. (42) are directions normal to the characteristic cone at its vertex (see Fig. 7). Equation (42) decomposes as follows and The local characteristic cone with vertex in (x,ȳ,t) consists of the line through (x,ȳ,t) and parallel to the vector (u, v, 1) and the cone whose equation is in fact, on the cone surface, the gradient of the left-hand side of (45) satisfies (44). Note that, whereas the first part of the characteristic cone depends only on the fluid velocity u and v, the second part, which is defined by Eq. (44), depends also upon the celerity √ g H. Note also that the term g H in Eq. (42) arises from the off diagonal terms g and H in the matrices A and B. These are the coefficient of η x in the first momentum equation, the coefficient of η y in the second momentum equation, and the coefficient of u x and v y in the continuity equation of system (40). Consequently, these derivatives must be discretized implicitly in order to obtain a method whose stability is independent of the celerity [9].
Based on the above analysis, the semi-implicit method derived in Sect. 3 will be extended to the following 2D xy model: A staggered mesh of N x × N y rectangular cells of length x and width y is introduced. Each cell is numbered at its center with indices i and j. The discrete u velocity is then defined at half integer i and integer j; v is defined at integer i and half integer j; η is defined at integer i and integer j (see Fig. 8). Additionally, assuming that the bottom profile h(x, y) is prescribed everywhere, the discrete total water depths at the grid locations are enforced to be nonnegative and are taken to be where η n i+ 1 2 , j and η n i, j+ 1 2 are determined from the nearest grid values by taking, e.g., the average, the upwind or the maximum. Then, a stable semi-implicit scheme for equations (46), is derived by using an implicit discretization for the gradient of surface elevation in the momentum equations and for the velocity in the continuity equation. Moreover, for stability, the friction terms in the momentum equations are discretized implicitly, but the friction coefficients γ and γ T will be evaluated explicitly. The remaining terms are discretized explicitly [8,9]. The basic semi-implicit algorithm is u n+1 where H n = 0 is assumed.
In Eqs. (50)-(51) F is an explicit, finite difference operator corresponding to the explicit spatial discretization of the advective and the viscous terms. A particular form for F can be chosen in a variety of ways thus, for example, by using an explicit upwind formula, or a more accurate and stable Eulerian-Lagrangian approach.
For any structure given to F Eqs. (50)-(52) constitute a piecewise linear system of 3N x N y equations and 3N x N y unknowns u n+1  Equations (53) constitute a well posed piecewise linear system of N x N y equations with unknowns η n+1 i, j . This system can be written in a more compact matrix form as where ζ is a vector containing the unknowns η n+1 i, j ; the piecewise linear vector function H(ζ ) represents the corresponding total water depths H i, j (ζ i, j ); T is the five-diagonal, symmetric and at least positive semi-definite matrix that arises from the linear terms in Eq. (53); and b is a known vector with the right-hand sides b n i, j of Eq. (53). A converging Newton type method for solving system (54) is formulated as follows (see [12] for details): Guess ζ 0 i, j > −h i, j for all i and for all j. Then, until convergence, compute where the superscript denotes the iteration step (not to be confused with the time level n) and D is a diagonal matrix whose diagonal entries are

Numerical methods for the 3D model
The semi-implicit finite difference methods described in Sections 4 and 5 can be further generalized to yield a practical algorithm for solving the fully three-dimensional equations (10) [4][5][6][7]. These equations can be written as The spatial mesh consists of boxes with length x , width y and height z k . Each box is numbered at its center with indices i, j and k. The discrete u velocity is then defined at half integer i and integers j and k; v is defined at integers i, k and half integer j; w is defined at integers i, j and half integer k. Finally, η is defined at integers i, j.
Then, on wet vertical faces perpendicular to the x-axis the discrete u momentum equation in matrix form can be written as system of N x N y equations for η n+1 i, j which has the same size and structure as the corresponding system that was derived in Sect. 5 for the 2D xy model. Thus, its unique solution can be determined by a converging Newton type method as described in Sect. 5.
Once the new free-surface location has been computed, equations (57) and (58) constitute a set of simple 2N x N y linear, tridiagonal systems with N z equations each. All these systems are independent from each other, symmetric and positive definite. Thus, they can be conveniently solved by a direct method.
Finally, by discretizing the continuity equation in (56), the new vertical component of the velocity in each water column can be computed recursively from bottom to top as Thus every time step can be summarized as a simple algorithm which is as follows (a) set up and solve the piecewise linear system (60) to obtain simultaneously the new water levels η n+1 and the resulting total water depths H (η n+1 ) throughout the horizontal flow domain; (b) solve the linear tridiagonal systems (57) and (58) to determine the new horizontal velocities u n+1 and v n+1 ; (c) diagnostically obtain the new vertical velocity w n+1 from (61).
Interestingly enough, if the vertical spacing z is taken to be so large that both the bottom and the free-surface always fall within one vertical layer, this three-dimensional algorithm reduces to the 2D xy model described in Sect. 5. In this particular case, in fact, the k index is always 1 and therefore it can be omitted from the notation. Moreover, Z reduces everywhere to H and the discrete momentum equations (57)-(58) simplify to (50)-(51). Additionally, the discrete free-surface equation (59) simplifies to (52).
The above three-dimensional method also reduces to the two-dimensional 2D xz model derived in Sect. 4 when it is applied to simulate the flow into a deep and narrow estuary. In this case, if y denotes the estuary width, the j index is always 1 and therefore can be omitted from the notation. Moreover, by assuming zero water depth H n i, j± 1 2 = 0 everywhere along the left and the right walls, respectively, one has V n+1 i, j± 1 2 = Z n i, j± 1 2 = 0, and Eqs. (57) and (59) reduce to (36) and (37) already derived in Section 4 for the 2D xz model.
In turn, the one-dimensional model derived in Sect. 3 can be regarded either as a particular case of the 2D xz model where N z = 1 is assumed, or as a particular case of the 2D xy model where N y = 1 is assumed (see Fig. 9). The fact that one-and two-dimensional methods can be derived from the present 3D model as particular cases is an interesting feature of the present formulation. This property leads to a single, and yet very efficient multidimensional computer model, that applies to arbitrarily complex flow regions.
Funding Open access funding provided by Università degli Studi di Trento within the CRUI-CARE Agreement. Open access funding provided by Universitá degli Studi di Trento within the CRUI-CARE Agreement.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.