Modeling the Impact of Sea Level Rise on Maximum Water Elevation During Storm Surge Events; A Closer Look at Coastal Embayments

Storm-surge models are commonly used to assess the impacts of hurricanes and coastal storms in coastal areas. Including the impact of the projected future sea level rise (SLR) in these models is a necessary step for a realistic flood risk assessment. Commonly, SLR is superimposed linearly on the estimated water elevation. This approach, while efficient, may lead to inaccuracies. Further, developing a new model with updated data that include the impacts of SLR (i.e., nonlinear approach) is time consuming. We compare the linear and nonlinear approaches to include the effect of SLR to predict Maximum Water/Flood Elevations (MWE) as a result of storm surge and SLR. After a simplified theoretical analysis, a number of idealized cases based on the typical coastal bodies of water are modeled to assess the impact of SLR on MWE using the linear superposition and nonlinear approaches. Additionally, two case studies are carried out: Narragansett Bay, RI and Long Island Sound, CT (USA). Results show that for the idealized cases with variable depth, in general, the linear superposition of SLR to MWE is conservative (i.e., predicts a larger flood elevation) relative to the nonlinear approach. However, if a constant depth is considered, results are not consistent (i.e. linear superposition can overestimate or underestimate MWE, and the results depended on the geometry). The simulated MWE from the Narragansett Bay simulation confirms the outcome of idealized cases showing that linear assumption is conservative up to 10% relative to the nonlinear approach. For this study, Hurricane Sandy and a Synthetic Storm from the US Army Corps of Engineers North Atlantic Comprehensive Coastal Study (NACCS) dataset are simulated. Long Island Sound model results are also consistent with the idealized case. In general, based on the results of the idealized and real studies, a discrepancy of up to 10% between the linear and nonlinear approaches is expected in estimation of MWE which can be underor over-estimation of flood elevation. Address(es) of author(s) should be given

Keywords Sea Level Rise · storm surge · water elevation · idealized case · hydrodynamic model · ADCIRC

