A continuous/discontinuous Galerkin solution of the shallow water equations with dynamic viscosity, highorder wetting and drying, and implicit time integration
Abstract
The highorder numerical solution of the nonlinear shallow water equations (and of hyperbolic systems in general) is susceptible to unphysical Gibbs oscillations that form in the proximity of strong gradients. The solution to this problem is still an active field of research as no general cure has been found yet. In this paper, we tackle this issue by presenting a dynamically adaptive viscosity based on a residualbased subgrid scale model that has the following properties: it removes the spurious oscillations in the proximity of strong wave fronts while preserving the overall accuracy and sharpness of the solution. This is possible because of the residualbased definition of the dynamic diffusion coefficient. For coarse grids, it prevents energy from building up at small wavenumbers. The model has no tunable parameter.
Our interest in the shallow water equations is tied to the simulation of coastal inundation, where a careful handling of the transition between dry and wet surfaces is particularly challenging for highorder Galerkin approximations. In this paper, we extend to a unified continuous/discontinuous Galerkin (CG/DG) framework a very simple, yet effective wetting and drying algorithm originally designed for DG [Xing, Zhang, Shu (2010) Positivitypreserving high order wellbalanced discontinuous Galerkin methods for the shallow water equations Adv. Water Res. 33:14761493]. We show its effectiveness for problems in one and two dimensions on domains of increasing characteristic lengths varying from centimeters to kilometers.
Finally, to overcome the timestep restriction incurred by the highorder Galerkin approximation, we advance the equations forward in time via a three stage, second order explicitfirststage, singly diagonally implicit RungeKutta (ESDIRK) time integration scheme. Via ESDIRK, we are able to preserve numerical stability for an advective CFL number 10 times larger than its explicit counterpart.
1 Introduction
The shallow water equations (SaintVenant SaintVenant1871 ()) are a common approximation to the –dimensional NavierStokes equations to model incompressible, free surface flows. Due to the ability of highorder Galerkin methods to keep dissipation and dispersion errors low (Ainsworth et al. ainsworthEtAl2006 ()) and their flexibility with arbitrary geometries and hpadaptivity, these methods are proving their mettle for solving the shallow water equations (SW) in the modeling of nonlinear waves in different geophysical flows (Ma ma1993 (), Iskandarani iskandaraniEtAl1995 (), Taylor et al. taylorTribbia1997 (), Giraldo giraldo2001 (), Giraldo et al. giraldoHesthavenWarburton2002 ()), Dawson and Aizinger dawsonAizinger2005 (), Kubatko et al. kubatkoEtAl2006 (), Nair et al. nairThomasLoft2007 (), Giraldo and Restelli giraldoRestelli2010 (), Xing et al. xingZhangShu2010 (), Kärnä et al. karnaEtAl2011 (), Hendricks et al. hendricksKoperaGiraldo2015 (), Hood karolineTHESIS2016 ()). Although it is typically assumed that highorder Galerkin methods are not strictly necessary, they do offer many advantages over their loworder counterparts. Examples include their ability to resolve fine scale structures and to do so with fewer degrees of freedom, as well as their strong scaling properties on massively parallel computers (Müller et al. mullerEtAl2016 (), Abdi et al. abdiEtAl2016 (), Gandhem et al. gandhamEtAl2015 ()). Highorder methods are often attributed with some disadvantages as well. For example, they are constrained by small timesteps. To overcome this restriction, we follow Kärnä et al. karnaEtAl2011 () and implement an implicit RungeKutta scheme based on Giraldo et al. giraldoEtAl2013 (). Furthermore, the numerical approximation of nonlinear hyperbolic systems via highorder methods is susceptible to unphysical Gibbs oscillations that form in the proximity of strong gradients such as those of sharp wave fronts (e.g. bores). Filters such as Vandeven’s and Boyd’s vandeven1991 (); boyd1996 () are the most common tools to handle this problem, as well as artificial diffusion of some sort. However, we noticed that filtering may not be sufficient as the flow strengthens and the wave sharpness intensifies. We have recently shown this issue with the solution of the nonlinear Burgers’ equation by high order spectral elements in (marrasEtAl2015b, , §5) and will show how a proper dynamic viscosity does indeed improve on filters for the shallow water system as well. To preserve numerical stability without compromising the overall quality of the solution, Pham Van et al. phanVanEtAl2014 () and Rakowsky et al. rakowskyEtAl2013 () utilized a LillySmagorinsky model lilly1962 (); smagorinsky1963 (). To account for subgrid scale effects, viscosity was also utilized in the DG model described by Gourgue et al. gourgueEtal2009 () to improve their inviscid simulations. Recently, Pasquetti et al. pasquettiGuermondPopoICOSAHOM2014 () stabilized the highorder spectral element solution of the onedimensional SaintVenant equations via the entropy viscosity method. Michoski et al. michoskiEtAl2016 () compared artificial viscosity, limiters, and filters for the (modal) DG solution of SW, concluding that a dynamically adaptive diffusion may be the most effective means of regularization at higher orders. The diffusion model that we present in this work stems from the subgrid scale (SGS) model that was proposed by Nazarov and Hoffman nazarovHoffman2013 () to stabilize the linear finite element solution of compressible flows with shock waves. This was later applied to stabilize highorder continuous and discontinuous Galerkin (CG/DG) in the context of stratified, low Mach number atmospheric flows by Marras et al. in marrasEtAl2015b (). This approach is based on the ideas of scales splitting in large eddy simulation, where the scales of physical importance are split into grid resolvable and unresolvable. The unresolved scales are then parameterized via the SGS model at hand (DynSGS). DynSGS is applicable to any numerical method and is not limited to Galerkin approximations; its implementation in existing numerical models is straightforward. Stabilization becomes more important as the grid resolution is decreased. The metric that we utilize to evaluate its importance is the spectrum of kinetic energy across wave numbers. To understand the impact of diffusion with respect to the energetic behavior of the numerical solution, we compare the spectra obtained for the stabilized continuous and discontinuous Galerkin solutions against their inviscid counterparts. As viscosity is accounted for, the solution improves at coarser resolutions for both CG and DG. In the case of DG, however, dissipation is already builtin by the flux computation across elements. For this reason, the amount of artificial diffusion that DG requires for stabilization is smaller than the one required by CG. This property of DG has been exploited to build implicit Large Eddy Simulation (ILES) models (e.g., Uranga et al. urangaEtAl2011 ()). ILES relies on the numerical dissipation of the discretization method to model the effects on the unresolved, subgrid scale eddies rather than using an explicit subgrid scale (SGS) model, under the requirement of highresolution rider2006 (); drikakisEtAl2007 (). The MPDATA method of Smolarkiewicz smolarkiewicz1983 (); smolarkiewicz1984 () is an example of such an ILES model used for modeling geophysical flows.
Finally, there is a known difficulty in including wetting and drying algorithms (typically designed for loworder methods) while preserving highorder accuracy. The application of wetting/drying with discontinuous Galerkin using loworder Lagrange polynomials can be found in, e.g., Bunya et al. bunyaEtAl2009 (), Vater et al. vaterEtAl2015 (), Kärnä et al. karnaEtAl2011 (), or Gourgue et al . citegourgueEtal2009, and using Bernstein polynomials up to order three in Beisiegel and Behrens beisiegelBehrens2015 (). The positivity preserving limiter of Xing et al. xingZhangShu2010 () (XingZhangShu limiter from now on) was designed for highorder discontinuous Galerkin to solve this problem in particular. Because it is mass conservative and preserves the global high order accuracy of the solution, we implemented it in our model and extended it to continuous Galerkin as well. This limiter is guaranteed to work if and only if the mean water height is positive. Its designers recommend using the total variation bounding (TVB) limiter of Shu shu1987 (); cockburnShu1989 () to make the solution positive definite. Because the TVB limiter requires a bound on the derivative of the solution, which we do not know a priori, we rely on the artificial dissipation introduced above instead.
Looking at the larger picture of things, the ever increasing interest in the highorder accuracy of inundation models stems from a long list of coastal disasters in the last 10 years alone. From the 2004 tsunami in the Indian Ocean that caused 280,000 deaths, to the 2011 Tohoku tsunami, followed by Superstorm Sandy in 2012, the devastating 2013 typhoons in the Philippines and India, and the 2014 hurricane Odile; the list seems to be getting longer rapidly as the global climate patterns are changing. To study the impact of similar events in the future, coastal planners around the world are more and more relying on numerical tools such as the one that we present in this paper.
The rest of the paper is organized as follows. The shallow water equations are defined in §2. Their numerical solution is described in §3. The inundation model is described in §5 and is followed by the derivation of the SGS model in §4. Numerical results are reported in §6. The conclusions are drawn in §7.
2 Equations
Let be a fixed domain of space dimension with boundary and Cartesian coordinates in 1D and in 2D and let identify time. Let always identify the direction of gravity. The absolute free surface level and bathymetry are identified by the symbols and, respectively, , so that . The viscous shallow water equations are then given as:
(1a)  
(1b) 
where is the acceleration of gravity, is the identity matrix, and is the dynamic dissipation coefficient to be defined shortly. In (1a), the coefficient defines whether viscosity is turned on () or off () in the continuity equation. We will return to this topic shortly.
3 Space and time discretization
Highorder spectral element (SEM or CG, for continuous Galerkin) and discontinuous Galerkin (DG) approximations on quadrilateral elements are used to discretize Eq. (1).
The numerical model used in this paper is in fact the twodimensional version of the NUMA model described in Giraldo et al. giraldoHesthavenWarburton2002 () and in Abdi and Giraldo abdiGiraldo2016 (). Furthermore, the model is derived from the model described in Kopera and Giraldo koperaGiraldo2013a (); koperaGiraldo2013b () which is here used for the shallow water equations. The solution is advanced in time using a fully implicit RungeKutta scheme (see §3.2).
3.1 Spectral element and discontinuous Galerkin approximations
We leave the details of the discretization to the work of Giraldo et al. giraldoHesthavenWarburton2002 () and Kopera and Giraldo koperaGiraldo2013a (); koperaGiraldo2013b (). Nonetheless, we introduce some notation that we are going to use later in the paper. To solve the shallow water equations by elementbased Galerkin methods on a domain , we proceed by defining the weak form of (1) that we first recast in compact notation as
(2) 
where is the transposed array of the solution variables and and are the flux and source terms.
In the case of spectral elements, the space discretization yields the semidiscrete matrix problem
(3) 
where, for the global mass and differentiation matrices, and , we have that . The global matrices are obtained from their local counterparts, and , by direct stiffness summation, which maps the local degrees of freedom of an element to the corresponding global degrees of freedom in , and adds the element values in the global system. By construction, is diagonal (assuming inexact integration), with an obvious advantage if explicit time integration is used.
In the discontinuous Galerkin approximation, the problem at hand is solved only locally, and unlike the case of spectral elements, the flux integral that stems from the integration byparts must be discretized as well. The elementwise counterpart of the matrix problem (3) is then written as:
(4) 
where we obtain from the element boundary matrix, , and the element mass matrix, . Out of various possible choices for the definition of the numerical flux in Eq. (4), we adopted the Rusanov flux. The Laplace operator of viscosity is approximated using the Symmetric Interior Penalty method (SIP; the reader is referred to Arnold arnold1982 () for details on its definition).
3.2 Time integration
Equation (3) is integrated in time by an implicit RungeKutta (RK) scheme that corresponds to the implicit part of the implicitexplicit scheme used in giraldoEtAl2013 () (see also Butcher_1999 ()). The method coefficients in standard () tableaux form are the following
(5) 
Scheme (5) is a three stage second order explicitfirststage singly diagonally implicit RK (ESDIRK) scheme. This scheme has desirable accuracy and stability properties: all stages are second order accurate, it is stiffly accurate and stable, and it is strongstabilitypreserving (SSP) gottlieb2001strong () with SSP coefficient of 2. These properties allow us to take large timesteps with high accuracy as well as alleviate the instability issues associated with sharp solution gradients gottlieb2001strong (). The twodimensional tests presented later in this paper showed to be the most demanding in terms of stability constraints. Method (5) allows us to gain up to one order of magnitude in terms of maximum admissible advective CFL when compared to an explicit method (explicit part of ARK3, Kennedy_2001 ()). In particular, the explicit fourstage RungeKutta solution of the solitary wave against one isolated obstacle described in §6.5 preserved stability for up to CFL=0.21 using both CG and DG approximations. Although we were not able to use arbitrarily large timesteps with the ESDIRK with the current implementation (we will address this issue in a future work), we resolved the same problem at CFL = 1.8. Schemes with a subset of these properties are employed by Kärnä et al. karnaEtAl2011 () and shown to be robust in this context. Method (5) used in this study satisfies all properties ().
Computationally, at each of the two implicit stages we have to solve a nonlinear equation , where are the stage values, . We do so by using Newton iterations with a stopping criterion based on the relative decrease in the residual; that is, stop at iteration if . At each Newton iteration we have to solve a linear system , where is the Jacobian matrix of . We approximate the Jacobian using directional finite differences and iterate with the generalized minimal residual (GMRES) method, which is effectively a Jacobianfree Newton–Krylov method knoll2004jacobian (). The GMRES stopping criterion is also based on the relative residual. The first stage is explicit and equal to the last stage of the previous step, effectively making it a twostage method, which saves some computational time.
4 The DynSGS model for the shallow water equations
There are different ways to derive the viscous model described by Eq. (1) from the inviscid SaintVenant equations. Analogous to our previous work on the large eddy simulation of stratified atmospheric flows marrasEtAl2015b (), the current model builds on the separation between grid resolved (indicated as for any quantity ) and unresolved (subgrid) scales. The unresolved scales are modeled via the eddy viscosity scheme described in this paper. Given an element of order and with side lengths of comparable orders of magnitude, we define the following characteristic length (and hence filter width sagautBook ()):
The value of sets the size of the smallest resolvable scales in the same way as cutoff filters do in large eddy simulation models.
The application of this filter to the continuity equation (1a) results in the presence of an additional subgrid term on the righthand side. It is often debated whether diffusion should be applied to the continuity equation pasquettiGuermondPopoICOSAHOM2014 (); gerbeauEtAl2001 (); Guermond_Popov_2014 (); clearly, should the discrete viscous operator not be conservative, mass dissipation should not be used. However, by relying on spectral elements with integration by parts of the secondorder diffusion operator, the discrete viscous operator is conservative, as shown in gubaEtAl2014 (). To get a sense of how necessary a stabilized continuity equation may be, we will show a few results for both conditions in §6.5.
Scale separation in the momentum equation yields a new equation that includes the gradient of the quantity
(6) 
where the indicates subgrid velocity. The coefficient is defined elementwise and is given as:
(7) 
where
(8) 
and
(9) 
In (8, 9), indicates the spatially averaged value of the quantity at hand over the global domain , the norms at the denominator are used to preserve the physical dimension of the resulting equation, and are the residuals of the inviscid governing equations. At each timestep, the residuals are known. The presence of makes the artificial diffusion mathematically consistent^{1}^{1}1In the context of stabilization, mathematical consistency is the property of artificial diffusion such that it vanishes when the residual is zero.. The quantity in (8, 9) is the maximum wave speed.
Remark 1: physical dimensions
We would like to emphasize the necessity for the physically correct dimensions of diffusion. This is an important issue that is often underestimated and not accounted for in the design of artificial diffusion methods for stabilization purposes.
5 Wetting and drying interface
The dry regions are handled by adding an infinitesimal layer of water on the dry surfaces. We utilized a threshold value of 1e3 m^{2}^{2}2One millimeter of water is physically negligible. for all the tests that we ran. Furthermore, the limiter of Xing et al. xingZhangShu2010 () is applied on the velocity and water depth. Even at high order, this approach is suitable for wetting and drying problems solved via both continuous and discontinuous Galerkin methods.
Albeit simple, the XingZhangShu limiter works well. However, it is guaranteed to work if and only if the mean water height is positive. Xing et al. xingZhangShu2010 () recommend using the total variation bounding (TVB) limiter of Shu shu1987 (); cockburnShu1989 () to make the solution positive definite. We rely on artificial dissipation to achieve this as the TVB limiter requires a bound on the derivative of the solution, which we do not know a priori.
6 Numerical tests
We verify the correctness of our numerical model through five standard benchmarks in both one and twodimensions. In 1D, given the Nwave by Carrier et al. carrierEtAl2003 () that mimics a wave generated by an offshore submarine landslide, we compute the solution of the tsunami approaching a sloping beach. See §6.1. The numerical solution is evaluated against a set of tabulated surface height and momentum data provided by tsunamiWorkshopData (). The next 1D problem consists of the oscillation of a flat lake in a parabolic bowl gallardoEtAl2007 (); xingZhangShu2010 (). The existence of an analytic solution to this problem allows us to measure the accuracy of the numerical solutions. See §6.2. In addition, various 1D wetting and drying test cases were analyzed in the thesis of Hood karolineTHESIS2016 () to assess our model.
In 2D, we analyzed the results for the following tests. The vertical oscillation of a column of water in an axisymmetric paraboloid with analytic solution thacker1981 (). We describe this problem in §6.3 to test the inundation against the twodimensional effects of a varying topography. The second twodimensional test involves the flooding in a closed channel as presented by Gallardo et al. gallardoEtAl2007 () and by Xing and Zhang xingZhang2013 (). No exact solution is presented for this test, but the results can be compared against previous studies. This is shown in §6.4. We conclude the set of 2D problems with the simulation of the flow in a twodimensional channel with an isolated obstacle synolakis1987 (). We use this test to gain insight into the practical – other than theoretical as stated in gerbeauEtAl2001 () – need to include viscosity by analyzing the energetics of the solutions. The test is run during a sufficiently long time to assess the robustness of the inundation and viscosity schemes to handle fast motion with interacting waves that impinge against a steep obstacle. See §6.5. The description and analysis of more tests can be found in the collection recently compiled by Delestre et al. delestreEtAl2014 ().
6.1 1D tsunami runup over a sloping beach
The runup of a long wave on a uniformly sloping beach was originally proposed at the third international workshop on longwave runup models tsunamiWorkshopData (). The onedimensional computational domain is defined as m. The dry initial beach is 500 m long. The initial waveform was defined by Carrier et al. in carrierEtAl2003 () for an m domain as:
(10) 
with constants . To scale the wave to the m long domain used for the current test, the scaling factor is introduced and Eq. (10) is reexpressed with respect to and , with larger amplitudes . The initial wave is plotted in Fig. 1.
The solutions at times s are plotted in Fig. 2 and are compared against the tabulated data available in tsunamiWorkshopData (). Fig. 2 shows that the effect of diffusion on the water surface solution is clearly negligible. This can be explained by looking at the structure and values of in Fig. 3. With a water surface that is smooth almost everywhere, the dynamic diffusion coefficient is so small that its effect becomes minimal. We will see later that this will not be the case in problems with a greater degree of irregularity of the surface.
We plot the spacetime evolution of the shoreline in Fig. 4 where the wave elevation and total depth are plotted in the proximity of the shore. The dashed red curve represents the tabulated shoreline tsunamiWorkshopData (). By direct comparison with Carrier’s results carrierEtAl2003 (), the patterns of the water surface elevation () and total water depth are in full agreement.
6.2 Grid convergence rate
To measure the convergence rate of the model, we compare the computed solutions against the analytic solution of the flow in a onedimensional parabolic bowl gallardoEtAl2007 (); thacker1981 (). The parabolic topography is defined as:
where m and in m. The initial velocity is and the water surface begins to oscillate due to gravity only. The solution is computed using 8, 16, 32, 64, and 128 elements of order 4. The CG solution using 128 elements is plotted in Fig. 5, where we compare the inviscid solution against the viscous computation at s. The solution preserves its smoothness at all times as long as the resolution is sufficiently fine. At high resolution, it is evident that the CG approximation to this problem does not require stabilization. As the resolution is decreased by a factor of 4 (Fig. 6), the inviscid solution begins to develop oscillations although stability is not compromised. Unlike the inviscid one, the stabilized solution preserves the surface flatness although it moves slower due to a viscous source term in the momentum equation. This fact is also reflected in the computation of the normalized error norms plotted in Fig. 7. The analytic solution against which the error is calculated is defined for an inviscid flow. A difference between the viscous and inviscid computation is hence expected. The same observation applies to the CG and DG computations alike. The DG solutions to the same problem are plotted in Fig. 8 and 9.
6.3 2D oscillation in a symmetrical paraboloid
This test was defined by Thacker thacker1981 () who computed the exact solution of the problem. The paraboloid of revolution is defined as :
in the domain , with m. The initial condition is a reversed paraboloid (See top row of Fig. 11). The variation of the bed slope and the radial symmetry are a perfect test to assess the inundation scheme with two dimensional effects. The inviscid CG and DG solutions are plotted in Figs. 10 and 11, respectively, between and s. The radial symmetry of the solution is preserved throughout the simulation, and no spurious modes can be observed in the proximity of the moving shoreline. Although it is not fully visible from the plots, we observed a small discrepancy between the exact and computed solutions that occurs at the very center of the bowl at the latest time. More specifically, the numerical solution tends to move slightly slower than the exact solution along the bowl centerline. This occurs for CG and DG alike.
6.4 2D flooding problem in a channel with three mounds
This test kawaharaUmetsu1986 (); brufauEtAl2002 () was used by Xing and Zhang xingZhang2013 () and Gallardo et al. gallardoEtAl2007 () to assess, on unstructured grids, the same positivitypreserving method that we are applying here xingZhangShu2010 (). We reproduce their results here to verify its correct implementation in our code. The three mounds in the channel , are defined by the function where
(11)  
The flooding is triggered by a dam break at s. The evolution of the flooding from to s is shown in Fig. 12. We computed the solution using a relatively coarse grid made of elements of order 4. The results are in agreement with gallardoEtAl2007 () and xingZhang2013 (): the flow separation is well resolved and occurs after 10 seconds from the dam breaking. After 15 seconds, the first water front has reached the back wall and is fully reflecting back by 20 seconds. By 40 seconds, the flow has almost reached full rest.
6.5 2D solitary wave runup and rundown on a circular island
A solitary wave runup on a circular island was studied in briggsEtAl1995 () at the Waterways Experiment Station of the US Army Corps of Engineers. In this example, the initial wave is modeled via the following analytic definition by Synolakis synolakis1987 ():
(12) 
where m is the wave amplitude, m, m is the initial still water level, and
The island is a cone given as
where , m, and is centered at m. The cone is installed on a flat bathymetry. The fluid is confined within four solid walls.
To understand how the proposed diffusion and numerical approximation depend on the grid, we ran the simulation at the four resolutions m. Fig. 13 shows that the stabilized DG solution is converging to the same solution and the main features of the interacting waves are reproduced almost equally across the four grids. Certainly, the 40 cm grid spacing is the most dissipative, although it is encouraging to see how the important features resolved at 5 cm are still well represented on the coarsest grid. The same observation applies to the CG solution (plot not shown).
In Figs. 1417 we plot the projection of the 2D solution on the plane . The spurious modes that characterize the water surface in the proximity of the sharpest wave front are fully removed by DynSGS (Fig. 14) without weakening the front sharpness. We observe the same in the velocity field plotted in Fig. 15. This is in full agreement with the application of DynSGS to nonlinear wave problems with strong discontinuities, as previously shown in (marrasEtAl2015b, , Figs. 16, 17) for the solution of the Burgers’ equation. We briefly mentioned above how DG already has builtin dissipation. This is clearly visible from the plot of Fig. 16; the unstabilized DG solution shows no oscillations except for, at most, some minimal under and overshooting. This implies that its residual is so small that the effect of SGS eddy viscosity reduces to a minimum value. This explains why the inviscid and viscous DG solutions look similar.
It is important to notice the same space distribution of the CG and DG velocity waves shown in Figs. 15 and 17. This is a visual confirmation that the two solutions present the same dispersion properties, which is to be expected from two analogous numerical approximations to the same problem. This is indicative of a correct implementation of the unified CG/DG framework.
The solution of the top plot in Fig. 18, was computed with stabilization applied to the continuity ( in Eq. (1)) and momentum equations. When compared against the unstabilized solution (bottom plot in the same figure), we notice that the features of the fronts of the interacting waves are fully preserved. Furthermore, the fronts are not excessively smeared out as the high frequency modes are removed. The effect of momentum stabilization on the velocity field is shown in Fig. 19.
The instantaneous energy spectra of the viscous and inviscid CG and DG solutions at s are plotted in Fig. 20. The difference between the CG and DG curves is striking. The viscous and inviscid DG spectra overlap almost fully and show approximately the same decay across the entire spectrum, from a 5/3 slope in the inertial subrange to a 3 slope in the dissipation wave numbers (refer to boffettaEcke2012 () for a review on twodimensional flows and their energetics). This is only true as long as the resolution is not too coarse, especially so in the case of CG. At very coarse resolutions ( m), neither CG or DG can avoid energy from building up in the highest modes unless artificial diffusion is used. The inherent dissipation of DG is no longer sufficient.
We stated above that is only active where the equation residuals (i.e. gradients) are important. In the case of water waves, this occurs in the proximity of the wave fronts. In Fig. 21, we plot to show its spatial structure and its evolution between and seconds. This plots clearly show how diffusion is equally zero away from the fronts and only activates where really necessary. It may not be so obvious to achieve this by using an artificial diffusion that is not residualbased. To provide a visual correlation between and the wave features, in Fig. 22 we plot the stabilized spectral element solution of the water surface at a grid resolution m.
7 Conclusions
We presented the numerical solution of the shallow water equations via continuous and discontinuous Galerkin (CG/DG) methods to model problems involving inundation. Careful handling of the transition between dry and wet surfaces is particularly challenging for highorder methods. The most adopted solution in these cases consists of lowering the approximation order in either the transition zone alone, or in the whole domain. By extending to CG a simple slope limiter originally designed for DG xingZhangShu2010 (), and by combining it with a thin water layer always present in dry regions, we showed that the nominal highorder accuracy of the underlying space approximations was preserved and that the solution smoothness in the proximity of moving boundaries was maintained in one and twodimensions. The simplicity of this approach is effective for CG as well as DG. More important, we demonstrated its effectiveness using higher order elements. This limiter, unless embedded with a positivitypreserving scheme for water height, does not prevent the highorder solution from triggering unphysical Gibbs oscillations in the proximity of strong gradients. To overcome this problem, we presented a dynamically adaptive dissipation based on a residualbased subgrid scale eddy viscosity model (DynSGS). By numerical examples, we demonstrated the following properties of this model when solving the shallow water equations: It removes the Gibbs oscillations that form in the proximity of sharp wave fronts while preserving the overall accuracy and sharpness of the solution everywhere else. This is possible thanks to the residualbased definition of the dynamic diffusion coefficient. For coarse grids, it prevents energy from building up at small wavenumbers; this is very important to preserve numerical stability in the flow regimes we are interested in. The model has no user tunable parameter, which is of great advantage when the model is to be used by an external user. It is important to underline that the natural, builtin dissipation of DG be may be large enough that the contribution of DynSGS becomes irrelevant. When this happens, the dynamic dissipation detects it from the residual, and hence limits its own strength. Finally, a three stage, second order explicitfirststage, singly diagonally implicit RungeKutta (ESDIRK) time integration scheme was implemented to overcome the small timestep restriction incurred by highorder Galerkin approximations.
8 Acknowledgments
The authors would like to acknowledge the contribution of Haley Lane, who implemented the onedimensional version of the wetting and drying algorithm used in this work. The authors would also like to acknowledge Karoline Hood who tested the correctness of the implicit solver in her NPS Master’s thesis (karolineTHESIS2016 ()). The authors are also thankful to Oliver Fringer and the Rogers for discussions regarding coastal flows, and to Stephen R. Guimond for providing his MATLAB functions to compute the energy spectra. FXG acknowledges the support of the ONR Computational Mathematics program, and FXG and EMC acknowledge the support of AFOSR Computational Mathematics.
References
 (1) The third international workshop on longwave runup models: Benchmark problem #1: Tsunami runup onto a plane beach. http://isec.nacse.org/workshop/2004_cornell/bmark1.html. 2004.

