An adaptive discontinuous Galerkin method for the simulation of hurricane storm surge

Numerical simulations based on solving the 2D shallow water equations using a discontinuous Galerkin (DG) discretisation have evolved to be a viable tool for many geophysical applications. In the context of flood modelling, however, they have not yet been methodologically studied to a large extent. Systematic model testing is non-trivial as no comprehensive collection of numerical test cases exists to ensure the correctness of the implementation. Hence, the first part of this manuscript aims at collecting test cases from the literature that are generally useful for storm surge modellers and can be used to benchmark codes. On geographic scale, hurricane storm surge can be interpreted as a localised phenomenon making it ideally suited for adaptive mesh refinement (AMR). Past studies employing dynamic AMR have exclusively focused on nested meshes. For that reason, we have developed a DG storm surge model on a triangular and dynamically adaptive mesh. In order to increase computational efficiency, the refinement is driven by physics-based refinement indicators capturing major model sensitivities. Using idealised numerical test cases, we demonstrate the model’s ability to correctly represent all source terms and reproduce known variability of coastal flooding with respect to hurricane characteristics such as size and approach speed. Finally, the adaptive mesh significantly reduces computing time with no effect on storm waves measured at discrete wave gauges just off the coast which shows the model’s potential for use as a robust simulation tool for real-time predictions.


Introduction
Extreme storms and floods are among the most devastating natural disasters worldwide faced by inhabitants of coastal areas and those whose livelihood depends on marine resources alike (see, e.g. (Lin et al. 2010)). Their accurate simulation and prediction are vital for hazard assessment and the implementation of mitigation strategies such as the opening and closing of flood barriers in order to protect local populations and property. Due to the time critical nature of forecasts, current operational simulation environments that form the basis of flood forecasts are chosen for their robustness and optimised run time. However, this efficiency comes at the cost of having to use reduced hydrodynamical model approximations for the underlying processes and simplified parameterisations for physical forcing terms. This means that operational models often do not model the full physics (Tull 2018). Instead, such simple inundation models simply raise water levels to a given constant elevation. Linear approaches like this are computationally inexpensive and work well for generalised scenarios and at broad geographic scale, but as we become more specific and local in the attempted projections, the complexities increase, and non-linearities become more significant in the behaviour of the modelled physical processes. Therefore, we cannot use linear mathematical approaches to precisely model these situations. Non-linear models, on the other hand, can be costly, and the high-order information they contain might require advanced filtering to guarantee numerical robustness and efficiency. These more accurate models, however, are of critical importance. As was shown in (Barnard et al. 2019), the accuracy of projections for future events, in particular with respect to our changing climate, depends on the accurate modelling of small scale physical dynamics at large scale.
In this study, we focus on discontinuous Galerkin (DG) methods to solve the 2D non-linear shallow water equations. These depth integrated equations are computationally less demanding than more accurate 3D Navier-Stokes equations and have been shown to yield good results for coastal applications. Since realistic measurement data often lacks high-order smoothness, and we aim to build a model that will eventually be useful for practical applications, we furthermore concentrate on piecewise linear polynomials to represent our quantities of interest. This leads to a formally second order in space accurate numerical scheme which we believe to be sufficient to be useful. DG methods have recently gained a lot of attention for geophysical applications because of their local mass conservation and geometric flexibility. Moreover, numerical computations are performed locally since elements only communicate over element interfaces (edges) through the computation of numerical fluxes. This is a particular advantage for current state-of-the-art computational architectures and although beyond the scope of this study, we remark that the local nature of this discretisation makes it ideal for parallelisation as shown for example in (Abdi et al. 2019).
Due to their relative computational expense, current operational flood models do not employ DG discretisations; they are based on other numerical techniques (see also Danilov 2013). A common choice is continuous finite elements (Piggott et al. 2008a), because they work well with nonuniform meshes. However, they only yield global mass conservation and parallelisation is more difficult because large linear systems often need to be inverted. The majority of computer models used in practice rely on finite difference discretisations. Those have the advantage of discretisations for parameterisations of source terms such as bottom friction being existent and well established. We remark, though, that in theory, the implementation of parameterisations in a DG framework would be easy as it offers the opportunity to control in cell values as opposed to point values only and it does not formally require differentiability of the source term. In addition to representing complex coastlines by the use of curvilinear meshes-an approach that has been successfully used in finite difference frameworks-DG methods offer the possibility to represent varying geometries directly through highorder approximations within elements. Another common choice for flood modelling are finite volume (FV) methods which are conservative and work well even on unstructured grids as shown in (Danilov 2013). First-order DG methods are equivalent to FV methods, so that DG methods can be seen as one possibility to extend FV methods to higher-order accuracy without requiring high-order reconstructions that might be computationally expensive.
A higher computational efficiency of mesh-based numerical methods such as DG methods can be achieved by using nonuniform meshes (Piggott et al. 2008b) or dynamical adaptive mesh refinement as in , leading to an improved performance on current computing architectures. DG methods have been successfully employed in combination with a non-uniform but static mesh to hindcast the coastal flooding caused by hurricane Ike (Dawson et al. 2011). Ideally, and to save computational time, the use of dynamically adaptive nonuniform meshes would be preferable. These meshes are usually driven by either heuristic (or physics-based) refinement indicators or error estimators that are based on model sensitivities which in turn can be derived from solving adjoint equations, as in (Farrell et al. 2013) for a finite element framework. Solving adjoint equations, however, can be computationally expensive as it requires the solution of a different set of equations backwards in time with coefficients that result from the forward solution of the system. This increases memory requirements significantly. Additionally, Beckers et al. (2019) show that the discontinuous nature of DG methods poses systematic difficulties on deriving appropriate adjoint equations. The additional flexibility of a dynamically adaptive mesh would decrease the dependence of a high level of mesh optimisation that is particular for the geographic region of interest. An approach using finite volumes on a quadrilateral nested mesh has been used in one previous study (Mandli and Dawson 2014), but the full potential of dynamically adaptive mesh refinement using physics-based refinement indicators for fully adaptive meshes is yet to be explored. Hence, this study aims at investigating DG methods for flood applications on a dynamically adaptive triangular mesh.
The accurate modelling of inundation using the shallow water equations is mathematically challenging as the interface between water and land becomes a moving boundary, and the theoretical validity of the underlying equations breaks down in near-dry regions. Recently developed numerical methods show improved robustness due to improved limiting and filtering techniques (Vater et al. 2015;Vater et al. 2019). We will adopt the novel limiting strategy presented in (Vater et al. 2019) which features a velocity-based reconstruction of the momentum and allows us to compute meaningful velocities even close to the wet/dry interface.
The major driving forces of storm surges are extreme winds and pressure gradients. This study considers hurricane storm surges in Sect. 4.1 and employs the cyclonic wind model by (Holland 1980) to compute continuous wind fields as well as corresponding atmospheric pressures. Finally, we implement all source terms that are relevant for coastal flooding and storm waves and show that in combination with the dynamically adaptive mesh refinement as proposed in , we built a new model that is suitable for coastal storm surge modelling.
In summary, in this study, we combine a DG model with a dynamically adaptive mesh that is driven by physics-based refinement indicators. We show that a recently developed advanced limiting strategy to treat wetting and drying as well as the discrete implementation of all relevant source terms lead to robust and accurate results on this mesh. This paper is organised as follows. Section 2 outlines the numerical model and summarises the implementation of all relevant source terms, and Sect. 2.2 gives a brief introduction to the adaptive mesh refinement used later on. The results section is separated into two sections: Section 3 presents a number of numerical test cases ranging from analytical to idealised yet realistic scenarios that showcase the model's inundation stability, conservation properties, as well as the robustness of the wind forcing. In Sect. 4, we then investigate the described model using the adaptive mesh refinement from Sect. 2.2 to demonstrate its suitability for idealised storm surge modelling on a dynamically adaptive mesh. The Sect. 5 gives an outlook for future work and discusses implications, shortfalls and limitations of the study.