Introduction
The rate of sea level rise (SLR) has been accelerating in the recent century (Church and White, 2006;Douglas, 1992;Hinkel et al., 2018) due to a warming climate which has led to melting of ice, and thermal expansion. While ocean surface waves, tsunamis, storm surges, tides, seasonal cycles, and ocean circulations may result in a relative rise in sea level (Sweet et al., 2017;Williams et al., 2008), the average SLR is associated with warming of climate (Moon et al., 2018). SLR is not uniform over the world (Church et al., 2004). Global mean sea level rise is estimated in the range of 0.3 -2.5 meters until 2100 (Parris et al., 2012). It is expected that the US northeast will experience a high rate of SLR in the order of 1 to 3.5 m (extreme highs; Sweet et al. (2017)). For instance, it is predicted that Boston, MA will have an intermediate SLR scenario of 1.2 m and an extreme SLR scenario of 3.3 m by 2100 (Sweet et al., 2017).
The effect of SLR should be considered for coastal flood risk assessments (Lawrence et al., 2018). The total water level during a flood event is caused by the combined action of tides, storm surge, and wave-setup (Wolf, 2009). Therefore, it is important to understand how SLR can affect these processes.
Several studies have investigated the impact of SLR on astronomical tides at various scales. (Devlin et al., 2014(Devlin et al., , 2017a(Devlin et al., ,b, 2019Friedrichs et al., 1990;Kresning et al., 2019;Schindelegger et al., 2018;Talke and Jay, 2020). For instance, based on a global tidal model, it has been shown that the M2 component (lunar semidiurnal constituent) is sensitive to SLR in the North Atlantic Ocean and many marginal seas. It can change with relative magnitudes of 1-5% per century assuming the adopted SLR scenarios (Schindelegger et al., 2018). In the Pacific Ocean, it has been shown that short-term (monthly to inter-annual) mean sea level variability can lead to inter-annual changes (increase/decrease) in tidal constituents (Devlin et al., 2017b). Consequently, these changes can lead to higher predicted flooding compared with predictions that are based on the stationary tidal components. In addition to tides, nearshore waves and consequently wave-setup in shallow waters can be affected by SLR (Atkinson et al., 2013).
During a storm, a regional ocean model is commonly forced by a combined signal of astronomical tide and storm surge at the open boundaries as well as wind shear stress and atmospheric pressure over the domain. Consequently, the interaction of SLR and MWE is more complex (e.g., very nonlinear signal at open boundaries and wind stress over the domain) during storm conditions compared with calm periods when MWE is caused by astronomical tide alone: the interaction of astronomical tides and SLR has been studied in detail by representing the tidal signal as a time series of harmonic constituents with predefined periods. For instance, the change in the resonance properties of a bay or an estuary due to SLR (i.e., depth change) may explain the nonlinear interaction of MWE and SLR when only astronomical tides are considered. (see Friedrichs (2010), Talke and Jay (2020), and Section 3.3.2). Therefore, more analytical and numerical research is necessary to investigate the interaction of SLR and MWE when storm surge is added to the tidal signal..
For practical flood risk assessments, if the impacts of the projected SLR on storm surge and tides are neglected, SLR can be linearly superimposed to the current estimated MWE. Consequently, the effect of SLR on MWE can be conveniently estimated for flood risk assessment. The linear superposition approach is attractive since it does not require additional simulations and is simple/fast to implement . In contrast, another numerical model that includes changes in the bathymetry and/or boundary conditions due to SLR can developed (i.e., nonlinear approach). Both linear and nonlinear approaches are common in flood risk assessments. For instance, the hydrodynamic response of estuaries and coasts (i.e., storm surge and waves) under combined SLR and storm events was studied in Ding et al. (2013). It was shown that for the Gulf of Mexico under concurrent Hurricane Gustav and hypothetical SLR scenarios, changes in the MWE were linearly proportional to the SLR for most of the regions. In other studies, it was shown that it is not appropriate to linearly add SLR to the design water levels in southeastern Louisiana (Smith et al., 2010). They showed that in the areas with maximum surge, the MWE generated under the SLR scenarios can be estimated by linear superposition (i.e., adding SLR). However, in the areas with moderate peak surges and in isolated areas, the storm surge increased much more than those from linearly adding the SLR.
The North Atlantic Coast Comprehensive Study (NACCS) simulated winds, waves, and water levels for 100 real extratropical and 1050 synthetic tropical storm events for the US East-coast from Virginia to Maine (Cialone et al., 2015). SLR scenarios were considered in the NACCS study.
In this research, the impact of SLR on MWE is investigated by comparing the linear superposition and nonlinear approaches (in the linear approach, the SLR is simply added to estimated MWE in the current condition whereas in the nonlinear approach, simulations are carried out based on the changed model input data such as bathymetry assuming a SLR scenario). While both methods are subjected to model uncertainty, we assume that nonlinear method is more accurate as it is subjected to less model error. The main objective of this study is to estimate the increased level of uncertainty or error when the linear superposition method is applied. Moreover, the focus of this study is storm surge; therefore, it is assumed that the tidal components remain stationary at the model domain boundaries under SLR scenarios; this assumption can lead to additional uncertainty. At first, a theoretical analysis based on the simplified shallow water equations will be presented. MWE will be calculated for a few typical idealized cases using the two approaches. Two case studies were also considered and discussed: Narragansett Bay, RI and Long Island Sound, CT (USA). Finally, discussion and conclusion have been provided as a summary of analytical, idealized, and real case studies.

Simplified theoretical analysis
The effect of SLR on water elevation can be theoretically studied using a simplified shallow water or 1-D Saint-Venant equations. These equations are commonly used to simulate long waves (for instance, storm surge or tides) and in flood simulation studies of rivers and coastal zones. Shallow water equations are derived from the Navier-Stokes equations by integrating over the depth assuming that the vertical acceleration is small, and flow is nearly horizontal. A number of mathematical representations can be used for shallow water equations based on different hydrodynamic parameters (for example, velocity or flow rate as the state variable; Lai et al. (2002); Montes (1998)). The simplified continuity and momentum equations for a wide rectangular channel (i.e., channel width is much larger than the water depth) and using water surface elevation and velocity as state variables may be written as (Friedrichs, 2010;Lai et al., 2002), where η is the water level above still water level (see Fig. 1), A is the cross sectional area, U is the depth averaged velocity in the x direction, B is the width of the channel, h is the channel depth, n is the Manning coefficient, ρ is the density of water, ρ a is the air density, C is wind drag coefficient, V w is the wind velocity in the x direction, Θ = θ w (x, t) θ c (x), θ w is wind direction, and θ c is the direction of the channel centreline at a particular point x. The terms of Equation 2 are the inertia force ( ∂U ∂t ), the pressure force (g ∂η ∂x ), the friction force (g n 2 U |U | The convective acceleration is ignored for simplicity. This is a reasonable approximation in many storm surge modeling and flood mapping applications, for instance, in Sea, Lake, and Overland Surges from Hurricanes (SLOSH) model (Jelesnianski et al., 1992). It has also been shown that the bottom friction is usually much more significant (sometimes 10 times more) than the convective acceleration in the momentum equation for simulation of long waves (Jones and Davies, 2008). The convective terms are significant when there is a sharp curvature in the geometry or bathymetry of the area (Li et al., 2008). As the analysis is 1D, the Corilols forces are not considered. Applying SLR to the equations of motion leads to a change in the still water depth h. Here, we investigate which terms of these equations are affected and whether these effects are linear or nonlinear. Therefore, to account for the SLR in Equations 1 and 2, water depth was increased by h 1 (SLR) (i.e., the mean water or still water depth is increased to h + h 1 (see Fig. 1). Here, it is assumed that SLR, h 1 , is uniform across the domain which may not be always the case. (Familkhalili and Talke, 2016). It is also assumed that the width of the channel, B, is much larger than the total depth and variations of width with the total depth is neglected. Consequently, Eq. 1 and 2 will be changed to, in which prime notation " 0 " indicates the values for the nonlinear approach. η 0 shows the calculated storm surge assuming a nonlinear approach. It is clear that if η 0 = η, the linear and nonlinear approaches result in the same water elevation (i.e, h + η + h 1 ). To simplify the new nonlinear friction and wind stress terms in Eq. 2, the Taylor series expansion was applied to the denominator of the each term as follows, Expanding Eq. 1 and substituting Eq. 5 and Eq. 6 into Eq. 4, the modified continuity and momentum equations can be written as, The modified equations of the motion include four additional terms: two terms in the continuity equation and two terms in the momentum equation. The term I in Eq. 7 is proportional to the convective acceleration which is not that significant in many storm surge modeling applications. Nevertheless, sharp gradients in water depth can lead to a significant spatial gradient in velocity or ∂U 0 ∂x .T e r mII in Eq. 7 is generated by a change in the geometry of the channel. For instance, in a bay with a variable width (funnel shape) this term can be significant. When the width of the water body (B) is constant, this term is zero. If the friction and convective terms are neglected in the momentum equation (e.g., in deep waters), the geometry is the main source of discrepancy between linear and nonlinear method. Many estuaries or bays have variable width. In the following sections, using idealized and real case studies, more complex water bodies will be assessed to investigate the impact of these terms as sources of errors..
Referring to Eq. 8, terms III and IV are potentially decreasing the bottom friction and wind force terms due to SLR (while the wind stress does not change by water depth the wind stress term depends on water depth). To better quantify this term, ξ and ζ are defined as the relative magnitude of these additional terms to the friction and wind terms, Fig. 2 shows how ξ and ζ vary with SLR scenarios and water depth. For instance, in 10 m water depth and assuming 1 m SLR, these additional terms lead to around 14% and 10% reduction in friction and wind-induced terms, respectively. As storm surge is important in shallow water zones, it is clear that these terms can be significant in propagation of storm surge in flood zones. Term III also leads to changes in the velocity field (e.g., change in current speed due to reduced bottom friction). While this simplified theoretical analysis provides some hypotheses about the discrepancy between linear and nonlinear approaches, as discussed above, more detailed analysis will be provided in the idealized and real case studies in the following sections. In the next section, three idealized cases will be introduced. The SLR will be applied to each case in order to estimate the water elevations using the linear superposition of SLR and the nonlinear method. The presented results will show how the potential sources of nonlinearity derived in this section, will affect the MWE in each scenario.

Model description
To simulate the MWE for the idealized and real cases, a popular numerical model (Luettich Jr et al., 1992), ADvanced CIRCulation (ADCIRC) was implemented. ADCIRC 2D solves the depth averaged shallow water equations for the free surface circulation based on the hydrostatic pressure and the Boussinesq approximations. ADCIRC uses the Galerkin finite-element method for spatial discretization and a three-level finite-difference method to estimate the temporal derivatives. In ADCIRC, the surface water elevation is computed using the continuity equation. The depth-integrated current velocities are computed from the momentum equations. ADCIRC is usually forced along the open ocean boundaries by water elevation/tides and surface wind stress and pressure over the model domain.

Idealized cases
Study of idealized cases (Huang et al., 2008;Pandoe and Edge, 2004) helps better understand and generalize the impact of SLR on the MWE for similar geometries. For idealized cases a constant 2 m SLR was applied. This SLR approximately corresponds to 1.9 m SLR for the Intermediate-High SLR Scenario (Sweet et al., 2017) in Boston, MA (USA).
Three idealized cases were simulated. In terms of dimensions and the geometry, these cases were generated to mimic the realistic dimensions of coastal embayments. However, there is a large variation among water bodies and coastal embayments around the world, and these selections represent a sample. Fig.  3 show the details of the idealized cases, which will be referred as type I, II, and III.
For all three cases, two scenarios for the depth were considered: a) variable depth of 100 m at the ocean boundary reducing to zero at the coastline assuming a constant gradient; b) constant depth of 20 m. In the constant water depth scenarios, we assumed that the depth linearly decreases to zero at the inland boundary in the cross shore direction and applied wetting and drying boundary condition. Further, for modeling an estuary that is connected to a river, it is more realistic to apply a radiation open boundary condition (IBTYPE(k)=30 in ADCIRC), which allows a long wave to propagate out across the intersection of the river and the estuary (although for the idealized case type III, considered in this study the results were less than 2% sensitive to changing the boundary condition from wetting and drying to radiation.)

Model setting
The time series of water elevation associated with Hurricane Sandy was used to force the model at the ocean boundary. Hurricane Sandy (Oct. 2012) was selected because it was a slow moving storm that severely impacted the US east coast from Florida to Maine with a long duration (interacted with several tidal cycles). The sample time series of the water elevations for this model was extracted from a previous study (Hagen et al., 2006) 140 km offshore of Atlantic City (the landfall of Hurricane Sandy), where the depth of water is 80 m. Fig. 3d shows the time series of water elevations used to force the ocean boundary nodes from Oct. 28 to Oct. 31, 2012. The three idealized cases were run during the period of Oct. 28 to Oct. 31, 2012 for three days. The model ramp up time was 12 hours, and the time step was 0.5 seconds. The bottom friction for all models was set based on the Manning coefficient n of 0.018 s/m 1/3 for sandy beds (Arcement and Schneider, 1989). ADCIRC model can use both Manning coefficient and quadratic friction drag coefficients: C d = gn 2 / 3 p h). The weighting factor (Tau0) in the Generalized Wave-Continuity Equation (GWCE), which weighs the relative contribution of the primitive and wave portion of the GWCE was set to 0.005 in the deep water and 0.02-0.03 in shallower waters according to the ADCIRC's User Manual recommendation. Typical values of Tau0 are between 0.005 to 0.1 (Luettich and Westerink, 2007). The advective terms in the momentum equation and the wetting/drying (Dietrich et al., 2004) options were included in the simulations.

Results of the idealized cases
The linear superposition and the nonlinear methods were used to compute the MWE for the idealized cases. As mentioned, in the linear approach, the SLR is simply added to MWE. In the nonlinear approach simulations were carried out based on the changed bathymetry (i.e. the depths were uniformly increased over the entire mesh) which includes SLR. For each scenario, two bottom friction formulations were implemented: a constant bottom drag coefficient and the Manning formulation. Due to the similarities of the results, only the results that are based on the Manning coefficient are presented. Fig. 4 shows the percentage difference of the MWE between the linear and the nonlinear approaches. The percentage of difference, or error of the linear superposition method, is calculated as follows, Er = (max water elevations) L (max water elevations) NL (max water elevations) NL ⇥ 100 (11) in which the subscripts L and NL represent the linear and nonlinear methods. It should be noted that the Er only estimates the water elevation error and does not necessarily represents the error in flood inundation area. Depending on the topography, a small change in elevation can lead to very large changes in the flood extent or inundated areas (Orton et al., 2015(Orton et al., , 2019. As Fig. 4 shows, for the idealized cases type I and II with constant depth (Fig. 4a,c), the MWE based on the nonlinear method near the coastline are up to 7% higher than those from the linear method (negative error). In contrast, for the idealized case type III (Fig. 4e), the MWE based on the linear superposition are higher up to 3.5% (positive value). Considering a variable depth, Fig. 4b for type I, the MWE estimated from the linear superposition of SLR are up to 10% higher than those from the nonlinear approach (in contrast to Fig. 4a). Referring to Section 2, this contrast shows the relative importance of change in width and bathymetry of the domain. For type II case with a variable depth (4d), the MWE based on the linear method are up to 5% lower near the coast consistent with the constant depth case. However, the linear approach overestimates the MWE in other parts of the domain. Finally, for type III case with variable depth the error oscillates between positive and negative values. Based on these results, for the variable depth scenarios, the linear method is a conservative approach for the idealized case type I while it is underestimating the MWE for type II and II cases by 6%. Further, it was shown that the maximum expected error of the linear superposition method is about 7% in nearshore areas which is in the range of the accuracy of the storm surge models (Zhong et al., 2010).
In addition, the effect of wind stress was investigated for the idealized case type I. The model was forced by 10 m/s and 20 m/s wind velocities. Wind velocities were uniform over the entire mesh in the direction shown in the Fig. 5. Similar to the case with no wind stress in Fig. 4(b), results show the overestimation of the MWE in the linear method up to 10% for 10 and 20 m/s wind velocities (Fig. 5a, 5b) It should be mentioned that the effect of bottom friction term is more important than that of the wind stress term, which is consistent with the outcome of theoretical analysis.

Narragansett Bay
Rhode Island has been impacted by several hurricanes over the past decades including Hurricane Bob (1991), Hurricane Irene (2011), and Hurricane Sandy (2012). To protect the downtown of Providence from the storm surge, a hurricane barrier was constructed in 1960s on the Providence River (Morang, 2016). Therefore, prediction of flood risk with inclusion of impacts of the SLR is of interest of coastal communities and coastal managers in this region as reported in several other studies (Grilli et al., 2017;Spaulding et al., 2016).
For this case study, a real historical storm and a synthetic storm were selected. Hurricane Sandy in 2012 led to major economic loss in the Northeast region although the its track was farther away in Atlantic City, NJ. (Xian et al., 2015) The Synthetic Storm 492 from the NACCS database was selected as it approximately produces the same water level as a 100-yr event at the available NOAA (National Oceanic and Atmospheric Administration) water level stations in Newport and Providence. Several data sources were used to simulate this case study in the ADCIRC model. The North East Coastal Ocean Forecasting System (NECOFS) for the northeast of the US is used as the base for the model domain (Chen et al., 2003;Torres et al., 2018). The NECOFS mesh (Fig. 6a) was enhanced in the Narragansett Bay, Rhode Island to improve the accuracy of the simulation within the bay, and used as the current ADCIRC mesh. It includes 105,560 nodes with the finest resolution of 20 m in the nearshore zones. The bathymetric and topographic data with resolution of 30 m and 1m, respectively, were  More details about the ADCIRC model setting including forcing is provided here and in the previous published paper (Torres et al., 2018): The model was forced by M2, S2, N2, K1, and O1 tidal components at the ocean boundaries and hurricane wind fields over the domain. The model was run under two scenarios: during Hurricane Sandy from Oct. 28 to Oct. 31, 2012, and during a Synthetic Storm 492 from the NACCS dataset (Nadal-Caraballo et al., 2015). The time step was 0.5 s for both models. The Synthetic Storm 492 was selected as it approximately produces the same water level as a 100-yr event (at the upper 95% confidence limit) at the available NOAA (National Oceanic and Atmospheric Administration) water level stations in Newport and Providence. For instance, at the Newport Station 8452660, the water elevation for the upper 95% limit of the 100-yr storm is estimated 2.9 m, and the MWE from the NACCS Storm 492 at the same location is 3.2 m. To generate the wind field for this synthetic storm, the ADCIRC parametric asymmetric Holland model was used (Dietrich et al., 2017). The Synthetic Storm 492 parameters (track, intensity, and forward speed) were extracted from the NACCS dataset. The wind forcing for Hurricane Sandy was obtained from the meteorological datasets of the NECOFS Weather Research and Forecasting (WRF). For each hurricane, two cases were simulated: 1) using the current bathymetry of the domain and 2) adding 2 m to the bathymetry throughout the domain. The former scenario was used to simulate MWE based on the linear superposition of 2 m SLR, and the latter was used to produce the nonlinear results. Fig. 7 shows a map for the error of the linear superposition approach (Eq. 11) in the Narragansett Bay for Hurricane Sandy (Fig.7a) and the Synthetic Storm 492 (Fig.7b). Areas with blue color are regions where linear superposition of SLR is conservative (i.e. positive error, linear results higher than nonlinear results). The ADCIRC results for both Hurricane Sandy and Synthetic Storm 492 simulations show 4 -10 % overestimation using the linear superposition method. The results are consistent with the type I idealized case shown in Fig. 4b. Both indicate that linear approach is conservative and overestimates the MWE by 10 %.

Long Island Sound
Long Island Sound is a large estuary in the US Northeast between Connecticut State to the north and Long Island to the south. The total area of the sound is 3350 km 2 with 960 km of coastline. More than 90 % of the Connecticut population live along the coast of this sound (Soundwater, 2019). New York City is located at the western end of the sound. . Long Island Sound has experienced several devastating storms including Hurricane Sandy (Georgas et al., 2014) and Hurricane Carol in 1954(Colle et al., 2008. Tides in Long Island Sound are dominated by the M2 constituent (semidiurnal lunar with a period of 12 hr and 25 minutes). Assuming that the width and depth of a channel are constant and friction ignored, tidal resonance occurs when the length of a channel, L, is close to the quarter of tidal wavelength(λ); i.e., L = λ M 2 /4=T M 2 p gd/4w h e r e p gh is the tidal wave speed (Talke and Jay, 2020). Therefore, SLR (i.e., change in depth) can affect the resonance properties of a basin. It has been shown in various studies that Long Island Sound is a near resonance system with nearly fourfold amplification of the semi-diurnal (M2 and S2) tides within the sound (Godin, 1988;Wong, 1990). Several studies have assessed changes in the resonance properties of Long Island Sound due to SLR. It is estimated that tides in Long Island Sound are amplifying by SLR at the rate of 10% per meter (Kemp et al., 2017;Talke and Jay, 2020). This effect will further amplify the impact of SLR on MWE. Similar to the previous case study in the Narragansett Bay, the NECOFS domain was re-meshed with a better resolution along Long Island Sound. The updated mesh contains 139,351 nodes. The bathymetric data were obtained from the interpolated bathymetric/topographic data of 1/9 arc second (3 m) resolution from NOAA's Continuously Updated Digital Elevation Model (Eakins et al., 2015). Fig. 6 shows the domain of the simulation in which the finest resolution of the grid is 200 m in the shallow regions.
The ADCIRC model was forced by M2, N2, S2, K1, and O1 tidal constituents at the boundary. Manning coefficient was assumed 0.018 s/m 1/3 for deeper water and 0.05 s/m 1/3 in the shallower zones. The wind forcing for hurricane Sandy was extracted from the European Center for Medium-Range Weather Forecasts (ECMWF) model (Tomassini et al., 1998). The ECMWF implements a number of physical parameterizations in a combined data-assimilation and circulation model for real time climate analysis. It should be noted that similar to Narragansett Bay case study, the NECOFS wind was initially used for simulation of Hurricane Sandy. However, as NECOFS system is focused on the US northeast, it was shown that ECMWF wind (Torres et al., 2018) resulted in better model performance (Fig. 8). The ADCIRC model was run under two scenarios: 1) Hurricane Sandy from Oct, 28 to Oct. 31, 2012 and 2) Synthetic Storm 492 from NACCS for an arbitrary time period of three days from July 12th to July 15th, 2000 with a time step of 0.5 s for both models. The Synthetic Storm 492 generated approximately 100-yr storm for this case study. At the New Haven Station 8465705, the water elevation for the upper 95% limit of the 100-yr storm is estimated 2.8 m, and the MWE from the NACCS Storm 492 at the same location is 3.3 m.
Four NOAA water level gauges at Kings Point, NY (8516945), Bridgeport, CT (8467150), New Haven, CT (8465705), and New London CT (8461490) that are located inside the sound were used for model validation (Fig. 6).
The time series of the observed water level were compared to those from the ADCIRC model during Hurricane Sandy and are shown in Fig. 8. The comparison of observed and simulated water levels (based on ECMWF) shows a relatively moderate model performance with some discrepancies as follows: there are similar phase errors between model and observation at peak water elevations in all four graphs. These discrepancies can be associated with errors in the wind forcing. Also, model results show errors in prediction of tidal variations, which is clearer during the first cycles of the simulations. Errors in estimation of the tidal component can affect the assessment of linear and nonlinear approaches, which is the objective of the simulations. The model demonstrates more accuracy in terms of the magnitude of MWE (regardless of timing) and in general moderately captures the physics of tides and storm surge in the domain. Therefore, its performance was considered acceptable for the purpose of this study. Fig. 9 shows the map of the error of the linear superposition method (i.e. Eq. 11) in Long Island Sound. Results show that the linear superposition leads to slightly higher MWE near the western end of the sound (i.e. close to NY City) for hurricane Sandy. In contrast, the linear superposition underestimated the surge up to 10% in the eastern part of the sound along the coastlines of CT for both Hurricane Sandy and Synthetic Storm 492. Nevertheless, the difference is less than 5% in the majority of the study area. As the water depth does not change significantly within Long Island Sound (Fig. 6), this case study can be compared with the first or second idealized case with a constant depth (Fig. 4a/c). It can be observed the error of the linear superposition method for MWE is less that 4% (underestimation) for both cases consistent with the real case study.

Discussion
Referring to Eq. 7 and 8, it was shown that in addition to a nonlinear term in the momentum equation (i.e., friction term), two nonlinear terms were generated in the continuity equations which lead to a discrepancy between the linear superposition and nonlinear methods. Term I is a function of velocity gradient. Fig. 10 shows the contour map of the velocity gradient for type I case with variable depth. This figure shows that there is a sharp gradient of velocity in the shallow water zone. Therefore, it is expected to see more error in the linear superposition method in shallower areas and flood zones which is confirmed by Fig. 4b compared with Fig 4a. Additionally, velocity gradient leads to another nonlinearity in the momentum equation due to convective acceleration which was not considered in the theoretical analysis for simplicity (i.e., u ∂u ∂x ). Based on the three idealized cases it was shown that the error of linear superposition was up to 10%. In general, the error in models with a constant depth were less (5%) compared with variable depth cases. It was also shown that the linear superposition method can overestimate or underestimate the nonlinear method depending on the geometry and bathymetry of a case study.
Both case studies in Narragansett Bay and Long Island Sound showed consistency with the results of the idealized cases in terms of the error of the linear superposition method. Simulation of a historical Hurricane Sandy storm  and Synthetic Storm showed that linear method overestimated MWE by 10%. For Long Island Sound, the linear method underestimated MWE by around 5%. It was shown that the shape of water body and its bathymetry leads to these difference.
In the two real case studies, the numerical models can be further improved to reduce model uncertainty. The freshwater riverine discharges have been neglected at the model boundaries. Wet hurricanes lead to compound flooding where some areas are flooded by the combined effect of riverine and coastal flooding. Particularly, for simulation of street-level flooding near estuaries, the freshwater input boundary condition should be included (Liu et al., 2020). However, in regional storm surge modeling studies, the riverine input discharge are sometimes neglected (Marsooli and Lin, 2018). For the Narragansett Bay case study, a sensitivity analysis was carried out in which an extreme input freshwater discharge of 500 m 3 /s (around a 500-year event; Kouhi et al. (2020)) were implemented at the Pawtuxet River boundary. While the level of flooding around the estuary changed, no bay wide effects at the Newport or Providence water level stations were observed. Therefore, for the purpose of this study, which is the comparison of linear and nonlinear approaches at a regional scale, this simplification in the models should not significantly affect the results.
Nevertheless, neglecting the freshwater riverine input is a source of error in some areas.
In addition to SLR, morphological changes in bays and estuaries will impact MWE. This becomes even more complex as SLR plays an important role in morphological changes of estuaries. It has been shown that shallow and deep estuaries will respond in different ways (Dissanayake et al., 2012;Leuven et al., 2019;Van Goor et al., 2003). While this study ignored changes in the bathymetry over long terms, this should be considered as another source of uncertainty for MWE predictions in future.
Over the past decades, tidal dynamics in estuaries have been extensively studied using theoretical and numerical methods. These studies have further analyzed the nonlinear interactions of the terms discussed in Section 2 (Friedrichs, 2010). For instance, the distortion of tidal signal is a result of the interaction of friction term and the geometry which can eventually impact MWE (Friedrichs and Aubrey, 1988). Further theoretical/mathematical analysis can be carried out to investigate the nonlinear impacts of tidal distortion, storm surge, and SLR.
Based on the result of the idealized and real case studies, it was shown that the error of the linear superposition method in estimating MWE were generally up to 10%. Nevertheless, it should be emphasized that in localized flood impact studies in which wave-surge-structure interaction is simulated (Park et al., 2018), the error can be much larger and the effect of SLR on the base flood elevation should be considered in numerical simulations.