(2)
D. Abdi, L. C. Wilcox, T. Warburton, and F. X. Giraldo.
A GPU accelerated continuous and discontinuous Galerkin
nonhydrostatic atmospheric model.
Submitted. Preprint available at
www.researchgate.net/publication/
290193830_A_GPU_
Accelerated_Continuous_and_Discontinuous_Galerkin_Nonhydrostatic_Atmospheric_Model, 2016.  (3) D. S. Abdi and F. X. Giraldo. Efficient construction of unified continuous and discontinuous Galerkin formulations for the 3D Euler equations. J. Comput. Phys., 320:46–68, 2016.
 (4) M. Ainsworth, P. Monk, and W. Muniz. Dispersive and dissipative properties of discontinuous Galerkin finite element methods for the secondorder wave equation. J. Sci. Comput., 27:5–40, 2006.
 (5) D. Arnold. An interior penalty finite element method with discontinuous elements. SIAM J. Numer. Anal., 19:742–760, 1982.
 (6) N. Beisiegel and B. Behrens. Quasinodal thirdorder Bernstein polynomioals in a discontinuous Galerkin model for flooding and drying. Environ. Earth Sci., 74:7275–7284, 2015.
 (7) G. Boffetta and E. Ecke. Twodimensional turbulence. Annu. Rev. Fluid Mech., 44:427–451, 2012.
 (8) J. P. Boyd. The erfclog filter and the asymptotics of the Euler and Vandeven sequence accelerations. In A.V. Ilin, L.R. Scott (Eds.), Proceedings of the Third International Conference on Spectral and High Order Methods, Houston Journal of Mathematics, pages 267–276, 1996.
 (9) M. J. Briggs, C. E. Synolakis, G. S. Harkins, and D. R. Green. Laboratory experiments of tsunami runup on a circular islad. Pageoph., 144:569–593, 1995.
 (10) P. Brufau, M. E. VázquezCendón, and P. GarcíaNavarro. A numerical model for the flooding and drying of irregular domains. Int. J. Numer. Methods Fluids, 39:247–275, 2002.
 (11) S. Bunya, E. J. Kubatko, J. J. Westerink, and C. Dawson. A wetting and drying treatment for the RungeKutta discontinuous Galerkin solution to the shallow water equations. Comput. Methods Appl. Mech. Engrg., 198:1548–1562, 2009.
 (12) J.C. Butcher and P. Chartier. The effective order of singlyimplicit RungeKutta methods. Numerical Algorithms, 20(4):269–284, 1999.
 (13) G. F. Carrier, T. T. Wu, and H. Yeh. Tsunami runup and drawdown on a plane beach. J. Fluid Mech., 475:79–99, 2003.
 (14) B. Cockburn and C W. Shu. TVB RungeKutta local projection discontinuous Galerkin finite element method for conservation laws II: general framework. Math. Comput., 52:411–435, 1989.
 (15) C Dawson and V. Aizinger. A discontinuous Galerkin method for threedimensional shallow water equations. J. Sci. Comp., 22:245–267, 2995.
 (16) A. J. C. de SaintVenant. Théorie du mouvement nonpermanent de eaux, avec application aux crues de rivières et à l’introduction des marées dans leur lit. C. R. Acad. Sc. Paris, 73:147–154, 1871.
 (17) O. Delestre, C. Lucas, P. A. Ksinant, F. Darboux, C. Laguerre, T. N. T. Vo, F. James, and S. Corier. SWASHES: a compilation of Shallow Water Analytic Solutions for Hydraulic and Environmental Studies. arXiv, 1110.0288v6 [math.NA], 2014.
 (18) D. Drikakis, H. Marco, F. F. Grinstein, C. R. DeVore, C. Fureby, M. Liefvendahl, and D. L. Youngs. Numerics for ILES: limiting algorithms. In F. F. Grinstein, L. G. Margolin, and Rider W. J., editors, Implicit LargeEddy Simulation: Computing Turbulent Flow Dynamics, pages 94–129. Cambridge University Press, 2007.
 (19) J. M. Gallardo, C. Parés, and M. Castro. On a wellbalanced highorder finite volume scheme for shallow water equations with topography and dry areas. J. Comput. Phys., 227:574–601, 2007.
 (20) R. Gandham, D. Medina, and T. Warburton. GPU accelerated discontinuous galerkin methods for shallow water equations. Commun. Comput. Phys., 18:37–64, 2015.
 (21) J. F. Gerbeau and B. Perthame. Derivation of the viscous SaintVenant system for laminar shallow water; numerical validation. Discrete Contin. Dyn. Syst. Ser. B, 1:89–102, 2001.
 (22) F. X. Giraldo. A spectral element shallow water model on spherical geodesic grids. Int. J. Num. Meth. Fluids, 35:869–901, 2001.
 (23) F. X. Giraldo, J. S. Hesthaven, and T. Warburton. Nodal highorder discontinuous Galerkin methods for spherical shallow water equations. J. Comput. Phys., 181:499–525, 2002.
 (24) F X. Giraldo, J F.. Kelly, and E. Constantinescu. Implicitexplicit formulations of a threedimensional Nonhydrostatic Unified Model of the Atmosphere (NUMA). SIAM J. Sci. Comput., 35:1162–1194, 2013.
 (25) F. X. Giraldo and M. Restelli. Highorder semiimplicit timeintegrators for a triangular discontinuous Galerkin oceanic shallow water model. Int. J. Numer. Methods Fluids, 63:1077–1102, 2010.
 (26) S. Gottlieb, C.W. Shu, and E. Tadmor. Strong stabilitypreserving highorder time discretization methods. SIAM review, 43(1):89–112, 2001.
 (27) R. Gourgue, O. Comblen, J. Lambrechts, T. Kärnä, V. Legat, and E. Deleersnijder. A fluxlimiting wettingdrying method for finiteelement shallowwater models, with application to the Scheldt Estuary. Adv. Water Resourc., 32:1726–1739, 2009.
 (28) O. Guba, M. A. Taylor, P. A. Ullrich, J. R. Overfelt, and M. N. Levy. The spectral element method on variable resolution grids: evaluating grid sensitivity and resolutionaware numerical viscosity. Geosci. Model Dev. Discuss., 7:4081–4117, 2014.
 (29) J L. Guermond and B. Popov. Viscous regularization of the Euler equations and entropy principles. SIAM J. Appl. Math., 74(2):284–305, 2014.
 (30) E. Hendricks, M. Kopera, and F. X. Giraldo. Evaluation of the utility of static and adaptive mesh refinement for idealized tropical cyclone problems in a spectral element shallow water model. Mon. Wea. Rev. (in press), 2016.
 (31) K. Hood. Modeling storm surge using discontinuous Galerkin methods. Masters Thesis, Dept. of Applied Maths, Naval Postgraduate School, Monterey, CA, U.S.A., 2016.
 (32) M. Iskandarani, D .B. Haidvogel, and J. P. Boyd. A staggered spectral element models with application to the oceanic shallow water equations. Int. J. Numer. Methods Fluids, 20:393–414, 1995.
 (33) T. Kärnä, B. de Brye, O. Gourgue, J. Lambrechts, R. Comblen, V. Legat, and E. Deleersnijder. A fully implicit wettingdrying method for DGFEM shallow water models, with an application to the Scheldt Estuary. Comp. Methods Appl. Mech. Engrg., 200:509–524, 2011.
 (34) M. Kawahara and T. Umetsu. Finite element method for moving boundary problems in river flow. Int. J. Numer. Meth. Fluids, 6:365–386, 1986.
 (35) C.A. Kennedy and M.H. Carpenter. Additive RungeKutta schemes for convectiondiffusionreaction equations. Appl. Numer. Math., 44(12):139–181, 2003.
 (36) D.A. Knoll and D.E. Keyes. Jacobianfree Newton–Krylov methods: a survey of approaches and applications. Journal of Computational Physics, 193(2):357–397, 2004.
 (37) M. A. Kopera and F. X. Giraldo. Analysis of adaptive mesh refinement for IMEX discontinuous Galerkin solutions of the compressible Euler equations with application to atmospheric simulations. J. Comput. Phys., 92117:275, 2014.
 (38) M. A. Kopera and F. X. Giraldo. Mass conservation of the unified continuous and discontinuous elementbased galerkin methods on dynamically adaptive grids with application to atmospheric simulations. J. Comput. Phys., 297:90–103, 2015.
 (39) E. J. Kubatko, J. J. Westerink, and C. Dawson. hp discontinuous Galerkin methods for advection dominated problems in shallow water flow. Comp. Methods Appl. Mech. Engrg., 196:437–451, 2006.
 (40) D. K. Lilly. On the numerical simulation of buoyant convection. Tellus, 14:148–172, 1962.
 (41) H. Ma. A spectral element basin model for the shallow water equations. J. Comput. Phys., 109:133–149, 1993.
 (42) S. Marras, M. Nazarov, and F. X. Giraldo. Stabilized highorder Galerkin methods based on a parameterfree dynamic SGS model for LES. J. Comput. Phys., 301:77–101, 2015.
 (43) C. Michoski, C. Dawson, E. J. Kubatko, D. Wirasaet, S. Brus, and J. J. Westerink. A comparison of artificial viscosity, limiters, and filters, for high order discontinuous Galerkin solutions in nonlinear settings. J. Sci. Comput., 66:406–434, 2016.