A discontinuous Galerkin discretisation of the full 2D non-linear shallow water equations
The system of equations of interest for this study are depthintegrated shallow water equations in two dimensions which can be written in flux form where the prognostic variables are U = (h, hu) ⊤ : the water depth h and the 2D momentum hu with a 2D velocity u = (u, v) ⊤ defined on Ω × T, with Ω ⊂ ℝ 2 and T a finite time interval. Spatial coordinates are denoted as x = (x, y) ⊤ ∈ Ω. The partial temporal derivative is denoted by ∂U ∂t and ∇⋅:¼ ∂ ∂x ; ∂ ∂y ⋅ is the 2D divergence operator. The flux F and source term S are defined as where g = 9.81 m s −2 is the acceleration due to gravity and I 2 is the 2 × 2 identity matrix. We denote the temporally constant bathymetry by b = b(x) and introduce the notation H(x, t) = h(x, t) + b(x) for the total height (see also Fig. 1). In addition to the influence of bathymetry, the source term S comprises a Coriolis forcing τ C , a vector-valued wind stress τ W and the water density ρ, which we will assume to be constant throughout this study, as well as the atmospheric pressure p A and a bottom friction τ B which are described in more detail in Sect. 2.1. Throughout this paper, vector-valued quantities are indicated by a bold print while all other quantities are assumed to be scalar. It is well understood that the shallow water approximation is somewhat over simplified for applications of coastal inundation, where non-hydrostatic effects may very well play a significant role. On the other hand, many practical applications involve large uncertainties in data and parameterisations of small-scale effects can be calibrated to yield practically useful results. Furthermore, a very effective extension of the shallow water solver by nonhydrostatic correction following Casulli and Stelling (1998) can be applied to our discretisation (see Jeschke et al. 2017).
We discretise Eq. (1) using a discontinuous Galerkin approach comprised of (a) decomposing the domain Ω = ∑ i Ω i into triangles, (b) approximating U = ∑ k U k (t)ϕ k (x) by linear Lagrange polynomials locally in each triangle and (c) integrating locally in space against test functions. Our test functions are linear Lagrange polynomials, so that the resulting semidiscrete system reads In Eq. (3), F * is a numerical Rusanov flux at the cell interfaces defined as where U L and U R are the prognostic variables on the left and right of the respective edge and λ max ¼ max jujþ ð ffiffiffiffiffi gh p ; jvj þ ffiffiffiffiffi gh p Þ is the maximum local directional shallow water wave speed, see also (Toro 2009). As demonstrated in (Toro 2009), other choices to compute F * are available. In this study, we chose the Rusanov flux for its computational simplicity and its slightly higher dissipation compared to the computation of the exact solution of the Riemann problem (see (LeVeque 2002)) which we believe adds to the robustness of the presented DG method. Note that we integrated the flux integral by parts twice to obtain the often called strong form (Hesthaven and Warburton 2008). This form has desirable properties with respect to well-balancing as shown in (Beisiegel 2014) and elaborated in more detail in Sect. 2.1.3. We solve the integrals in Eq. (3) with an exact, i.e. third-order, interpolatory Gauss quadrature with corresponding Gauss-Legendre points in order to not introduce lower order numerical errors into the scheme. This allows to extract the time derivative and, after re-organisation, write the system (Eq. (3)) as a system of equations in t of the form where H denotes the discretised flux and source terms. This system (Eq. (4)) can be solved using a strong-stability preserving (SSP) multi-stage Runge-Kutta method provided suitable boundary conditions. For hyperbolic problems, this type of time integrator is preferred because it maintains stability properties of the continuous equations such as total variation stability; hence, they prevent spurious energy build-up even in the presence of shocks (see for example (Gottlieb et al. 2011)).
Due to superior stability properties demonstrated in (Kubatko et al. 2008), we used a three-stage, second-order Runge-Kutta method (RK23) for this study: where a slope limiter is applied to the prognostic variables after every intermediate step. The slope limiter we are using is velocity-based and non-destructive with respect to wellbalancing and non-negativity preservation. It modifies (h, hu) ⊤ in elements where artificial gradients occur and sets h = 0 for small water depth h < ϵ tol m, provided a userdefined and test case-dependent parameter ϵ tol . This will remove oscillations within elements that are typically treated by using dynamic viscosity as in (Marras et al. 2016); hence, an additional viscosity source term is not needed for the presented test cases. More detail can be found in (Vater et al. 2019) where the limiter was successfully applied to tsunami benchmark test cases and to model flood scenarios. In the present study, we show that it can be used as well for accurate shore line modelling and the reduction of spurious oscillations on adaptive meshes even when a variety of source terms are involved. Hyperbolic problems are time step restricted and numerical stability of methods for their explicit solution depends on the Courant-Friedrich-Levy (CFL) condition (see Courant et al. 1928). Using a maximum possible value for a stable CFL number-we found cfl max = 0.3 a good choice for linear DG elements-and using the three-stage Runge-Kutta method described above, we can compute a stable time step by where c max ¼ max is the global maximum shallow water velocity, with k running over all degrees of freedom. The quantity Δx min ¼ min k Δx min k is an estimate of the shortest height Δx min k of a triangle that is computed using the volume formula for triangles at the current time step leading to a global time step Δt. This is done before each time step allowing for a temporally changing and maximum possible global time step.

Implementation of source terms
Coastal flooding is a problem involving the interplay between many source terms. The source terms directly impact the momentum and thus the free surface. For reasons of comprehensiveness, we will give detail of their implementation in the following subsections.

Earth's rotation
Hurricanes can cover large areas up to 100s or 1000s of km 2 . On large geographic scale, Earth's rotation has a nonnegligible influence on water circulation. In the presented model, this Coriolis forcing τ C is of the form τ C = f(−hv, hu) ⊤ and f = 2ω sin (φ) with φ the latitude and ω = 7.2921 · 10 −5 rad s −1 the rotation rate of the Earth. A common approximation of the Coriolis force is the β-plane approximation, i.e. we linearly approximate f = f 0 + βy, where f 0 ∈ ℝ, y is a planar coordinate and β≈ 2ω r E is a constant that depends on the Earth's rotation and radius r E . We will use both parameterisations, the conventional and the β-plane approximation in this study.

Bottom friction
The bottom friction τ B is assumed to take on the form of a quadratic Manning law τ B ¼ gn 2 ∥u∥ 2 h 7=3 hu, where n is a dimensionless roughness parameter that is directly related to the nature of the bed and will take on values between 0.001 and 0.01 depending on the specifics of the test case. The bottom friction depends on the prognostic variables U. In coastal areas, numerical flow directions might be reversed during wave run up due to small fluid depths and a resulting large friction term. To ensure stability of the numerical result, we follow the split-implicit time stepping procedure in (Kesserwani and Liang 2012) to ensure that flow directions do not reverse within one time step. This means that friction terms for h < ε tol m and ∥hu∥ 2 < ε tol2 m 2 s −1 are neglected, where ε tol is the wet/dry tolerance of the slope limiter mentioned in Sect. 2 and ε tol2 is another small parameter. We chose ε tol ∈ [10 −4 , 10 −2 ], and ε tol2 = 10 −8 throughout this study. In the first momentum equation for hu, we then compute for every degree of freedom Then, if ΔtSf/D < | hu|, the flow will not reverse within that time step and we can compute the friction term as ΔtSf/D; otherwise, the friction is set to 0.
Analogously, the bottom friction can be computed for the second momentum equation by interchanging the roles of u and v.

Bathymetry and well-balancing
Non-linear interaction with bathymetry plays a crucial role in wave amplification. Special care has to be taken to prevent spurious waves caused by erroneous numerical approximations of bathymetry gradients. In this study, we solve the strong form (Eq. (3)) of the equations, which plays a crucial role in the following discussion. Well-balancing, i.e. preservation of the steady state at rest, is a desirable property of numerical discretisations. The steady state at rest assumes bathymetry gradients to be the only present source term. While well-balancing is achieved by setting gravity to 0 in cells adjacent to the wet/dry interface (Vater et al. 2019), the strong form is well-balanced on a discrete level without any modification when the flux divergence ∇ · F(U) is discretised after differentiation, since in general where F div (U k ) and F(U k ) are the nodal coefficients of the divergence of the flux and the flux, respectively. Since edgebased terms are always balanced, we show that for every wet element Ω i we obtain for u = 0. Equation (7) is fulfilled as long as the flux divergence and the source evaluated at the Lagrange points are balanced. This will always be the case in fully wet cells for the water at rest where u = 0, as always holds true for all degrees of freedom k. This balance is not achieved if the order of differentiation and discretisation of ∇ · F(U) are reversed. This can be easily shown by a counter example on a master element with edges (−1, −1), (−1, 1), (1, −1) and using linear Lagrange polynomials.

Wind fields and wind drag coefficients
The wind stress is τ W = γ τ τ with a wind friction γ τ ∈ ℝ + that models the energy transfer from the atmosphere to the ocean.
Given an external wind field v, the wind stress can be computed as where c d ∈ ℝ is a drag coefficient and ρ a is the air density which we assume to be 1.15 kg m −3 . The dimensionless drag coefficient c d = c d (v) depends on wind speeds. For hurricane models, several different wind drag parameterisations have been explored and commonly used ones are depicted in Fig.  2 (left). In Garrat (1977) (blue line), observations from the past 10 years are used to show that for absolute wind speeds |v| within a range of 4 to 21 m s −1 the drag coefficient fulfils a linear relationship c d · 10 3 = 0.75 + 0.067 |v| or a power law of the form 0.51|v| 0.46 . We note that the understanding is that the drag coefficient is dimensionless and units are removed from occurring velocities |v|. In Weisberg and Zheng (2006) (gray line), a drag coefficient as in (Large and Pond 1981) is used to study sensitivities with respect to approach speed, direction of approach and landfall location (see also Sect. 4.1 of this manuscript). Finally, Powell (2007) discusses wind drags for more extreme winds |v| > 50 m s −1 and numerically computes drag coefficients that exceed those known in the literature. Their new drag coefficients show improved results for practical applications.

Atmospheric pressure
The atmospheric pressure gradient ∇p A models part of the influence of the atmosphere onto the water column. In areas of relatively low air pressure in comparison to the ambient air pressure, this causes the water surface to slightly bulge upwards, increasing h in this area. The influence of the pressure gradient is non-negligible: Numerical observations show that about 10% of the resulting surge can be attributed to the pressure gradient. It is implemented in a straight forward manner, utilising the local derivatives of the Lagrange polynomial expansion. The contribution of the gradient of the atmospheric pressure to the source integral can be written as where ξ k are the Gauss-Legendre quadrature points, x l are Lagrange interpolation points and the matrix multiplication inside the parentheses projects the derivatives of the basis functions ϕ l onto quadrature points ξ k . Note that the model projects every triangle onto a master triangle, so that computationally, these matrix entries are constants and the matrix multiplication only has to be done once.

Adaptive mesh refinement
The computational model as described at the beginning of this Sect. 2 uses the mesh generator amatos  to create dynamically adaptive and conforming triangular meshes. Smaller triangles are obtained by bisection (Rivara 1984). While the numerical method works equally well on truly unstructured meshes (as demonstrated in (Vater et al. 2019)), the more structured hierarchical mesh refinement approach has a number of computational advantages that are of interest here. A depth-first traversal of the refinement tree gives rise to a space filling curve that allows for optimisation of the memory layout as well as computational organisation for optimised parallel computations (see, e.g. Behrens 2005; Behrens and Bader 2009). The dynamic grid manipulation involves problemdependent refinement indicators η Ω i for each element Ω i (Behrens 2006) to control the element-wise refinement or coarsening. A commonly used example for a useful indicator is the gradient of the total water height at time t: For applications such as idealised hurricane storm surge, see also Sect. 4.1 of this manuscript, a composite refinement indicator might be necessary that refines areas of bathymetry gradients as well as areas of strong winds. The indicator we used for this purpose is defined as with u = (u, v) ⊤ the horizontal velocity vector and Á j j Ω i ;1 ; Á j j Ω i ;2 the discrete L 1 , and L 2 norm on element Ω i respectively.
The width of elements is controlled by user-defined mesh levels λ ref and λ crs with λ crs ≤ λ ref . Starting from a macrotriangulation, the mesh is uniformly refined λ crs times until the coarse mesh level is reached. Subsequently, the refinement indicator is computed to determine elements that need to be flagged for refinement or coarsening. Moreover, user-defined tolerances 0 ≤ θ crs < θ ref ≤ 1 determine the fraction of the domain to be modified as follows: η Ω i t ð Þ the maximum value of the refinement indicator over all elements at time t. This manipulation strategy is carried out in a loop until the desired finest mesh level λ ref is reached. This approach will produce meshes that consist of small elements on mesh level λ ref in the area of interest with a continuous but narrow transition zone to coarse elements on mesh level λ crs outside the area of interest.
The node values are then interpolated or restricted after modification using the known Lagrange basis functions for each element. In regions close to the waterline, we keep the mesh relatively fine (see also Sect. 4.1) to avoid well-balancing issues. The refinement process has the further advantage of using a cache-efficient space-filling curve-ordering of elements (Behrens and Bader 2009), which allows fast access of neighbouring elements: This feature is particularly beneficial for local numerical discretisations such as discontinuous Galerkin method since elements only communicate over edges.
For convenience, the meshes are kept conforming, i.e. free of hanging nodes, throughout the simulation. We stress that

Results I: Source terms
In this section, we present test cases to show the major functionalities of the model. Note that physical dimensions are not explicitly given, but we assume throughout this paper that we are using SI base units. In some cases, we chose to not use SI base units to increase readability. In those cases, the dimensions are explicitly mentioned. The following subsections will in particular focus on the aspects of: -Inundation stability and conservation properties (Sect. 3.1) -Balanced wind and pressure forcing (Sect. 3.2) -Robustness of wind forcing (Sect. 3.3)

Tidal flow in a symmetric embayment
Here, we show one configuration of the simulation of an embayment first presented in Ip et al. (1998). Similar simulations have been reported on in Luettich and Westerink (1995) using the finite element model ADCIRC as well as in Dawson et al. (2011) using another discontinuous Galerkin model and a modified bathymetry. This simulation is to show idealised estuarine flooding by an in-and outgoing tide to show inundation stability and conservation properties of the model.
On a square domain Ω = [−L, L] 2 with L = 1500, a symmetric bathymetry as shown in Fig. 3 where we use the definition x b ¼ 10 −3 x; b y ¼ 10 −3 y: Moreover, the water surface is assumed to be in a steady state at rest. On the right boundary of the domain (at x = L), a tidal forcing is prescribed as  (Ip et al. 1998) and h de denotes the water depth in the deepest end of the domain. All other boundaries are set to be impermeable. Friction is assumed to be negligible, i.e. n = 0. For the initial water surface which is at rest, the water surface is levelled with the highest point of the bathymetry b max = 3.0.
To avoid polluted results caused by the definition of the boundary condition, we artificially increased the domain by L to the right, i.e. for b x ≥ 1:5, and assume the bathymetry in this area to be identically 0, hence adding a discontinuity to the bathymetry. As we will see, this did not pose a problem to our robust numerical method. We ran five full tidal cycles with a time step of Δt ≈ 4, a spatial resolution of Δx = 93.75 and a CFL number of 0.3 which we found to be the largest possible CFL number for this particular test case that allowed for a stable computation and which is similar to the maximal CFL numbers found in (Xu et al. 2014) for similar applications. The determination of theoretical bounds for CFL numbers for multi-dimensional problems is difficult, and most results are restricted to one space dimension (Kubatko et al. 2014). Snapshots of the numerical solutions over a cross section at y = 500 are found in Fig. 4 which shows that wetting and drying is occurring on the left boundary of the domain and, hence, that the inundation scheme is stable. We observe-see left display in Fig. 5-global mass conservation in agreement with (Ip et al. 1998) showing that the limiter presented in (Vater et al. 2019) is non-destructive with respect to global mass conservation and that the variation of total fluid volume due to the tidal in/outlet is of the order of O 10 7 À Á . We note that the model used in (Ip et al. 1998) achieves a much larger time step of Δt = 111.78 at a spatial resolution of Δx = 75. We attribute this to the fully implicit evolution scheme that their model uses.
Furthermore, we numerically computed the integral of the flux Hu over cross sections at x ∈ {±1200, ±600, 0}: The results are shown in the right display of Fig. 5. We see that after a short calibration time, the results match the results in (Ip et al. 1998) well and yield values for the cumulative transport of the order 10 3 . We note that especially at b x ¼ −1:20, we still get meaningful and non-zero results. This is notable because that part of the domain is close to the left hand boundary and therewith exposed to wetting and drying (see Fig. 4 for a close-up) and shows that the slope limiting strategy is gentle enough to reduce spurious oscillations while also yielding realistic values for velocities. Finally, we point out that the results are found to not be sensitive to the parameter ϵ tol that determines a cut-off value for small water depths. We ran the simulation with 10 −6 ≤ ϵ tol ≤ 10 −3 and did not find different results.

Steady-State Wind Test Case
This test case was described in (Davis and Sheng 2003). In a rectangular basin, Ω = [0, L] × [0, D] with dimensions L = 21000 and D = 5000 with constant bathymetry b(x) = 0 and impermeable boundaries on all four edges; a linear water surface h d + ξ is described as shown in Fig. 6. Here, h d = 5 and with g the acceleration due to gravity, ρ = 10 3 kg m −3 , the water density, and absent initial velocities u = 0. A constant wind stress is assumed as τ = (τ 1 , τ 2 ) ⊤ = (0.1,0) ⊤ and γ τ = 1 (for the definition of γ τ see also 2.1.4). In previous studies (Davis and Sheng 2003;Qiang et al. 2016), the steady state is described as shown in Table 1. The main balance here is to be achieved in the x-momentum equation between the flux divergence ∇ · F(U) and the wind stress, i.e.
on a discrete level. Here we used the symbol ¼ ! to indicate that the balance is to be enforced. The partial derivative of ξ can be computed as ∂ ∂x ξ x ð Þ ¼ τ 1 gh d ρ ð Þ −1 as ξ is a linear function of x. Substituting in the derivative of ξ and τ 1 = 0.1 in Eq. (10), we get after division by g: This equality is approximately satisfied because ξ ≪ h d by definition, i.e. h d þξ h d ≈1. For the simulation, we used a variable time step with a CFL number of 0.3 which is close to the theoretical maximum and was found to lead to stable results and a Runge-Kutta time stepping scheme of second order (RK23) until T end = 10 4 with  Table 1 and Fig. 7. In Table 1, we compare point values at t = T end with the steady state solution h ana and find them to be identical. Figure 7 considers the relative errors in surface height h at the numerical wave gauges at x ∈ {500, 10,500, 20,500} over time in more detail. We observe that numerical errors are of the order O 10 −5 À Á and that they stay bounded over time. We remark that although possible in this DG method, we did not consider physical viscosity in these simulations. This is in contrast to the findings in (Davis and Sheng 2003) who solved RANS equations with horizontal and vertical eddy viscosities to achieve an exact solution. The reason for this is that their model allowed the definition of an eddy viscosity of the form A H Δu and a boundary condition of the form ∂ z A v ∂ z u at the free surface. By setting A v u z = τ/ρ, they achieve an exact balance of the wind field and the gradient leading to an exact solution. Overall, we see that the steady-state values from Table 1 are matched and the wind and gravity forces are balanced. Our observed relative errors are of the order O 10 −5 À Á and are decreasing with increasing spatial resolution. Figure 8 furthermore shows snapshots of the numerical solution over the cross section at y = 2500 at discrete time steps. Shown are the fluid height and velocities for increments of 500 in light grey and for t = 0, T end /2, T end in colour. We observe that due to the boundary conditions, small meridional velocities u = (u, v) ⊤ form at the walls at (x = 0 and x = 21000). However, they stay bounded throughout the simulation, do not influence the surface height h significantly and do not destroy the steadystate balance. We attribute this to enough numerical diffusion being present in the model to disperse energy and prevent build-up.

Wind-induced circulation in a semi-enclosed homogeneous, rotating basin
This test case is described in (Sanay and Valle-Levinson 2005). In a semi-enclosed rectangular domain Ω = [0, 10D] × [0, D] with D = 10,000, a piecewise linear bathymetry (see also the sketch in Fig. 9) of depth h max ∈ ℝ is prescribed as follows: Here, x = (x, y) ⊤ is the spatial coordinate and h 0 = 3.0 is the minimum water depth. The water surface is at rest at time t = 0 and a constant wind stress τ = (τ C , 0) ⊤ , τ C ∈ ℝ, aligned with the xaxis is prescribed and linearly ramped up over a period of 6 hours. Six different configurations of the parameters, which are given in Table 2, were tested to assess how the maximal occurring velocities are impacted by rotation, strength of the wind stress τ, as well as depth of the basin h max .
We show simulation results on a uniform grid using 81,920 triangular linear elements with radii of the inscribed circles of about Δr = 38.5. This corresponds to a spatial distance of Δx = 781.25 and Δy = 156.25 between grid points. On the left most edge of the basin, we used transmissive boundary conditions. Impermeable boundary conditions are prescribed on all other boundaries. A quadratic Manning law was used with n = 0.0025 as was suggested in a previous study (Sanay and Valle-Levinson 2005) to parameterise the bottom friction and a β-approximation to model the Coriolis forcing as described in Sect. 2.1.1 was used with f 0 = 0.001.
The results are found in Fig. 10. For all six configurations the magnitudes of the velocities are plotted in colour as well as the velocity vectors for every 100th point. Note that configurations 5 and 6 are scaled by factors of 10 and 0.2 respectively to improve readability. Figure 11 shows cross sections of velocities (scaled with the maximum occurring velocity) at x = 50,000 (mid basin) and x = 98,000 (close to the head of the channel located at x = 100,000). In general, we observe a re-circulation zone of below 10,000 from the head which is in line with the observations reported in Fig. 9 Wind-induced circulation in a semi-enclosed homogeneous, rotating basin: sketch of initial condition. Depicted are the bathymetry and initial water surface (left) and the wind field (right)  (Sanay and Valle-Levinson 2005). The major characteristics of this test case are reproduced by the simulation results: We find a positive correlation between wind strength and velocities (two bottom panels of Fig. 10) as well as a negative correlation between water depth h and the magnitude of the velocities |u|. For the nonrotational case (configuration 1, top panel in Fig. 10), we agree with the observations from (Sanay and Valle-Levinson 2005) and find a symmetric profile of the meridional velocities. This changes for all rotational test cases in which we find asymmetries near the head (at x = 98,000-see Fig. 11). We furthermore find (see Fig.  11) steeper gradients on the left shoal for configurations 2 and 3. For configuration 4, Sanay and Valle-Levinson (2005) find a symmetric velocity profile-almost as for configuration 1-this, we cannot confirm. However, we would like to point out that they used different equations to model this problem. Hence, an explanation might be found, again, in the characteristics of the shallow water equations. They are depth-integrated, so that the Coriolis force effects the entire water column. The model employed in (Sanay and Valle-Levinson 2005) used between 10 and 30 vertical layers in their model which we believe to add a dissipation that we are unable to reproduce. Resolution appears to be a critical issue for this test case. With a spatial resolution of only Δr = 154 (a total of 5120 triangles in Ω), with Δr is the radius of the inscribed circle, we observe that after long integration times instabilities develop in the form of configurations (2-6), indicating that the resolution is not sufficient for a realistic and physically correct solution. We attribute this effect to the depth-integrated character of the shallow water equations as well as the occurrence of a geostrophic imbalance. As opposed to the model used in (Sanay and Valle-Levinson 2005) that included several vertical layers and with that the ability to dissipate energy in the vertical dimension resulting in a rotation of the fluid in the (y-z) plane, the depth integrated shallow water equations cannot take vertical motion into account.

