Enriched Galerkin method for the shallow-water equations

This work presents an enriched Galerkin (EG) discretization for the two-dimensional shallow-water equations. The EG finite element spaces are obtained by extending the approximation spaces of the classical finite elements by discontinuous functions supported on elements. The simplest EG space is constructed by enriching the piecewise linear continuous Galerkin space with discontinuous, element-wise constant functions. Similar to discontinuous Galerkin (DG) discretizations, the EG scheme is locally conservative, while, in multiple space dimensions, the EG space is significantly smaller than that of the DG method. This implies a lower number of degrees of freedom compared to the DG method. The EG discretization presented for the shallowwater equations is well-balanced, in the sense that it preserves lake-at-rest configurations. We evaluate the method’s robustness and accuracy using various analytical and realistic benchmarks and compare the results to those obtained using the DG method. Finally, we briefly discuss implementation aspects of the EG method within our MATLAB/GNU Octave framework FESTUNG. ? Corresponding author. Moritz Hauck, Florian Frank Department of Mathematics, Friedrich-Alexander University Erlangen-Nürnberg, Cauerstraße 11, 91058 Erlangen, Germany Vadym Aizinger Chair of Scientific Computing, University of Bayreuth, Universitätsstraße 30, 95447 Bayreuth, Germany E-mail: vadym.aizinger@uni-bayreuth.de Hennes Hajduk Institute of Applied Mathematics (LS III), TU Dortmund University, Vogelpothsweg 87, 44227 Dortmund, Germany Andreas Rupp Interdisciplinary Center for Scientific Computing (IWR), Heidelberg University, Im Neuenheimer Feld 205, 69120 Heidelberg, Germany 2 Moritz Hauck et al.

of key requirements for SWE models).
The aforementioned issues led to a large number of studies dedicated to the development, analysis, and practical evaluation of various numerical techniques for solving the SWE. Historically, the first models used finite differences on structured grids, but, with the emergence of unstructured-mesh models 15 (e.g. TELEMAC [2] or ADCIRC [3]), finite elements and finite volumes became the de-facto standard. A big advantage of the finite element approach is its potential to naturally accommodate higher-order approximations on unstructured meshes; in this vein, various discretizations based on the continuous (CG) and discontinuous Galerkin (DG) approximation spaces (or mixtures of 20 both) have been described and compared in the literature [4,5]. The results of these comparisons can be summarized (in a somewhat oversimplified fashion) as follows: -Using continuous (e.g. piecewise linear) finite elements for elevation and velocity is computationally very efficient (at least the low-order schemes) 25 but tends to have stability issues usually represented by spurious elevation or velocity modes that arise from the LBB condition or from too large elevation spaces; -Using discontinuous Galerkin (e.g. piecewise linear) spaces for both unknowns is robust and needs no additional stabilization but is computation- 30 ally more expensive (up to a factor of four to five in serial execution [6]). In addition, discontinuous Galerkin schemes are locally conservative, better suited for adaptive and non-conforming meshes, and can partially offset their higher computational costs by more efficient parallel scaling [7]; some combinations of continuous and discontinuous spaces such as the 35 lowest order Raviart-Thomas spaces tend to perform better [4], but are rather complicated to extend to higher orders on triangular meshes.
The idea of the enriched Galerkin (EG) method is to enhance the approximation space of the continuous finite elements by adding element-local discontinuous functions and, by using the solution procedure nearly identical to that of the DG method (i.e. integrals of edge fluxes, Riemann solvers, etc.), produce a robust, locally conservative discretization with significantly fewer unknowns than a DG discretization of the same order.
In its original form, the EG method (also considered in this paper) adds a piecewise constant DG component to a continuous piecewise linear or multi- 45 linear CG discretization. This method was introduced for the linear advectiondiffusion-reaction equation in [8] and was proved there to be stable and converge at the same rate as the discontinuous Galerkin method for piecewise linear discretization. The piecewise constant enrichment of the linear continuous Galerkin approximation is intrinsically stable and provides the local conserva-50 tion property as was shown in [9].
Since then, the EG method has been generalized to the classical continuous Galerkin function space of the polynomial order k augmented by discontinuous, element-wise constant functions [9,10] and even to arbitrary enrichment with polynomials of degree m with −1 ≤ m ≤ k (where enrichment of degree m = 55 −1 represents no enrichment at all) as discussed in [11]. Arbitrary enrichment allows to consider EG as generalization of both CG and DG since EG uses the same bilinear and linear forms as DG.
Thus, EG inherits many advantages of DG. However, in multiple dimensions, EG has fewer degrees of freedom (DOF) than the DG method if m < k.

60
The standard EG method (with m = 0) has been developed to solve general elliptic and parabolic problems with dynamic mesh adaptivity [12,13,14,15] and was extended to address multiphase fluid flow problems [15]. Recently, the EG method has been applied to solve the nonlinear poroelastic problem [16,17], and its performance has been compared to other two-and three-field formulation methods [18].
In the context of convection-dominated problems, additional stabilization techniques are often required. Classical stabilizations include the streamline upwind Petrov-Galerkin (SUPG) method, weighted essentially non oscillatory (WENO) schemes, as well as provably bound-preserving alternatives (e.g. [19, 70 20, 21, 22, 23]). Popular methods designed particularly for DG discretizations include the edge-based Barth-Jesperson limiter [24], and its vertex-based counterparts [25,26]. Limiters for an EG discretization of the linear advection equation have recently been proposed in [27]. The approach therein deviates from classical DG slope limiters but rather fits in the framework of algebraic flux 75 correction [28], which only recently has been extended to the DG setting [29,30]. Numerical solutions based on the methods in [27] can be proven to satisfy discrete maximum principles under CFL-like time step restrictions, which makes the approach superior to geometrical slope limiting.
The main focus of the present work is to formulate and evaluate an EG 80 scheme for the SWE and to compare the quality of the EG and DG discretizations. The method is implemented in the FESTUNG framework [31,32,33,34] by modifying our DG implementation for the SWE presented in [35,36]. The same scheme was initially introduced in our UTBEST solver [37,38] and later extended to three dimensions in UTBEST3D [39,40].
The paper is structured as follows. The mathematical model and its discretization using the EG method are presented in Sec. 2, a brief description of the implementation using our MATLAB / GNU Octave framework FESTUNG is the subject of Sec. 3. In Sec 4, we demonstrate the accuracy and robustness of our EG scheme using an analytical convergence test, a supercritical 90 flow example with discontinuous solution, and a realistic tidal flow scenario for Bahamas islands. A short Conclusions section wraps up this work.

Governing equations
The SWE in conservative form are given by They are considered on a two-dimensional, polygonally-bounded domain Ω and finite time interval (t 0 , t end ). By ξ, we denote the free surface elevation of the water body with respect to a certain zero level (e.g., the mean sea level).
The quantity H = ξ − z b represents the total fluid depth with z b denoting the bathymetry. q := (U, V ) T is the depth integrated horizontal velocity field, f c the Coriolis coefficient, g the gravitational acceleration, and τ bf the bottom friction coefficient. Wind stress, the atmospheric pressure gradient, and tidal potential are combined in the body force term F := (F x , F y ) T . Defining c := (ξ, U, V ) T , system (1) can be rewritten in the following compact form: with and In this work, we use several types of boundary conditions for the SWE (1). Following the standard approach, interior values are used as boundary values at parts of the boundary, on which the corresponding unknowns are not prescribed. By· , we denote prescribed boundary values of the respective unknowns. The following types of boundary conditions are used in this work: Dirichlet boundary: Here, all unknowns are specified: Land boundary: Denoting by n the exterior unit normal to ∂Ω, we set the normal flux to zero: q · n = 0.
Open sea boundary: Only the free surface elevation is prescribed for this boundary type: ξ =ξ.
Radiation boundary: No unknowns are specified, which corresponds to a 100 free outflow condition.
Finally, we need to prescribe initial conditions for the elevation and depth integrated velocity

Enriched Galerkin discretization
Let {T ∆ } ∆>0 be a simplicial, shape-regular, quasi uniform, geometrically conformal family of triangulations of Ω ⊂ R 2 with #T denoting the total number of elements of T ∆ . We obtain the local variational formulation of system (2) by multiplying with smooth test functions φ ∈ C ∞ (Ω) 3 and integrating by parts on each element T ∈ T ∆ yielding where we write (·, ·) T and ·, · ∂T for the L 2 -scalar products on elements and their boundaries, and denote by n = (n x , n y ) T an exterior unit normal to ∂T . Defining the broken polynomial spaces of order m ∈ N 0 as : v |T is a polynomial of degree at most m,∀T ∈ T ∆ and P −1 (T ∆ ) := {0}, we can modify them to obtain the EG test and trial spaces Here, '+' denotes the sum of subspaces which is not a direct sum if m = −1. Examples of spaces are given in Fig. 1. From (10) it follows immediately that This dimension formula will become handy in Sec. 3 for determining EG shape 105 functions. For more details about these spaces and their properties, the reader is referred to [11,Sect. 3].
To obtain the semi-discrete EG formulation of (9), c and φ are replaced by their discrete counterparts c ∆ , φ ∆ ∈ P k,m (T ∆ ) 3 . Since the values of a discontinuous function are not unique on element edges, we replace the boundary 110 term A(c) · n by a numerical fluxÂ(c ∆ , c + ∆ , n) that depends on the discontinuous values of the solution on element T (without superscript) and its edge neighbor (superscript + ). On domain boundaries, the specified boundary values of free surface elevation and velocity are utilized in place of c + ∆ for the flux computation.