(44)
A. Müller, M. A. Kopera, S. Marras, L .C. Wilcox, T. Isaac, and F. X.
Giraldo.
Strong Scaling for Numerical Weather Prediction at Petascale with
the Atmospheric Model NUMA.
Submitted to the IEEE International Parallel and Distributed
Processing Symposium. Preprint available at
www.researchgate.net/publication282887372_
Strong_Scaling_for_Numerical_
Weather_Prediction_at_Petascale_with_the_Atmospheric_Model_NUMA, 2016.  (45) R D. Nair, S J. Thomas, and R D. Loft. A discontinuous Galerkin global shallow water model. Mon. Wea. Rev., 133, 2007.
 (46) M. Nazarov and J. Hoffman. Residualbased artificial viscosity for simulation of turbulent compressible flow using adaptive finite element methods. Int. J. Numer. Methods Fluids, 71:339–357, 2013.
 (47) R. Pasquetti, J.L. Guermond, and B. Popov. Stabilized spectral element approximation of the Saint Venant system using the entropy viscosity technique. In Spectral and High Order Methods for Partial Differential Equations ICOSAHOM 2014, Lecture Notes in Computational Science and Engineering, pages 397–404, 2015.
 (48) C. Phan Van, E. Deleersnijder, D. Bousmar, and S. Frazão Soares. Simulation of flow in compound openchannel using a discontinuous Galerkin finiteelement method with Smagorinsky turbulence closure. J. Hydroenv. Research, 8:396–409, 2014.
 (49) N. Rakowsky, A. Androsov, A. Fuchs, S. Harig, A. Immerz, S. Danilov, W. Hiller, and J. Schröter. Operational tsunami modelling with TsunAWI  recent developments and applications. Nat. Haz. Earth Sys. Sci., 13:1629–1642, 2013.
 (50) W. J. Rider. The relationship of MPDATA to other highresolution methods. Int. J. Numer. Methods Fluids, 50:1145–1158, 2006.
 (51) P. Sagaut. Large Eddy Simulation for Incompressible flows. Springer, 3rd edition, 2006.
 (52) C W. Shu. TVB uniformly highorder schemes for conservation laws. Math. Comput., 49:105–121, 1987.
 (53) J. Smagorinsky. General circulation experiments with the primitive equations: I. the basic experiement. Mon. Wea. Rev., 91:99–164, 1963.
 (54) P. K. Smolarkiewicz. A simple positive definite advection scheme with small implicit diffusion. Mon. Wea. Rev., 111:479–486, 1983.
 (55) P. K. Smolarkiewicz. A fully multidimensional positive definite advection transport algorithm with small implicit diffusion. J. Comput. Phys., 54:325–362, 1984.
 (56) C. E. Synolakis. The runup of solitary waves. J. Fluid Mech., 185:523–545, 1987.
 (57) M. Taylor, J. Tribbia, and M. Iskandarani. The spectral element method for the shallow water equations on the sphere. J. Comput. Phys., 130:92–108, 1997.
 (58) W. C. Thacker. Some exact solutions to the nonlinear shallowwater wave equations. J. Fluid Mech., 107:499–508, 1981.
 (59) A. Uranga, P.O. Persson, M. Drela, and J. Peraire. Implicit Large Eddy Simulation of transition to turbulence at low Reynolds numbers using a discontinuous Galerkin method. Int. J. Numer. Methods Engrg., 87:232–261, 2011.
 (60) H. Vandeven. Family of spectral filters for discontinuous problems. J. Sci. Comp., 6:159–192, 1991.
 (61) S. Vater, N. Beisiegel, and J. Behrens. A limiterbased wellbalanced discontinuous Galerkin method for shallowwater flows with wetting and drying: one dimensional case. Adv. Water Res., 85:1–13, 2015.
 (62) Y. Xing and X. Zhang. Positivitypreserving wellbalanced discontinuous Galerkin methods for the shallow water equations on unstructured triangular meshes. J. Sci. Comput., 57:19–41, 2013.
 (63) Y. Xing, X. Zhang, and C. W. Shu. Positivitypreserving high order wellbalanced discontinuous Galerkin methods for the shallow water equations. Adv. Water Res., 33:1476–1493, 2010.