Results II: Adaptive mesh refinement
In this section, we show and comment on the use of adaptive mesh refinement in the presented model. To do this, we will focus on two test cases: -Idealised storm surge modelling and sensitivity analysis (Sect. 4.1) -Idealised dam break (Sect. 4.2)

Idealised hurricane approaching a linearly sloping coast
In order to study the viability of the current model for use in hindcasts, we implemented an idealised test which is similar to a test presented in (Mandli 2011). It is designed to reproduce observations of a published Fig. 11 Wind-induced circulation in a semi-enclosed homogeneous, rotating basin: scaled velocities at t = 20,000 sensitivity study (Weisberg and Zheng 2006). The latter showed that hurricane flood intensity is sensitive to the storm's approach speed, direction of approach θ and landfall location. Our idealised test set up is defined as follows: Let Ω = [−200,000, 500,000] × [−300,000, 300,000] be a rectangular basin with a transmissive boundary at x = − 200,000, impermeable boundaries otherwise and a bathymetry defined by the piece-wise linear function & where α b = 0.025 is the slope of the bathymetry and x = (x, y) ⊤ is the spatial coordinate (see also Fig. 12). The initial water surface is at rest and described by h(x, 0) = max (3000.0 − b(x), 0.0). We then initialise a cyclone at an initial position (see large black dot in Fig. 12) and an approach angle θ. The corresponding wind stress τ requires the computation of a continuous wind field. This can either be accomplished by using reanalysis data or a parameterised model that allows the derivation of a continuous wind field from a few discrete parameters; see (Beisiegel and Dias 2019) for a short discussion on using a parameterised model in combination with re-analysis data for the Republic of Ireland. In the present model, the wind field is computed using a parameterised model (Holland 1980), which we elaborate on in more detail in Sect. 4.1.1, with parameters p c = 95,000 Pa, p n = 100,500 Pa, A = 23 and B = 1.5 which are representative for the 2008 hurricane Ike. Six different configurations as in (Mandli 2011) are implemented (cf. Table 3) and final times T end are chosen such that the storm's landfall is captured. Note, that we use oceanographic conventions, i.e. 0 ∘ corresponds to travelling to the right. The boundary conditions are transmissive at x = − 200 × 10 3 and impermeable otherwise. Transmissive boundaries for this subcritical flow were implemented following (Antonopoulos and Dougalis 2016) using a standard approach based on Riemann invariants.

Holland's model to compute hurricane winds
The wind stress in Eq. (8) depends quadratically on the wind v. For hurricanes, winds can be computed using the model (Holland 1980): where r is the distance from the centre of the storm; t the tangent to the circle with radius r, A, B ∈ ℝ are shape parameters; p n and p c are the ambient and central pressure respectively; ρ the water density; and f the Coriolis parameter. The scaling parameters A and B are then obtained from the maximum wind speed as well as the radius of maximum winds (RMW): where e is Euler's number. An example for a normalised wind profile is found in Fig. 2. The wind model (Holland 1980) also gives a corresponding atmospheric pressure as is seen from Fig. 2 (right) which shows pressures (dashed line) normalised by p n − p c as well as wind speeds normalised by max |v| 2 .

Numerical results
We ran the simulation with a temporally changing time step, an explicit RK23 time integrator, and a CFL number of 0.2 with a spatial resolution of Δx = 10,937.5 a Manning parameter of n = 0.001, and the wind drag as in (Weisberg and Zheng 2006). We note that this CFL number differs from the one used Fig. 13 Idealised hurricane approaching a linearly sloping coast: waterfall plot of time series at even numbered wave gauges from bottom to top for all six configurations. Amplitudes are in m with an added offset of 10 · k for wave gauge G k for all k Table 4 Idealised hurricane approaching a linearly sloping coast: maximum wave height η max with respect to ramp up time at gauges G 10 , G 12 and G 20 Ramp up time 15 min 30 min 1 hour 2 hours 4 hours η max at G 10 41.4 × 10 −2 41.4 × 10 −2 41.7 × 10 −2 41.6 × 10 −2 41.7 × 10 −2 η max at G 12 6.5 × 10 −2 6.4 × 10 −2 6.5 × 10 −2 6.3 × 10 −2 6.5 × 10 −2 η max at G 20 2.3 × 10 −2 2.2 × 10 −2 2.2 × 10 −2 2.3 × 10 −2 2.4 × 10 −2 in all other experiments in this study but was found to be necessary to ensure numerical stability. We then compared the wave signal η(x, t) = h(x, t) − 3000.0 at numerical wave gauges G k located at G k · 10 −5 = (4.5, −2.5 + 0.5k) ⊤ for k = {0, 2, 4, …, 20} (see left display of Fig. 12) with the findings in (Mandli 2011) and found overall good agreement. The results are found in Fig. 13. It shows the water wave signal for all six configurations obtained at the numerical wave gauges G k for all even numbered wave gauges with the amplitudes plotted in metre with a vertical offset of 10 · k for gauge G k to increase readability. In agreement with the literature (Weisberg and Zheng 2006), we find that the observed flooding is sensitive to the approaching angle. The plots for configurations 2 and 3 in Fig. 13 show significantly different signals at the wave gauges to the left and right of the wave gauge at which the storm made landfall. They differ in shape and arrival time. The general N shape of the largest waves as seen in (Mandli 2011), however, could not be reproduced. Since higher resolution simulations with a halved Δx showed the same behaviour, we attribute this effect to the implementation of impermeable boundary conditions in this test case. The rotational direction of the wind velocity, and the impermeable boundaries at x = ± 300 × 10 3 are thought to be responsible for the only approximate symmetries between configurations 2 and 3, as well as 4 and 5. In the absence of Coriolis forcing, the rotational direction impacts the flow, such that the plot of configuration 2 in Fig. 13 is not exactly the same as the plot for configuration 3 mirrored at y = 0. We furthermore remark that all configurations show small oscillations right after the storm made landfall. This is due to the treatment of the wind stress which was switched off after the storm reached the beach.

Influence of ramp up times and robustness
Numerical models often require a gentle ramping up of source terms in order to reduce spurious oscillations and to allow for a robust computation. In our simulations, we ramped up the wind stress τ and pressure p using an exponential blending in time. For early times, t ≤ t ru this filter F takes on the form with ϕ the quantity that is to be started and c f a tunable coefficient. The storm only starts travelling towards the coast with angle θ after the ramp up time t ru is reached. Before, it is kept at its starting position, so that the wind and pressure fields are slowly ramped up until they reach their full strength. We ran configuration 1 with five different ramp up times between 15 min and 4 hours (see also Table 4) to test the robustness of the results. Ideally, we would like t ru to be as small as possible to save computing time but large enough to not pollute numerical results. We observe that for all times between 15 min and 4 hours we get robust −10.0 × 10 −2 −9.9 × 10 −2 −9.8 × 10 −2 −9.5 × 10 −2 12 η max 6.3 × 10 −2 6.3 × 10 −2 6.5 × 10 −2 6.5 × 10 −2 η min −10.5 × 10 −2 −10.5 × 10 −2 −10.4 × 10 −2 −10.0 × 10 −2 20 η max 2.3 × 10 −2 2.3 × 10 −2 2.4 × 10 −2 2.4 × 10 −2 η min −8.1 × 10 −2 −8.0 × 10 −2 −7.9 × 10 −2 −7.6 × 10 −2 Fig. 14 Idealised hurricane approaching a linearly sloping coast: maximum wave height η max versus maximum wind speed, max v, for different drag coefficients (black squares-Weisberg and Zheng, red plusses-Powell, blue xs-Garrat and grey circles-constant. Because of the close agreement of obtained η max , markers overlay numerical results. Furthermore, as is seen from Table 4, the maximum wave height does not show a lot of variation depending on different ramp up times and the maximum variation is found to be 3 × 10 −3 . Wave gauge signals, however, detect a small wave at times t ≈ 1 h for ramp up times below 2 hours. Hence, we chose a ramp up time of t ru of 2 hours for this study. in order to determine their influence on maximum wave heights at wave gauges G k which we assume a good indicator for wave run-up at the coast. Table 5 shows maximum and minimum wave heights for the original configuration 1 at selected wave gauges closest to the storm's landfall. We repeated the simulation with all four different wind drag models and observe merely minor differences of the order of at most 0.5 × 10 −2 . The parameter c d = c d (v) depends on the wind speed. Therefore, we ran configuration 1 with varying maximum wind speeds max | v | ∈ {15, 25, 35, 45, 55} for every wind drag model. As shown in Fig. 14, we see that the differences in maximum wave heights η max are negligible. In fact, they are, again, of the order of at most 0.5 × 10 −2 . Hence, we conclude that in an idealised model such as the one presented in this manuscript different wind drag models do not lead to significantly different results. This can be explained by the form of the wind stress τ = ρ a c d ∥ v∥ 2 v. Using the selected parameterisations, c d will differ at most by a factor of 2 in very localised regions of the storm, which does not lead to a significant increase or decrease in the observed wave heights close to the coast.

Influence of storm size
The size of a hurricane plays an important role in the observed flooding. For reasonable storm sizes, a variability of about 30% in observed surge is reported (Irish et al. 2008). In the wind model (Holland 1980), the shape and size of the storm depend on the shape parameters A and B. These, in turn, depend on the radius of maximum winds (RMW), the difference between ambient and central pressure Δp = p n − p c , the air density ρ a and the maximum wind speeds as shown in Eq. (11). Assuming storm conditions that are representative for the 2017 hurricane Ophelia, we simulate configuration 1 as described above and vary the radius of maximum winds. According to (Hsu and Yana 1998), the average radius of maximum winds of hurricanes is 47 km, with a standard deviation of 27 km which is why we chose to run simulations with the radii stated in the Table 6. As can be seen, only parameter A varies with varying RMW if all other conditions are kept fixed as it describes the radial scaling on the RMW and the location of the maximum wind relative to the origin. We measured the maximum wave height at wave gauge G 10 -the location at which the synthetic storm made landfall-and will hence record the maximum surface elevation η max . The results are depicted in Table 7. We see that with this simple parameterisation, we achieve measured maximum wave heights with a variability of about 39%. Given that we tested with parameters resembling hurricane Ophelia for the most part, we conclude that we are within the range of variability that was found in (Irish et al. 2008).

Adaptive simulations
Dynamically changing non-uniform meshes as described in Sect. 2.2 are ideal for simulating localised phenomena at a reduced computational cost. Since storm wave heights are strongly influenced by changes in bathymetry as well as the size and strength of a storm, we define a refinement indicator η Ω i to take both of them into account for every element Ω i : Using the indicator (Eq. (12)), we achieve a refinement of the beach or bathymetry gradient as well as the storm position (see Figs. 15,16,17,18,19 and 20). Figures 15, 16, 17, 18, 19 and 20 show numerical results for configurations 1-6 on an adaptive and a uniform mesh. Plotted are the non-uniform meshes and currents at the initial time, about halfway to landfall and close to landfall. We can see that the adaptive and uniform simulations yield comparable results with small errors in areas of high resolution and slightly larger errors outside. The refined region in the adaptive simulation comprises the entire area close to the coast as well as the cyclone and with time we see that the refined area resolving the cyclone moves with the storm. Some numerical error is observed which we attribute to the choice of the heuristic refinement indicator which captures physical features that only indirectly correlate with numerical error and model sensitivities. This, Fig. 21 Idealised hurricane approaching a linearly sloping beach: Waterfall plot of wave gauges G k over time for uniform (black dashed line) and adaptive (blue solid line). Amplitudes are in m with an offset of 10 · k m for wave gauge G k Fig. 22 Idealised dam break: Domain description with inlet depicted by a diamond however, does not impact the measured wave heights near the coast. The maximum wave heights measured at gauges G k are depicted in Fig. 21 and show good agreement between the adaptive and uniform simulations. Hence, since our interest lies in the wave gauge signal, the refinement indicator seems suitable. Although dynamic mesh refinement adds computational overhead, this was not found to be significant. For configuration 1 (Fig. 15), the uniform simulation comprised 8192 elements, while the adaptive simulation only used on average: 2714 (and a maximum 3074) elements-a reduction to at least 37.52% of the elements (on average 33.13%). In terms of run time, the adaptive simulation used about 31.33% of the computational time by yielding quantitatively the same result in terms of measured wave heights. A similar behaviour can be observed for all configurations. with h 1 = 0.6 and h 2 = 0.05. We furthermore assume zero initial velocities u(x, 0) = 0. At time t = 0, a dam break is simulated by removing the wall between 3.95 ≤ y ≤ 4.35 as indicated by the diamond in Fig. 22.

Idealised dam break
To preserve the space-filling curve (SFC) ordering of the elements and to resolve the narrow inlet, we used the macrotriangulation depicted in Fig. 23 as an initial mesh for the uniform as well as adaptive mesh refinement. This comes at the expense of the time step being limited by the narrowest element which we will further comment on later in this paragraph.
We ran the uniform simulation with a CFL number of 0.3 which resulted in a time step of about Δt = 5 · 10 −4 and a total number of elements of 69,632 which corresponds to a 12 times uniform refinement of the mesh depicted in Fig.  23. Snapshots of the numerical solution at times t = 0, 2, and 4 on a uniform and adaptively refined mesh are found in Fig. 24. We observe that a large wave develops from the inlet and starts travelling across the shallow part of the domain. Using the fine resolution uniform simulation as a reference solution, we ran an adaptive simulation, refining according to momentum maxima, i.e. using a refinement indicator η Ω i t ð Þ ¼ max x∈Ω i ju x; t ð Þj. The adaptive mesh comprised a finest resolution identical to the uniform simulation but a coarsest resolution resulting from only 8 times refining the macro-triangulation. The results on the adaptive and uniform mesh are virtually indistinguishable. We can see that the fine mesh area follows the emerging wave and only regions are refined that experience fluid movement. Due to the not uniform macro-triangulation, the region around the inlet is highly refined, ensuring that flows are accurately captured. The adaptive mesh comprised on average 8437 elements which is about 12% of the number of elements of the uniform simulation. In terms of computing time, the adaptive simulation only took 30.5% of the computing time for the uniform simulation-leading to a cost reduction of almost 70%.

Conclusions and future work
In this study, we have developed a discontinuous Galerkin model on a dynamically adaptive triangular mesh that solves the fully 2D non-linear shallow water equations for the simulation of coastal flooding and idealised storm surges. Numerical test cases that we believe to be a good basis for a test suite that might be useful for storm surge modellers demonstrate that the obtained model is inundation stable due to advanced slope limiting techniques (Vater et al. 2019) and maintains important conservation properties such as mass as well as integrated fluxes as described in Sect. 3.1 for the simulation of a hypothetical embayment. Moreover, a steady state is achieved numerically in Sect. 3.2 in which a balance between pressure gradients and wind stress is simulated. In Sect. 3.3, we show the robustness of the wind forcing and the effect of Coriolis forcing on wind-induced circulation. We furthermore see that the wind forcing is robust with respect to ramp up times; hence, no spurious artefacts are introduced.
Finally, in Sect. 4.1, we show the capability of the model to simulate idealised hurricane storm surge using the wind parameterisation (Holland 1980), hence making it suitable for simulation of realistic hurricanes such as the 2008 Atlantic hurricane Ike or the 2017 hurricane Ophelia. A sensitivity analysis furthermore reveals that the model is not sensitive to the choice of wind drag parameterisation or storm ramp up time. The observed variability of maximum wave heights (and therewith wave run up) with varying RMW confirms previously published studies, underlining the capability of the model to yield realistic results. Most notably, using dynamically adaptive meshes, we obtain virtually the same signal at wave gauges close to the beach at significantly less computational cost: The reduction of computing time was measured to be up to 70%. Overall, this is to demonstrate that the developed model is suitable for the simulation of idealised hurricane storm surge and shows a satisfactory robustness and accuracy as well as adaptive mesh capabilities that help reduce computing costs significantly.
The same reduction of about 70% could also be shown for the idealised dam break problem in which we showed that the mesh is accurately following the emerging waves.
In this study, we have dealt with a number of idealised test cases and demonstrated the model's potential to use an adaptive mesh for the simulation of hurricane storm surge. The application of the presented model to more realistic data is beyond the scope of this paper and will be left for future research. The results presented, however, allow the conclusion that the combination of dynamically adaptive mesh refinement with a DG discretisation significantly increase the potential practicality of the model and can be seen as a first indication of DG methods being a useful tool for the application to hurricane storm surge modelling.