115
In conclusion, summing up over all elements T ∈ T ∆ , yields In this work, we use the Lax-Friedrichs flux combined with the Roe-Pike averaging [38, 41,42] defined aŝ where λ = λ(c ∆ , c + ∆ , n) is given by where we exploit the fact that the bathymetry z b is assumed to be continuous.
Discrete initial conditions are obtained using suitable projections of ξ 0 , and q 0 into the discrete ansatz space (10) (cf. [11, Sect. 5] for a comparison of reasonable choices). For the temporal discretization, we use strong stability 120 preserving (SSP) explicit Runge-Kutta methods [43].

Implementation aspects
Our implementation of the EG method is based on the DG code for the SWE realized within the FESTUNG framework [35,44]. To obtain the required EG operators, we use a simple strategy of modifying the DG code. This trick is 125 by no means restricted to the SWE system; it can be just as easily applied to any linear or nonlinear PDE. The main idea of our approach is to exploit the fact that EG (or, for that matter also CG) approximation spaces are embedded in the DG spaces of the same order. Therefore, all EG or CG basis functions can be represented as linear combinations of DG basis functions; 130 this representation can be written in the form of a (generally rectangular) system of linear equations, which can then be used to convert an available DG operator into the corresponding one for EG or CG spaces. The downside of this approach is that a full DG discretization must be assembled firstthus limiting the efficiency gain due to smaller approximation spaces of the 135 EG method. In the solution procedure, a linear system of equations with the mass matrix needs to be solved. We observed some speedup in the solution procedure for the EG method when choosing a diagonal mass matrix (e.g., by lumping) or semi-implicit time discretizations. This speedup is due to a smaller size of the linear equation system compared to DG.

140
In this work, we utilize finite element spaces of polynomial degree k ≤ 2. Let us first consider the DG space P k (T ∆ ) and denote by N := dim P k (T ∆ ) the number DG unknowns for a scalar quantity. We have dim P 1 (T ∆ ) = 3 #T and dim P 2 (T ∆ ) = 6 #T . In the following, we construct bases for various EG spaces from CG and DG basis functions. Here, one has to be somewhat careful since, in general, simply combining a CG basis with element-wise DG basis functions produces 150 a linearly-dependent set. For the admissible combinations of k and m with −1 ≤ m < k ≤ 2 considered in our work, the EG bases can be constructed as follows: k ∈ {1, 2}, m = −1: The space P k,−1 (T ∆ ) coincides with the CG ansatz space of order at most k. Therefore, we can simply choose the CG basis k = 2, m = 1: We use the standard linear DG basis and extend it by the nodal quadratic shape functions equal to 1 at one edge midpoint, but omit the ones corresponding to cell vertices. The functions are linearly inde-165 pendent and are contained in the space P 2,1 (T ∆ ). As the number of basis functions equals dim P 2,1 (T ∆ ) = 3 #T + #E, it must be a basis of P 2,1 (T ∆ ).
Here, #T and #E denote the number of triangles and edges of T ∆ , respectively. The dimension formula (11) ensures that this, in fact, constitutes a basis of P 2,1 (T ∆ ).

Assembling nonlinear EG operators from DG operators
Next, we show how to obtain an EG from the corresponding DG operator.
To simplify the presentation, we formulate our approach here for the scalar case noting that the generalization to vector fields is straightforward. DG discretizations of nonlinear PDEs such as the SWE feature nonlinear operators of the form where b : P k (T ∆ ) × P k (T ∆ ) → R is linear in the second argument. Our goal is to form the EG operator by making use of (14). Since P k,m (T ∆ ) ⊂ P k (T ∆ ), it is possible to express the EG basis functions as linear combinations of the DG basis Next, we insert (16) into (15) C jl y j ϕ l , ϕ k for i ∈ {1, . . . , M }. By defining C := (C kl ) kl , k ∈ {1, . . . , M }, l ∈ {1, . . . , N }, we can thus write Employing (17), we assemble the terms corresponding to (15) by modifying an existing DG discretization of the operator (14) and computing the matrix C. Due to the above choice of EG and DG basis functions, the matrix C can be determined from simple geometric considerations during preprocessing.
is a bilinear form, we can write its induced operator (mass matrix) as and can be obtained from (17), and, in operator form, is given by That is, for operators that are induced by bilinear forms, the assembly can be preprocessed.