Conclusion
Two methods were compared to account for the impact of the SLR on the estimation of maximum water elevations (MWE): linear superposition in which SLR is simply added to model results and nonlinear approach where a new model with updated data is used to include SLR. Based on a simplified theoretical analysis, it was shown that the nonlinear terms in the shallow water equations initiates from the gradient of current velocity (which depends on the gradient in the bathymetry), changes in the width of the water body, and the bottom friction. Therefore, it was expected that in a water body where the water depth is nearly constant or very large and the shape is regular, the linear superposition method does not lead to a significant error. These hypotheses were tested using numerical simulation.
Three idealized cases representing geometries of various coastal embayments were studied (Fig. 3). In all cases, two scenarios of constant and variable depth were investigated. Considering a variable depth for the idealized case type I, the MWE estimated from the linear superposition of SLR are up to 10% higher than those from the nonlinear approach. For the idealized case type II with a variable depth, the MWE based on the linear method are up to 5% lower near the coast consistent with the constant depth case. However, the linear approach overestimates the MWE in other parts of the domain. Fi-nally, for the idealized case type III with variable depth, the error oscillates between positive and negative values. Based on these results, for the variable depth scenarios, the linear method is a conservative approach for the idealized case type I while it is underestimating the MWE for the type II and type III cases by 6%. In general, for the constant depth scenario, the error of the linear superposition method was less than 5%. The error increased to 10% when the bathymetry were variable. The discrepancy were generally larger in shallow water and nearshore areas.
Additionally, two real case studies in Narragansett Bay and Long Island Sound (USA) were conducted. In Narragansett Bay, the storm surge was simulated for Hurricane Sandy and a Synthetic Storm which generated a 100-yr storm in RI waters. Results of the simulation showed that the linear superposition method overestimates the MWE by 10% for both Hurricanes. The results were consistent with those from the corresponding idealized case (type I). In Long Island Sound case study, Hurricane Sandy and a synthetic storm were simulated. The error of the linear method was up to 5% where linear method underestimated MWE.
It was shown that the difference of the linear and nonlinear methods to include the impact of SLR on MWE are up to 10% in regional storm surge models for all idealized and real cases. Depending on the shape of the water body, the linear superposition method can overestimate or underestimate MWE. Because this error is within the range of the uncertainty of regional storm surge models, linear superposition can be used for a rapid estimation of the impacts of SLR on storm surge risk. However, for localized studies in which the propagation of flood (surge and wave impacts) around and through coastal structures in flood zones are simulated, linear superposition method may lead to significant errors (e.g., depth induced wave breaking over structures).
More work in necessary to explore and understand the nature and sources of nonlinear interaction of tides, storm surge, and SLR.

Funding
This work was partially funded by the Northeaster Regional Association of Coastal Ocean Observing Systems (NERACOOS) Coastal Resilience Project (0005266). The authors gratefully acknowledge Changsheng Chen from UMASS Dartmouth for providing the NECOFS FVCOM GOM4 mesh and WRF wind data for Hurricane Sandy, and the National Oceanic and Atmospheric Administration (NOAA) for supplying hurricane and water level data.

Conflicts of interest
The authors declare that they have no conflict of interest.

Availability of data and material
The data that support the findings of this study are available from the corresponding references listed in the manuscript.

Code availability
Not applicable