Coastal Flood Modeling Challenges in Defended Urban Backshores

: Coastal ﬂooding is a signiﬁcant and increasing hazard. There are multiple drivers including rising coastal water levels, more intense hydrologic inputs, shoaling groundwater and urbanization. Accurate coastal ﬂood event prediction poses numerous challenges: representing boundary conditions, depicting terrain and hydraulic infrastructure, integrating spatially and temporally variable overtopping ﬂows, routing overland ﬂows and incorporating hydrologic signals. Tremendous advances in geospatial data quality, numerical modeling and overtopping estimation have signiﬁcantly improved ﬂood prediction; however, risk assessments do not typically consider the co-occurrence of multiple ﬂooding pathways. Compound ﬂooding refers to the combined effects of marine and hydrologic processes. Alternatively, multiple ﬂooding source–receptor pathways (e.g., groundwater–surface water, overtopping–overﬂow, surface–sewer ﬂow) may simultaneously amplify coastal hazard and vulnerability. Currently, there is no integrated framework considering compound and multi-pathway ﬂooding processes in a uniﬁed approach. State-of-the-art urban coastal ﬂood modeling methods and research directions critical to developing an integrated framework for explicitly resolving multiple ﬂooding pathways are presented.


Introduction
Urban coastal flooding is a global humanitarian and socioeconomic hazard with multiple drivers including rising coastal water levels, more intense hydrologic inputs, shoaling groundwater and increasing urbanization. Currently, over 20 million people reside below present day high tide levels, 200 million are vulnerable to storm flooding [1], and over 600 million individuals reside in the coastal zone [2]. Mean sea levels are expected to rise 0.28-0.98 m by 2100 [3]. These estimates are, however, likely under-representative of potential sea level rise rates [4,5]. Regional trends show significant variability [6]. Relatively modest sea level rise (i.e., 0.50 m) will significantly increase flood frequencies [7]. Sweet and Park [8] showed that "tipping points", i.e., flooding over 30 days per year, will be reached by 2050 and flood frequency will increase drastically (e.g., near daily flooding under RCP4.5 scenario) by 2100 for many locations.
Urban flood events are the most significant contributor to the overall flood risk [9,10]. Hanson et al. [11] suggested a threefold increase in coastal population exposure by the 2070s which will be exacerbated as low-lying areas are urbanized. From an economic perspective, assets exposed to the 100-year coastal flood are expected to increase dramatically in the coming decades. Hallegatte et al. [12] suggested a near order of magnitude increase in flood damages by 2050.