Numerical results
In this section, we investigate the performance of the EG method using artifi-180 cial and realistic test problems for the SWE. The main goals of these numerical studies can be summarized as follows: -Verify the expected rate of convergence against a manufactured analytical solution; evaluate the solution quality for realistic benchmarks;

185
compare the EG and DG methods in terms of accuracy, stability, and robustness.
We denote the results obtained for different combinations of k and m by the corresponding finite element spaces. To simplify notation, we omit their dependency on T ∆ . Hence, we write P k,m instead of P k,m (T ∆ ) and use the 190 convention that P k,k is the DG space of polynomial degree k.
As temporal discretization, we use an SSP Runge-Kutta time stepping scheme described in [45] with s = k + 1 stages. The time step size depends on the specific test problem.

Analytical convergence test 195
For our first numerical experiment, we approximate a smooth solution of the SWE on the square domain Ω := (0, 1000) × (0, 1000). The coarsest mesh (corresponding to level one) consists of 16 triangular elements and is shown in Fig. 2 (left). To investigate the convergence behavior, we consider a total of five meshes obtained from the coarsest mesh by uniform refinement via edge 200 bisection.
Parameters τ bf and f c are set to zero. Substituting the following analytical solution ξ(x, y, t) : V (x, y, t) := C 1 + C 2 C 3 sin π(x + y + C 3 t) 600 . of ∆t is sufficiently small to make temporal discretization errors negligible compared to spatial approximation errors. Solution parameters are chosen as The projected initial surface elevation for P 1,0 on the coarsest mesh is shown in Fig. 2  convergence can be observed for k = 2 and m ∈ {1, 2}, but not for the CG approximation k = 2, m = −1. The results for k = 2 and m = 0 are slightly better than without the piecewise constant enrichment, but do not exhibit third order accuracy. In conclusion, the EG scheme converges almost exactly 220 as well as the DG method if m = k − 1, k ∈ {1, 2}, although the absolute errors tend to be somewhat larger than for their DG counterparts. This is to be expected because EG has fewer unknowns, and even the projected exact solution in general becomes less accurate if the number of unknowns is reduced.   Table 1 Analytical convergence test: L 2 (Ω) errors Err(·) at the final time and experimental orders of convergence EOC(·) for the surface elevation and depth integrated velocity components.
In Tab. 2, we list the total numbers of degrees of freedom for all configu-225 rations. Note that, the DOF for the EG method with element-wise constant enrichment (k ∈ {1, 2}, m = 0) has approximately half as many DOF as the DG method of order k. EG with k = 2, m = 1 has approximately three quarters of the DOF for the quadratic DG method, thus indicating a potential gain in computational cost for EG compared to DG. In conclusion, the EG method 230 performs similarly to DG for smooth solutions while having significantly fewer DOF.

Supercritical flow in a constricted channel
In order to demonstrate the stability and robustness of the EG method, we solve the supercritical flow problem proposed by Zienkiewicz and Ortiz in [46]. The computational domain is a channel whose lateral boundary walls are constricted on both sides with an angle of five degrees (cf. Fig. 4). This benchmark uses constant bathymetry z b ≡ −1 while parameters τ bf , and f c are once more set to zero. The following initial and boundary conditions are prescribed for this problem: Initially, we set the surface elevation and momentum to ξ 0 ≡ 0, q 0 = (1, 0) T , respectively. Land boundary conditions are imposed at the lateral (wall) boundaries, while at the inlet (left), free surface elevation as well as velocity are for all times specified to be identical to their corresponding initial values. Finally, at the outlet (right) radiation boundary conditions are used. Denoting by u and H the axial velocity and water depth at the inlet, respectively, the flow regime is set to be supercritical by choosing the inlet Froude number Fr which is equivalent to setting the gravitational constant g = 0.16 m s 2 in this artificial example.

235
The solution to this problem converges to a steady-state, for which an analytical solution is available (see [47]). Fig. 4 illustrates the computational domain along with the unstructured mesh used in all computations. We run all simulations to a steady-state using pseudo time stepping with a time step of ∆t = 1/10 and present the results for all schemes in Fig. 5.
The steady-state surface elevation shown in Fig. 5 (top left) is discontinuous and displays interactions of waves reflected from the channel constrictions. Fig. 5 (top right to bottom right), display the steady-state solutions for all considered EG and DG methods (0 ≤ m ≤ k ≤ 2). Since no limiter was used, all approximations exhibit spurious oscillations close to the discontinuities. We 245 expect the results to improve and the numerical solutions to become boundpreserving if a limiter, such as the one developed in [27], is utilized (at least) for the free surface elevation.
The type of enrichment appears to play a major role for the stability of EG schemes. Thus the results for enrichments with m = k − 1, k ∈ {1, 2} are 250 almost indistinguishable from the corresponding DG results. In particular, the characteristic wave features such as positions and magnitudes of discontinuities are in good agreement with the analytical solution. On the other hand, the EG solution for k = 2, m = 0 exhibits severe oscillations not only in the vicinity of the discontinuities but also in the remainder of the computational 255 domain. This phenomen raises the concern a piecewise constant enrichment of the quadratic CG space is not sufficient to obtain an intrinsically stable scheme. However, an enrichment by piecewise linear discontinuous functions seems to remedy this issue, which confirms the previous findings in 4. 1 and [11]. Our study of this benchmark suggests that optimal EG schemes (m = k − 1) 260 possess similar stability properties to their DG counterparts while offering potential advantages in computational efficiency.

Tidal flow at Bahamas Islands
The realistic benchmark considered in this work involves a tidal flow scenario around the Bahamas Islands. The domain geometry, bathymetry as well For the bottom friction coefficient we use the standard quadratic friction law τ bf = C f |q|/H 2 (see e.g. [48]) with coefficient C f = 0.009. The (constant) Coriolis parameter is set to 3.19 × 10 −5 s −1 . The following tidal forcing is prescribed at the open sea boundary:  with time t in hours. In real-life ocean simulations, the initial conditions are often unknown or very difficult to obtain, therefore a cold start initialization is performed: The flow domain is assumed to be at rest initially (ξ 0 ≡ 0, q 0 ≡ (0, 0) T ). Then, starting immediately at initial time t = 0, boundary 275 condition (18) is imposed as A tidal forcing driving the flow. The simulations are run for a total of twelve days using constant time step ∆t = 15 seconds for EG and DG discretizations corresponding to 0 ≤ m ≤ k ≤ 2. In Fig. 7-9, we compare the numerical solutions for all considered EG and DG methods at the four recording stations. 280 Fig. 7 demonstrates excellent agreement for the surface elevation for all EG and DG methods. The curves lie nearly on top of each other and no differences can be visually detected at this resolution. The differences in the surface elevation between the EG and DG approximations at the recording stations are on the order of 10 −4 meters.

285
In Fig. 8 and Fig. 9, we plot both velocity components at the recording stations. For the recording station two ( Fig. 8 and Fig. 9 top right), one can observe slight differences between the approximations of orders one and two. Such behavior is fully consistent with the station comparisons performed in [38]. The effect of approximation order has greater magnitude than the dif-290 ferences between the EG and DG discretizations. Small differences between EG and DG discretizations of the same order can be observed in the plots at some locations. The deviations at the recording stations are on the order of 10 −3 to 10 −4 m s −1 .
In addition to testing the accuracy and robustness of the EG method, we 295 also used the Bahamas example to verify the well-balancedness of the method. For this purpose, we ran our simulation for the lake-at-rest configuration with open sea boundary condition set to zero instead of using (18). Just as expected, no spurious circulation emerged for any of the EG or DG configurations.

300
In this work, we present the first enriched Galerkin discretization for the system of 2D shallow-water equations, evaluated its performance in analytical and realistic test problems, and compared the numerical results to those obtained using our discontinuous Galerkin solver. The results of our studies demonstrate that EG schems with enrichments using discontinuous spaces of order one less than the order of the continuous space display similar accuracy and robustness as the corresponding DG discretizations. Similarly to DG, EG discretizations guarantee local conservation of all conserved unknowns while the total number of degrees of freedom is substantially lower than that for the DG space of the same order. This makes the enriched Galerkin method a very attractive can-310 didate for solving the shallow-water equations and other nonlinear hyperbolic PDE systems. Several interesting avenues of future research concerning the development of the EG solver for the SWE present themselves at this time. Thus one could employ the quadrature-free methodology to improve the computational effi-315 ciency of the method-similarly to our recent work for the DG SWE solver in [49]. Also EG discretizations appears to hold promise (perhaps even more so than the DG method) for bathymetry reconstruction using modified shallowwater equations [50]. Finally, stabilization techniques, such as limiters, e.g., the one designed specifically for EG in [27], will be required for certain applications of the SWE such as dam break problems.