State of Knowledge
Scientific literature presents a general modeling methodology (e.g., [20,22,23,[27][28][29][30][31][32][33][34][35]); however, there has been a call for more rigorous regional flood modeling and improved methodologies (e.g., [36][37][38]). Urban flooding is the most significant contributor to flood risk [9,10], and poses a number of modeling challenges, including offshore boundary conditions, terrain and infrastructure depiction, representing hydrologic inputs, characterizing spatially and temporally variable wave runup and overtopping flows, inclusion of hydraulic infrastructure, routing overland flow and quantitative model evaluation. Overland flood prediction is presented as central to the modeling framework, with offshore and hydrological forcing considered as external boundary conditions ( Figure 1). Ideally, computationally efficient couplings that explicitly model all relevant flooding pathways will be developed.

Overland Flow Modeling
Overland flow may be routed via simplistic planar surface projections [37,39], mass conservation schemes (e.g., LISFLOOD-FP [27]), simple inertial formulations of the shallow water equations (e.g., [40]) or through more sophisticated hydrodynamic models (e.g., AdH, MIKE FLOOD, TUFLOW, DIVAST, BreZo, Delft3D, and Telemac) which solve the full 2D shallow water equations. The static method, also known as planar surface projection, equilibrium, or "bathtub" method, has drawn substantial criticism for poor predictive skill in protected urban backshores [20,41]. Poulter and Halpin [29] and Heberger et al. [37] presented a raster-based modeling approach where areas that fall below the water level are flooded, i.e., the assumption is that hydraulic flow paths exist and the flood is sustained sufficiently long to fill the impacted region up to the height of the embayment. Variants of this approach can be devised to account for protection by levees and seawalls, but all static models retain the fundamental assumption that flooding occurs instantaneously upon exceeding overflow thresholds. Static methods are an "all or nothing" approach. If the water level exceeds the threshold, all connected areas are predicted to be flooded; alternatively, if the water level is below the flooding threshold, all areas are presumed to be dry. This leads to extraordinary under-and over-prediction [20,21]. Since flow is not hydraulically routed, depth and velocity information are unavailable with static models. Figure 2A shows a static prediction and Figure 2B shows a hydrodynamic simulation along with the observed flood extent. Static models perform poorly and should be avoided for low-lying defended urban backshores.
An alternative method is to extend the hydraulic model domain to all areas subject to episodic inundation and to make flood mapping an integral part of the hydraulic model analysis. For example, Brown et al. [28] presented a coupled storm surge and overland flow model for Canvey Island located in the Thames Estuary. There are two primary benefits to this approach: hydrodynamic information is explicitly calculated (e.g., depth, velocity, arrival time) and predictions are substantially more accurate in urban backshores.¯Ā

Shallow Water Equations
Overland flow has been investigated using kinematic and diffusive wave models that do not consider inertial terms in the momentum balance. However, in urban environments where flow transitions are common and high resolution depth and velocity estimates are desirable, inertial terms improve predictions [42]. Hunter et al. [43,44] showed that diffusion wave models may be less computationally efficient than solving the full 2D nonlinear shallow water equations because the optimal time step is a quadratic function of the grid size. The efficiency of diffusive models is thus lost for fine-scaled grids required to resolve hydraulic pathways in urban areas. Bates et al. [40] proposed an "inertial" model removing the convective but retaining the local acceleration term which improved computational times and reasonably simulated depth and velocities. Neal et al. [45] showed that such an approach leads to reduced computational time compared to full shallow water models for gradually sloping, subcritical flows. Inertial models and their parallelized versions (e.g., [46]) thus offer a potential option for probabilistic and large-scale model applications which require computational efficiency [24]. These simplified models can accurately predict flood extent in slowly temporal-spatial varying applications, but momentum terms are critical to accurate velocity estimates [42,45].
Full 2D nonlinear shallow water models have proven to provide excellent depth and velocity estimates in urban areas (e.g., [47,48]). Particular attention should be paid to the numerical implementation and its suitability for rapidly variable flows. Discontinuous Galerkin finite element methods or Godunov-based finite volume schemes admit supercritical flows from abrupt elevation changes inherent to urban environments such as sea walls, streets, and curbs without case specific parameter tuning. Finite volume schemes are the most widely used 2D nonlinear shallow water solution method [49]. Numerous Godunov-type finite volume codes have been successfully implemented in coastal embayment modeling [20,[50][51][52] and urban flood simulations [47,48,[52][53][54][55].

Terrain and Infrastructure
Accurate topographic representation is fundamental to predicting coastal flooding. Sanders [56] showed high resolution LiDAR-derived digital elevation products significantly outperformed larger scale mapping products. For example, digital elevation models (DEMs) derived from airborne interferometric synthetic aperture radar (IfSAR) and shuttle radar topography mission (SRTM) data are surface products with significant noise and do not sufficiently resolve hydraulic topography, which negatively impacts flood predictions [56].
Multiple studies have noted that low relief areas are especially sensitive to terrain representation [57,58]. Generally, LiDAR-derived data with a spatial resolution of ∼1 m postings and vertical accuracy ∼10-15 cm are adequate for terrain resolution [48,56,58,59]. However, Néelz et al. [60] found significant LiDAR limitations for resolving built areas, while Webster et al. [61] and Gallien et al. [21] suggested LiDAR inadequately resolves abrupt elevation changes such as wharves or sea walls. Similarly, Gallien [34] showed highly-filtered bare earth DEM products may misrepresent topography, particularly between buildings or other tall structures. Numerous studies point to the importance of resolving key hydraulic features such as roads, walls, storm drainage and dunes which restrict or redistribute flow (e.g., [9,20,23,29,41,47,48]). Gallegos et al. [48] showed that street depressions must be reasonably resolved for hydraulic routing and suggested mesh resolutions of ∼5 m or three cells balances accuracy and computational effort. Gallien et al. [20] showed that, for accurate results of weir-like sea-wall overflow, wall elevations must be accurately surveyed and explicitly resolved in the hydraulic mesh ( Figure 3). These studies emphasize the need for high accuracy real time kinematic (RTK) or post-processed GPS surveys of hydraulically important urban topography, but only a limited number of studies have attempted to explicitly include infrastructure in hydrodynamic coastal flood modeling efforts (e.g., [20,21,23,28,34,[62][63][64]). Data preparation, whether the product is a digital surface model (DSM) including canopy, buildings, vehicles and other surface features or is a digital terrain model (DTM) representing bare earth topography is of key interest. Urban flood models fall in one of the following four categories depending on the method of representing buildings in the terrain [53]: building-hole models, where the buildings are replaced with building-perimeter polyline [47], and either a free-slip boundary [54] or other internal hydraulic conditions are imposed along the perimeter; building-resistance models, where a relatively large resistance parameter (bottom friction) is assigned to cells corresponding to building footprints [65] or developed parcels; building-porosity models, where buildings are represented through spatially-distributed porosity and drag coefficients, without resolving the exact geometry of the buildings [52,66]; and building-block models, where buildings with their representative shape and height are included in the terrain, essentially blocking the flow from passing through them (e.g., [67]). Each method requires a different mesh style. All methods reasonably predict flood hydrographs and extent when the model resolves important flow paths and storage functions of the land surface [53], but velocities depend on building representations [53,64]. The model of choice depends on availability of computational power, the importance of velocity predictions, and the availability of the necessary data coupled with the user's familiarity with the GIS techniques needed for each method implementation [53].

A B
In urban applications, significant portions of the flow may be conveyed by storm sewers. Henonin et al. [68] presented an urban flood model classification that represents the storm system and overland flow model dimensionality (e.g., 1D-1D, 1D-2D, etc). SWMM [69] is a widely used link-node approach for resolving subsurface flows for 1D-2D modeling [42]. In coastal applications, storm system drainage is critical to predicting flood extent for three reasons: (1) bay or ocean water levels may flood low protected backshores through the storm drain outfalls; (2) the storm drain system may redistribute water in a domain; and (3) overland flow volumes may be reduced by storm system drainage. Gallien et al. [20] showed that the storm drain system redistributes flood water even when the outlets are closed to prevent back-flooding from high bay levels. That is, water entering one curb inlet may cause flooding through subsurface hydraulics connections even if the surface topography is not hydraulically connected. This interaction cannot be resolved in models ignoring drainage, which challenges the traditional assumption that storm system flows can be ignored when the drainage system is operating at capacity [70,71].
Infrastructure failure should be explicitly incorporated into flood predictions e.g., [28,72]. However, generally the approach is "all or nothing". Infrastructure is either omitted from coastal flooding models (nothing) resulting in substantial vulnerability overestimation, or infrastructure is resolved (all), but potential failures are not considered, resulting in vulnerability underestimation. The uncertainty generated by failure of protective barriers such as levees can be examined by integrating scenarios of flood defense breaches into both static and hydrodynamic methods (e.g., [28]), however only the hydrodynamic approach accounts for storage and resistance effects, and provides detailed information about the velocity and depth distribution. In the case of dynamic protective barriers (e.g., sand dunes and berms), sediment transport models can be applied to explicitly account for dynamic breaching during a storm event. For example, Elsayed and Oumeraci [73] calibrated XBeach, applied it to predict a dam breach and subsequent overland flow, and suggested that flow and velocity were in good agreement with observed values.
"Soft" infrastructure such as beach and dune restoration and artificial dune building are often the first line defense in coastal management, but little is known about their performance and design. Often beaches are only crudely resolved from a single survey, or estimated from LiDAR data, which may under-or overestimate vulnerability. Similarly, morphological evolution during energetic wave events is often unresolved, which significantly alters overtopping flood volume estimates [74]. Artificial dunes (also referred to as "berms") provide critical backshore protection and are widely deployed as a first line of defense for large wave and/or high water level events [72,75]. These structures are often sacrificial, intended only to deflect specific high water or energetic wave events, and are often unresolved in individual LiDAR surveys. Temporary berming is a popular mitigation strategy for urban coastal flooding (e.g., [72,[76][77][78][79][80]). Gallien et al. [21] considered flood mitigation effects of artificial dunes and showed that the backshore flooding extent increased nearly 400% without a temporary dune ( Figure 4).

Boundary Conditions
Urban coastal flood models require characterizing both open coast and estuary water levels, hydrologic inputs and wave overtopping flows. Sea level rise will increase oceanic and bay water levels and allow wave energy to propagate landward, potentially increasing wave overtopping flooding.

Water Level
Water level boundary conditions may be represented by simple arbitrary values, output from multi-scale models, probabilistic or synthetic data. Sea level rise assessments commonly utilize a single water level value depicting a future change in sea level (e.g., [59,[81][82][83][84]). Alternatively, probabilistic water levels may be used to force an overland flow model [31]. Kopp et al. [85] determined probabilistic water levels for global tide gauges through the 22nd century. Sophisticated nested modeling systems (e.g., [33,35]) downscaled global climate models to consider regional water levels to predict regional water level impacts.
In open coastal areas, ocean water levels are sufficient for flood model boundary conditions, however, estuary and embayment water levels may not correspond to open coast water levels and require explicit embayment resolution [32]. For example, the southern lobe of San Francisco Bay experiences tidal amplification (∼1.4×) [86], and Lyddon et al. [87] showed significant spatial variability of water levels from the combined effects of tide and storm surge in hyper-tidal estuaries. Gallien et al. [20] suggested that even relatively small (∼3 cm) water surface elevations in bays can cause significant changes in protected urban flood prediction. More sophisticated hydrodynamic event simulation at the regional scale (<100 km) has been approached by establishing a simulation domain wherein hydraulic models are applied to simulate spatial and temporal changes in bay water levels in response to boundary forcing (e.g., [20,28,[33][34][35]). Critically, hydrologic inputs may superelevate bay and estuary water levels.

Hydrologic Impacts
Hydrologic impacts on coastal flooding may be through fluvial (river) or pluvial (runoff) pathways. From a fluvial perspective, previous research suggested river hydrographs peak subsequent to storm surge and may not affect coastal forcing (e.g., [16,17]). However, Svensson and Jones [88] suggested that steep catchments are particularly vulnerable to compound flooding because maximum surge and riverine discharges may be coincident. Similarly, Chen and Liu [89] showed the importance of storm surge and river discharge phasing, while Orton et al. [90] and Maskell et al. [91] showed that exclusion of freshwater inputs can lead to significantly lower projected bay and estuary water levels predictions. Orton et al. [90] noted that omitting freshwater inputs and density variations led to low bias, which is particularly concerning for urbanized areas with extensive subsurface infrastructure. Ward et al. [92] mapped sea level and fluvial input dependence and showed that over half the global sites exhibited significant dependence. Compound fluvial flooding is often approached on the local scale using deterministic hydrodynamic modeling (e.g., [89,90,93,94]), while larger scales typically consider probabilistic methods (e.g., [15,19]). For example, Olbert et al. [95] and Gallien et al. [62] used deterministic modeling to show that sheltered urbanized estuarine areas are particularly sensitive to streamflow. Moftakhari et al. [19] used a probabilistic approach to show that neglecting compound impacts significantly underestimates flooding hazards in a number of large urban areas and advocated for local-scale research to accurately quantify impacts. Numerous research groups have pointed to the importance of resolving compound fluvial coastal flows for accurate coastal flood predictions (e.g., [19,62,90,91,[95][96][97]). Critically, a paucity of observational data limits flood risk assessment, and the collection of high quality observations is advocated [93,98].
From a pluvial perspective, storm surge may saturate the area and rainfall becomes runoff and potentially, storm system flow. During high coastal water level events, storm systems are often closed to prevent back flooding, which exacerbates pluvial flooding. Although rainfall and drainage have recently been recognized as a significant issue [99,100], only limited efforts incorporate this key hazard. Cheng et al. [101] developed the COSM model, which couples nearshore (ADCIRC) and watershed (pWASH123D) processes, and suggested a tight two way coupling is necessary to accurately resolve nearshore-watershed interactions. Tang et al. [102] presented a hybrid modeling approach, where a 3D ocean model (FVCOM) is coupled to a 2D nonlinear shallow water model to explicitly resolve marine flooding, which is combined with spatial output from a probabilistic hydrologic model to resolve compound flooding from various storm events. Thompson and Frazier [103] estimated potential future flooding hazards using a 2D hurricane surge model (SLOSH) and inland precipitation using the Interconnected Channel and Pond Routing Model (IPCR) model, and showed that precipitation significantly increases the flood extent. Joyce et al. [94] used an ocean (ADCIRC) and a nearshore wave transformation model (SWAN) in a one way coupling to force a hydrological/hydraulic and storm water model (ICPR). With the exception of Cheng et al. [101], these analyses did not consider the subsurface groundwater flows.
From a flooding perspective, groundwater may inundate coastal areas [106,107]. High beach groundwater levels may reduce infiltration processes which have been associated with increased wave runup [113,114]. Areas with shallow aquifers are particularly vulnerable to emergence [107]. Underground infrastructure (e.g., pipes, tunnels, and basements) will become increasingly vulnerable to rising groundwater [105]. Coastal groundwater shoaling and emergence generally occur more frequently as the vertical unsaturated subsurface space narrows [115]. The inclusion of groundwater inundation in modeling has been shown to result in more than twice as much coastal flooding compared to marine inundation alone [106].

Wave Overtopping
Backshores historically unimpacted by incident wave overtopping, may become vulnerable with higher sea levels, and Arns et al. [116] suggested that future water levels will amplify wave runup design heights by ∼50% in currently impacted areas. Although wave overtopping has received sustained attention in the literature (e.g., [117][118][119][120]), wave runup and overtopping are considered significant deficiencies in current modeling methodologies [28,121,122], and have been recognized as a key future research area [123].
Wave forcing may be represented by known deep water wave characteristics (i.e., buoy data) or more sophisticated numerical wave modeling efforts. Flood vulnerability assessments often use simple empirical runup models (e.g., [124]) driven by deep water wave data. This method cannot capture nearshore bathymetry effects or key nearshore processes (e.g., wave-current interaction). An alternative is to propagate waves into the region of interest using SWAN [125] or other phase-averaged (spectral) wave models. From the spectral parameters, a random-phased wave field can be generated in a phase-resolving model (e.g., [126,127]). Stochastic uncertainty of the random phasing should be addressed when computing runup/overtopping metrics [128]. In areas of high urbanization with reflective armoured shores, significant inter-tidal or submerged infrastructure (e.g., harbors, jetties, and breakwalls), Boussinesq-type [129] or other non-hydrostatic phase-resolving models (e.g., [130]) have proven to adequately represent the nearshore wave field.
A simplistic method for depicting flooding from overtopping involves adding maximum wave runup to determine a total water level (e.g., [37,131]) and projecting this water level across the land surface. The total water level method is applied using a bathtub model and consequently suffers identical deficiencies shown in overland flow modeling. Overtopping time scales in episodic flooding events caused by a coincident large wave conditions and high tides range from minutes to a few hours, insufficient time to fill the backshore, and consequently total water level methods for wave overtopping have proven to significantly overpredict the inundation zone in defended urban backshores [62]. Overtopping is inherently a dynamic, temporally variable process which can be represented using empirical and numerical models.
Only a limited number of studies have attempted to include wave overtopping volumes using empirical estimates or numerical estimates [31,[132][133][134][135]. Numerical models represent the current state-of-the-art for simulating overtopping flows and theoretically, if the physics are well represented, could predict overtopping in an infinite number of dune, dike or wall configurations. However, field scale implementations have been challenged by computational effort, grid sensitivities, and boundary condition effects which have restricted most applications to numerical wave flumes (e.g., [136]) or analytical solutions and laboratory validation data (e.g., [121,[137][138][139]).

Empirical Wave Overtopping Models
Literature presents a number of empirical wave overtopping models (e.g., [140][141][142][143][144][145][146][147][148][149][150][151][152][153]). Typically models calculate a dimensionless overtopping rate that considers freeboard (crest elevation above still water level) and wave characteristics (significant wave height and peak period). EurOtop has been widely used in estimating overtopping volumes for coastal flooding applications (e.g., [21,34,123,132,154,155]). The general EurOtop average overtopping formula for coastal dikes and embankment seawalls is, q where H m0 is the significant wave height at the toe of the structure, α is the slope, g represents gravity, q is the mean overtopping rate per unit length, γ b is the berm influence factor, γ f is the roughness influence factor, γ β is the oblique wave attack factor, and γ v is the vertical wall influence factor [153].
The van der Meer [149] formulation relies on a breaker parameter ξ m−1,0 , which characterizes the wave breaking condition (i.e., breaking, non-breaking) and is given as, where L m−1,0 is a spectral deep-water wave length. For a complete discussion of the TAW and EurOtop equations, refer to the works of van der Meer [149], Pullen et al. [151] and van der Meer et al. [153]. It should be noted that EurOtop and other empirical wave overtopping formulae provide average overtopping rate estimates and do not resolve individual wave overtopping volumes or infragravity energy which are important to wave overtopping flooding.

Numerical Wave Overtopping Models
Numerical modeling is attractive for resolving overtopping flows because, theoretically, a validated model may be configured for any desired structure, bathymetry or wave train. Additionally, numerical models can simulate impulsive, temporally variable overtopping rates and facilitate dynamic (wave-by-wave) coupling with the overland flow model. A variety of methods to resolve wave overtopping have been investigated; the Navier-Stokes equations, Boussinesq-type models and the Nonlinear Shallow Water (NLSW) equations.

Navier-Stokes Equations
The full incompressible Navier-Stokes equations provide an excellent model for wave overtopping and can resolve complicated hydrodynamics such as wave breaking and overcome limitations associated with any particular wave theory [137]. Significant numerical modeling work resolving wave runup and breaking processes has been accomplished (e.g., [156][157][158][159][160][161][162]). A selection of these models were subsequently investigated as overtopping models. In the Eulerian framework, Liu et al. [138] presented a 2D Reynolds Averaged Navier-Stokes model based upon the Lin and Liu [159] RANS solver with a k − turbulence closure model. The model was named Cornell Breaking Waves and Structures (COBRAS) and was applied to regular wave overtopping of a porous structure. Comparison to laboratory data showed COBRAS correctly simulated the free surface and overtopping rates [138]. Losada et al. [137] employed an extended COBRAS model, COBRAS-UC, developed by the University of Cantabria (e.g., [163][164][165][166]), which uses a volume averaged Navier-Stokes turbulence model developed by Hsu et al. [158], and successfully simulated both regular and irregular wave overtopping of a rubble mound breakwater. Similarly, Soliman [167] and Reeve et al. [168] employed the Lin and Liu [159] RANS formulation to investigate combined overtopping and overflow of smooth impermeable sea walls with small and negative freeboards. Results generally agree with previous models; however, for combined overflow and wave simulations, significant overprediction is observed. Thompson et al. [169] used a RANS formulation to predict seawall overtopping and showed generally good agreement with EurOtop average overtopping rates. A three-dimensional (3D) LES model was proposed by Li et al. [170] to study overtopping of sloped and vertical structures and produced promising results when ran on high resolution grids; however, grid dependency was observed. Similarly, Okayasu et al. [139] extended a 3D LES model to investigate wave overtopping on gentle slopes and stepped sea walls, and results suggested a sensitivity to water depth at the structure toe. Lagrangian methods for solving the Navier-Stokes have gained some attention within the scientific community and are attractive because the method is meshless and requires no special surface tracking [157]. Overtopping has been considered using SPH in small scales (e.g., [171][172][173]) with generally good results. Although both the Lagrangian and Eulerian methods of solving the Navier-Stokes equations have proven successful for modeling wave overtopping volumes, computational effort is extremely high and currently impractical for field scale or long wave train (i.e., 1000 waves) simulations.

Boussinesq
The classical Boussinesq equations for variable depths were developed by Peregrine [174] and benefit from the ability to include frequency dispersion; however, these equations were valid only where depth to deep water wave length ratio is less than 0.2 [175]. Subsequently, considerable effort to extend the equations applicability into deeper water has occurred (e.g., [176][177][178][179][180]) and numerous studies have extended into the surf and swash zones (e.g., [129,[181][182][183][184][185]). Stansby [186], Stansby and Feng [187] and Orszaghova et al. [188] considered solitary wave overtopping, whereas Lynett et al. [134] considered 1D wave spectra. Sitanggang and Lynett [189] proposed a Boussinesq solver to propagate waves from intermediate depths to nearshore coupled to a Navier-Stokes solver to resolve wave breaking. Alternatively, Tonelli and Petti [190] suggested a Boussinesq model for propagating waves when nonlinearity and dispersion effects are of similar order and switching to a NLSW solver when nonlinearity dominates. The combined Boussinesq and nonlinear shallow water models capture wave dispersion and are attractive for resolving dynamic wave runup and overtopping.

Nonlinear Shallow Water Equations
The nonlinear shallow water equations are a common method of representing wave overtopping. Although the hydrostatic assumption is violated in circumstances where vertical velocity is a significant flow feature, the NLSW equations have proven to be a reasonable approximation in the mid to inner surf zone [191], and represent the majority of overtopping models [121]. Numerous NLSW models have investigated wave runup (e.g., [191][192][193][194]) and overtopping (e.g., [121,136,[195][196][197][198][199]). However, a sensitivity to boundary location was observed in multiple studies and careful placement of the boundary is advocated [34,[199][200][201][202]. XBeach [203] is a widely used 2D flow and sediment transport model that resolves both infragravity and incident wave forcing. XBeach-Surfbeat solves a time dependent wave action balance that forces a Generalized Lagrangian Mean (GLM) formulation of the nonlinear shallow water equations [203]. XBeach-Nonhydrostatic solves the NLSW equations, and a one-layer, nonhydrostatic pressure correction [204] enables short wave variation in intermediate to shallow depths. XBeach has been used primarily to model erosion, overwash or barrier breaching during storm events (e.g., [203,[205][206][207]) and more recently has been used to predict temporally and spatially variable overtopping flows (e.g., [34,73]). XBeach-G [208] is a variant of XBeach developed specifically for gravel beaches, which utilizes the nonhydrostatic flow solver with a groundwater model to simulate infiltration and exfiltration through permeable gravel beds. Similarly, SWASH [130] is a nonhydrostatic shallow water solver which admits multiple layers to account for vertical flow structure and has been used to model overtopping processes [209,210]. Generally, results suggest that NLSW models are capable of predicting mean overtopping rates; however, random wave phasing presents a particular challenge. A number of studies have suggested a sea state of 1000 random waves provides consistent results for a given water level [198,199,202]. Notably, the effect of converging overtopping rates with tidal fluctuations has not been considered.
Overtopping flows are typically transferred to an overland flow solver using an overtopping flow hydrograph from empirical or numerical estimates (e.g., [34,64,73,169]). Although mass is conserved in this coupling, momentum information is lost, potentially impacting the flood description near the overtopped feature and precluding structure wave loading analysis. Elsayed and Oumeraci [73] simulated overtopping, breaching and overflow processes using XBeach and HEC-RAS and reported that hydrograph-transferred flows, which conserve mass but omit momentum, produced higher water depths in the backshores.

Overtopping Flood Validation
Few studies have attempted to validate wave overtopping volumes in relation to flooding. Laudier et al. [132] presented a quantitative field scale validation of overtopping on a natural beach in Central California, and results suggested that empirical models overestimate overtopping rates. Cheung et al. [135] and Lynett et al. [134] presented numerical overtopping models along with qualitative validation data (e.g., high water marks or levee damage) and in the case of Lynett et al. [134], empirical and numerical estimates differed by a factor of ten. Smith et al. [63] considered a urban coastal flood event along the North Somerset coast in the UK and used point sources to introduce overtopping volumes to the flooding domain, however overtopping rates were not modeled in a prognostic manner, but rather from a post event analysis of the inundation area which revealed the flood volume. Moreover, the analysis suggested significant uncertainty in the overtopping estimate, and the study concluded overtopping volumes are a dominant source of uncertainty relative to flood extent prediction. Gallien [34] collected wave overtopping flood validation data and showed that empirical and numerical estimates performed similarly despite order of magnitude differences in total overtopped volume, suggesting that a few large waves in succession (likely resulting from the combined effects of incident and infragravity waves) dominated backshore flooding. Indeed, wave overtopping is considered a significant deficiency in the current modeling methodology [28,121,122], and multiple studies stress the need for field validation data [29,32,62,168,211,212].

Flood Model Evaluation
Generally, in coastal flood modeling systems, individual model elements such as water elevation or wave characteristics are validated against a known event (e.g., [33,64,169]). These evaluations, although useful to accurately represent boundary conditions, do not provide flood prediction information. Flood extent, depth, velocity, arrival time and duration are often qualitatively evaluated. Recent technological advances in lower cost sensors, surveying, remote sensing and unmanned aerial observations systems (i.e., drones and UAV) present the opportunity to improve quantitative evaluation of flooding events. Traditionally, flood extent has been used to validate coastal flood models (e.g., [20,21,27,34,62,63,97]).

Flood Extent
Coastal flood extent model performance has been quantified using overlap [33] or a coefficient of areal correspondence [213]. Overlap represents the fraction of the observed flood that was correctly modeled and is described as, where A O and A M represent the observed and modeled flood extent, respectively. Overlap tests only whether observed flooded areas were predicted and ignores overprediction (e.g., [33]). A fit measure of one shows that all observed flooding was predicted, however does not penalize overprediction, i.e., if a model predicts extensive nonphysical flooding it could still achieve a 1 (i.e, 100%) overlap. A more robust flood prediction evaluation relies on the coefficient of areal correspondence [213], commonly referred to as fit agreement or F<1>, which represents the intersection of modeled and observed flood extents divided by the union of the predicted and observed flood extent (e.g., green in Figure 3) as follows, A fit measure of zero and unity corresponds to no agreement and complete agreement, respectively. Further evaluation of a model can consider a measure of underprediction, F UP , which characterizes the fraction of flooded area observed but not predicted, as follows, In this case, a fit measure of zero and unity correspond to no underprediction and complete underprediction, respectively. Lastly, a measure of over overprediction, F OP , characterizes the fraction of flooded area predicted but not observed, as follows, In this case, a fit measure of zero and unity corresponds to no overprediction and complete overprediction, respectively. Superior models will maximize F A while minimizing both F UP and F OP .

Hydrodynamic Validation
Few studies have attempted to hydrodynamically validate flooding. Laudier et al. [132] conducted field observations of a lagoon-filling event but has drawn criticism in the literature [214]. Cheung et al. [135] and Lynett et al. [134] presented numerical overtopping models along with qualitative validation data (e.g., high water marks or levee damage). Smith et al. [63] presented a model that was validated using photos derived of high water marks and eyewitness data. LeRoy et al. [64] used water depth, reported stagnation area and eyewitness reports. However, there are no known quantitative hydrodynamically field-validated coastal flooding models. Multiple studies stress the need for field validation data [21,29,32,168,211,212]. A paucity of in-situ hydrodynamic validation (depth, velocity, flood arrival time) inhibits accurate urban coastal flood modeling and prediction [34,169].

Multi-Pathway Flooding
Rising sea levels will challenge low-lying coastal communities necessitating increasingly accurate flood prediction. Neglecting compound impacts may significantly underestimate flood risk, particularly in urban areas [19]. One of the key issues in effective coastal flood prediction is representing the combined effects of oceanic and freshwater forcing along with urban infrastructure that constrain or redistribute water, and integrating structural breaching and failure processes. Careful consideration of future flooding source-receptor pathways is recommended, since mitigating one pathway may exacerbate another (e.g., [21]).
Sea level rise will propagate wave energy shoreward, potentially exposing barrier island communities and urban sand spits to wave overtopping flooding in addition to overflow from high water levels found in inlets and bays. Numerical models (e.g., nonlinear shallow water, Boussinesq, Navier Stokes) reasonably simulate temporally and spatially variable overtopping rates and may be coupled to hydrodynamic models that simulate weir-like overflow and propagate overland flooding (e.g., BreZo, Delft3D, and HEC-RAS). Particular attention should be paid to flow transfer when separate overtopping and overland flow models are used for computational efficiency. Ideally, model development will tightly couple these processes in a computationally efficient way.
In coastal applications, storm system drainage is critical to predicting flood extent for three reasons; (1) bay or ocean water levels may flood low protected backshores through the storm drain outfalls; (2) the storm drain system may redistribute water; and (3) overland flow volumes may be reduced by storm system drainage. Storm system and overland flow are commonly considered in 1D-2D coupling using 1D link-node models such as SWMM coupled to a 2D overland flow solver, but the uncertainty of 1D-2D models has not been well addressed [42].
Compared to marine flooding, groundwater inundation modeling is extremely limited [215,216]. Even in areas where sea level rise may not cause direct groundwater inundation, water tables at shallow and intermediate depths have the potential to jeopardize infrastructure [106], and at locations further inland than coastal surface flooding [105,217,218]. Elsayed and Oumeraci [219] presented a coupled modeling approach to resolve overtopping, overland flow propagation and infiltration impacts on coastal aquifers, and suggested further development is needed to resolve the unsaturated zone and, in turn, the surface-groundwater interaction.

Conclusions and Future Work
Nonlinear shallow water models employing shock capturing schemes (e.g., Galerkin and Godunov) have been shown to accurately route overland flow and handle critical flow transitions caused by urban features such as streets, curbs and walls; however, depth and velocity estimates are highly sensitive to terrain and building representation [53,64]. Static ("bathtub") models perform poorly in defended urban backshore and should be avoided.
High resolution LiDAR-derived DTMs generally represent surface topography, however to accurately route flood flows the modeler must be careful to include fine scale hydraulic features which are un-or under-resolved in the LiDAR data (e.g., seawalls and artificial dunes). Building representation significantly affects velocity estimates [53]. Appropriate site representation often requires extensive site knowledge and additional surveying. A paucity of observational data obstructs accurate flood risk assessment and high quality observational data collections are advocated [93,98]. Flood extent is the most common validation method; however, quantitative in-situ hydrodynamic validation (e.g., depth and velocity) data are urgently needed [34,169].
Wave overtopping is a significant deficiency in current modeling methodologies [28,121,122], and is identified as key future research area [123]. Often, the beach itself is poorly resolved or excludes temporary dune structures, both of which fundamentally alter backshore flood predictions. Static methods are sensitive to freeboard (a function of beach elevations) and cannot represent temporally variable overtopping flows. EurOtop and other empirical models, strictly intended for structures but often used on beaches, can provide temporally variable average overtopping estimates (for varying offshore boundary conditions), however do not resolve impulsive swash event volumes. Flood extent and hydrodynamic differences from average versus impulsive overtopping estimates are not considered in the literature and deserve attention. Numerical models can simulate impulsive, temporally and spatially variable overtopping rates and may be coupled to hydrodynamic models that simulate weir-like overflow and propagate overland flooding. Typically, flows are transferred through one-way coupling using a flow hydrograph which conserves mass but does not consider momentum. Depending on the problem of interest, the modeler is left to decide if a computationally efficient mass conservation one-way coupling is sufficient, or if higher-resolution tightly coupled models are needed. Notably, random wave phase uncertainty affects overtopping volume estimates (particularly with dynamic water levels). Wave overtopping flooding represents a significant modeling challenge and field observations are needed [34].
Compound and multi-pathway flooding are critical to accurately characterize future, non-stationary coastal vulnerability. In estuarine areas, water levels may vary significantly through tidal amplification, tide-surge interactions, and fluvial inputs and should be resolved in flood modeling efforts. Pluvial flows have received only limited attention but are a key component in assessing urban coastal flooding. Nested modeling is the current state-of-the-art for representing local water levels, wave overtopping, overland flow and sewer flows. Ideally, model development will incorporate multiple processes to minimize coupling challenges. Deterministic modeling is necessary to identify dominant flooding pathways and proves helpful in considering particular scenarios, but cannot capture the range of potential coastal flooding hazards. Probabilistic frameworks considering uncertainty and coastal process dependence are recommended. Critically, data and models must be synthesized to explicitly represent all potential flooding source-receptor pathways under future conditions.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript: