Investigating coastal backwater effects and flooding in the coastal zone using a global river transport model on an unstructured mesh

Coastal backwater effects are caused by the downstream water level increase as a result of elevated sea level, high river discharge and their compounding influence. Such effects have crucial impacts on floods in densely populated regions but have not been well represented in large-scale river models used in Earth system models (ESMs), partly due to model mesh deficiency and oversimplifications of river hydrodynamics. Using two mid-Atlantic river basins as a testbed, we perform the first attempt to simulate the backwater effects comprehensively over a coastal region using the MOSART river transport model under an ESM framework, i.e., Energy Exascale Earth System Model (E3SM) configured on a regionally refined unstructured mesh, with a focus on understanding the backwater drivers and their long-term variations. By including sea level variations at the river downstream boundary, the model performance in capturing backwaters is greatly improved. We also propose a new flood event selection scheme to facilitate the decomposition of backwater drivers into different components. Our results show that while storm surge is a key driver, the influence of extreme discharge cannot be neglected, particularly when the river drains to a narrow river-like estuary. Compound flooding, while not necessarily increasing the flood peaks, exacerbates the flood risk by extending the duration of multiple coastal and fluvial processes. Furthermore, our simulations and analysis highlight the increasing strength of backwater effects due to sea level rise and more frequent storm surge during 1990–2019. Thus, backwaters need to be properly represented in ESMs to improve the predictive understanding of coastal flooding.

1 Introduction

Backwater zones are regions at the river downstream sections, a fluvial–marine transition area between upstream flow and estuary or coastal river plume, where the river flow is affected by the coastal processes, such as sea level changes, tides and storm surge (Lamb et al., 2012), and can extend hundreds of kilometers upstream in low-lying watersheds (e.g., up to 500 km in the Mississippi River). Coastal backwaters are created by elevated sea level that can cause upstream propagation of flood waves and the attenuation of the spatial and temporal water stage fluctuations (Luo et al., 2017). These effects play a critical role in floodplain storage and river discharge (Paiva et al., 2013) and also have a key influence on the biogeochemistry and geomorphology at the terrestrial–aquatic interface (Dykstra and Dzwonkowski, 2020; Lamb et al., 2012; Ward et al., 2020). With population growth near coastal regions (Tellman et al., 2021), coastal backwaters are expected to exert a greater impact on human and natural systems.

The backwater zones usually face severe flood risks as a result of tide, storm surge, rainfall runoff and their combined effects. During a landfalling hurricane with strong winds and heavy rainfall, storm surge drives coastal waters to propagate into the river network and interact with high river discharge (Bilskie and Hagen, 2018). When multiple drivers occur simultaneously or in close successions, the flood event is referred to as “compound flooding” (Santiago-Collazo et al., 2019). Coastal backwater-induced floods have strong temporal and spatial variabilities (Hendry et al., 2019), depending largely on the local topography and storm characteristics (Gori et al., 2020). Due to climate warming, the frequency and intensity of such compound flooding have exhibited an increasing trend (Bates et al., 2021; Rahmstorf, 2017), as a result of intensified storm surge (Camelo et al., 2020; Marsooli et al., 2019), more frequent extreme precipitation (Alfieri et al., 2016) and accelerated sea level rise (SLR) (Kulp and Strauss, 2019; Orton et al., 2019). Although SLR and storm intensification are considered the most influential flooding drivers (Hwang et al., 2020), projected changes in river discharge also play an important role in modulating the flood potentials (Bermúdez et al., 2021).

Understanding the backwater drivers is prerequisite to mitigating the related flood risks. However, the interactions among the backwater drivers and their respective contributions through fluvial processes, storm and climate are not well understood (Dykstra and Dzwonkowski, 2021). Although there is extensive literature that addresses the storm surge-induced coastal inundation (or flooding on land) and the related impacts on flood risks in coastal cities (Hinkel et al., 2014; Ye et al., 2020), limited efforts are made toward understanding the extreme surge that propagates into the river network (Ikeuchi et al., 2017). The latter is more critical in low-lying mega-delta regions that reside over 0.5 billion people globally (Syvitski and Saito, 2007). Streamflow in the backwater zones is affected by river topology, upstream discharge, sea level variations and their interactive effects (Castelltort et al., 2020; Hellmers and Fröhle, 2022). Specifically, river topology is characterized by the river channel geometry, riverbed elevation and the river's receiving water body. Among these factors, riverbed elevation, due to its control on the backwater propagation extent, has been widely recognized in previous studies (Gori et al., 2020). By contrast, the river's receiving water body has not drawn much attention to date. Since rivers contribute to a variety of water bodies including deltaic floodplains, estuaries and coastal oceans (Mikhailov and Gorin, 2012), the interactive effects of river discharge and sea level vary substantially. For example, the impact of river discharge on the local sea level is much more intense in a narrow tidal-river estuary than in an open sea (Rayson et al., 2015; Chegini et al., 2022). In a narrow estuary, flood risks are further exacerbated because high discharge increases the local sea level, high sea level induced by storm surge impedes river discharge to the ocean, and the interaction of these two mechanisms intensifies the backwater effects (Eilander et al., 2020).

Large-scale river models are one of the major components of Earth system models (ESMs) that couple the atmosphere, land, river and ocean models to simulate the global water cycle (e.g., Golaz et al., 2019; Leung et al., 2020) and assess flood risks (Hirabayashi et al., 2013; Towner et al., 2019). Although hydraulic or hydrodynamic models were used more often in previous studies to simulate storm surge-induced coastal inundation (Bakhtyar et al., 2020; Muñoz et al., 2020), there have been growing applications of large-scale river models to assess the compound fluvial and coastal flooding at basin (Chen et al., 2013), regional (Ikeuchi et al., 2017; Yamazaki et al., 2012) and global scales (Eilander et al., 2020; Mao et al., 2019) because they are more computationally efficient for a large spatiotemporal assessment. The long-term evolution of flood drivers and risks can be quantified in the context of climate change. Moreover, such models, when directly coupled with other components of ESMs, can also provide estimations of energy, biogeochemical and sediment processes that are often neglected in pure flood inundation models (Li et al., 2022). However, several limitations in the current generation of ESMs impair the realistic representation of coastal backwaters and human–land–river–ocean interactions at the terrestrial–aquatic interface (Ward et al., 2020). First, most ESMs are configured with one-way coupled river and ocean models, in which water only flows from rivers to oceans and the impact of elevated sea levels on upstream river stage is ignored. Second, the meshes used in most ESMs are too coarse to represent backwater effects. For example, in the high-resolution configuration of the Energy Exascale Earth System Model (E3SM), a uniform resolution of 12.5 km is used for the river model (Caldwell et al., 2019). The resolutions of other widely used large-scale river models also only range from 5 to 25 km. Much higher spatial resolutions (∼ km) are required to resolve the smaller-scale topology near the coastline (Bates et al., 2021; Trigg et al., 2016) for coastal backwaters. Last, while most existing ESMs apply structured meshes in their river components, unstructured meshes are needed to achieve more flexible variable resolutions within areas of interest, such as high resolutions along the coastline, as well as to accommodate the high spatial variation of coastal processes.

Motivated by the increasing flood risks in a warming climate, this study is part of a larger effort to develop capabilities in representing land–river–ocean interactions in E3SM for modeling the changing compound flood risks in coastal regions and the potential implications for the regional and global water and biogeochemical cycles. More specifically, the objectives of this study are to (a) assess the capability of two-way coupled river and ocean models on a regionally refined mesh to capture coastal backwater effects, and (b) understand the major and interactive backwater drivers and their long-term variations under climate change in two contrasting coastal river basins. The backwater drivers are decomposed using a novel extreme flood event selection scheme. Each selected event is identified by the dominant flood drivers. In Sect. 2, we provide an overview of the study domain, the river model, the unstructured mesh, and the methods of extreme event selection and drivers decomposition. Model evaluation and analyses are provided in Sects. 3, 4 and 5. In Sect. 6, we discuss the findings and limitations. Finally, the conclusions are provided in Sect. 7.

2 Methods and data

This section describes the study domain of two mid-Atlantic river basins. A global river routing model on a regionally refined unstructured mesh with two-way river–ocean coupling physics is introduced. We also develop a method to select extreme events and decompose the flood drivers of the selected events. Two hurricane events are selected to demonstrate the applicability of the proposed methods.

2.1 Study domain

The mid-Atlantic region of the United States is exposed to frequent tropical cyclones that bring intense precipitation and storm surge (Sun et al., 2021). In this study, we define the mid-Atlantic region as Susquehanna River Basin (SRB) and Delaware River Basin (DRB) (Fig. 1). The Susquehanna River drains 71 228 km² to the northern end of Chesapeake Bay, contributing ∼ 50 % of freshwater inflow to the estuary (Leathers et al., 2008). Chesapeake Bay is the largest estuary in the United States with a surface area of 11 601 km² and a shoreline extending over 7000 km. Chesapeake Bay has varied tidal characteristics across the estuary, e.g., mixed tide in the northern portion and semidiurnal tide near the bay mouth. In addition to the Susquehanna River, several other large rivers also drain to this estuary. For Chesapeake Bay, the amount of the freshwater outflow is approximately the same as the seawater inflow from the mid-Atlantic coastal waters (Valle-Levinson, 1995). The Delaware River drains 35 070 km² to Delaware Bay and contributes 58 % of freshwater to the estuary (Whitney and Garvine, 2006). Delaware Bay has 2030 km² in surface area and is dominated by semidiurnal tide. The tidal range is 1.5 m at the bay mouth and increases toward Trenton (USGS gauge 01463500 in Fig. 1). Trenton is referred to as the downstream limit of freshwater (Sharp, 1984), as there is a hydraulic jump at 2.7 km downstream of the station due to an abrupt decrease in channel bathymetry (Zhang et al., 2020). Together, SRB and DRB have over 4 million residents and DRB provides drinking water to 6 % of the US population.

Figure 1 An overview of Susquehanna River Basin (SRB) and Delaware River Basin (DRB). This map is created using the free and open source QGIS on the world topographic map (ESRI, 2012).

The locations of in situ observations are provided in Fig. 1. Among over 100 USGS gauges in the mid-Atlantic region, we selected all USGS gauges in the main stem of the Susquehanna River and Delaware River, respectively, for simulated streamflow validation. The water level data at six NOAA tidal gauges are also selected for data analysis and/or model validation. While the coastal tidal gauge (8534720) is used for identifying storm surge, the two tidal gauges near the river mouths (8573364 and 8545240) provide the downstream boundary condition (BC) of the river model for the Susquehanna River and Delaware River, respectively. The four tidal gauges at the downstream section of the Delaware River (8545240, 8546252, 8539094 and 8548989) are used in the model evaluation.

2.2 Global river routing model

The Model for Scale Adaptive River Transport (MOSART) (Li et al., 2013, 2015), the river component in E3SM (Golaz et al., 2019), is used for river modeling. MOSART is a river routing model applicable across local, regional and global scales. The model is driven by runoff from a land surface model and simulates water flow from hillslopes to tributary subnetworks and to main channels. The routing schemes in MOSART include kinematic wave and diffusion wave equations, two simplified forms of the 1-dimensional Saint Venant equations. The routing of surface runoff in hillslopes and tributaries is represented using the kinematic wave method. The flow in the main channel is represented by the diffusive wave method. The momentum equation in the diffusive wave method is (Chow, 1988)

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial h}{\partial x}}-S_{\mathrm{0}}+S_{\mathrm{f}}=\mathrm{0},\end{array}$$

where h is the water depth in the channel, S₀ is the riverbed slope and S_f is the friction slope. Compared to the diffusive wave method, the kinematic wave method neglects the first term of Eq. (1) in its momentum equation. The flow velocity (v) is estimated using the Chezy–Manning equation:

$$\begin{array}{}\text{(2)}& v={\displaystyle \frac{\left|S_{\mathrm{f}}\right|}{S_{\mathrm{f}}}}{n}^{-\mathrm{1}}{R}^{\frac{\mathrm{2}}{\mathrm{3}}}{\left|S_{\mathrm{f}}\right|}^{\frac{\mathrm{1}}{\mathrm{2}}},\end{array}$$

in which Manning's n is used as the frictional coefficient and R is the hydraulic radius. The backwater effect can be represented in the diffusive wave method, as the flow velocity is determined by both the riverbed slope (S₀) and the water level variations (h) along the river channels (Luo et al., 2017). In extreme conditions, when the downstream water stage is higher than that of the current channel, S_f becomes negative, resulting in a backwater. The extreme reverse flow was recently observed in the Mississippi River during Hurricane Ida (Miller, 2021). In our study domain, the backwater processes in the Susquehanna River and Delaware River were reported by USGS (U.S. Geological Survey, 2016) and showed significant impacts during Hurricane Irene (Zhang et al., 2020).

The cross section of the main channel is specified as rectangular in MOSART when the channel water depth is not more than the bankfull depth (H). The channel width (W) and bankfull depth (H) are estimated from the total upstream drainage area (A_total) using empirical formulations (Bent and Waite, 2013):

$$\begin{array}{}\text{(3)}& {\displaystyle }W=a(A_{\mathrm{total}}{)}^b,\text{(4)}& {\displaystyle }H=a(A_{\mathrm{total}}{)}^b,\end{array}$$

where a and b are empirical parameters. When channel water depth exceeds H, an elevation profile is invoked to capture the elevation variation in the floodplain (Luo et al., 2017).

In this study, the runoff inputs for MOSART are obtained from the Global Reach-Level Flood Reanalysis (GRFR) (Yang et al., 2021), an offline simulation from a high-resolution land surface model that has been calibrated and bias corrected. The original configuration of the MOSART diffusive wave method applies a static coastal boundary condition (CBC), i.e., either normal depth or fixed mean sea level at the river mouth (Luo et al., 2017). The normal depth boundary assumes that the friction slope (S_f) is equal to the riverbed slope (S₀) at the river outlet cell. This simplification, while reasonable for global simulations in which the influence of coastal processes is limited, can be problematic in low-lying coastal regions. To represent the backwater effects induced by the dynamic sea level variation, this study introduces a new, dynamic CBC option in MOSART to read in time-varying water level data at a time interval consistent with the land–river coupling frequency in E3SM. The coupling time interval is set as 1 h in this study. This dynamic CBC is only configured for the rivers of interest, while the static CBC is used in all other outlet boundaries due to the limited data availability.

2.3 Coastal refined global unstructured mesh and flow direction map

A global unstructured mesh with a resolution of ∼ 100 km has been developed using the JIGSAW mesh library (Engwirda and Ivers, 2016; Engwirda, 2017), which enables (a) the flexibility of embedding high-resolution subdomains within the ESM's global configuration; (b) oversampled geometrical features (O(< 100 m)), e.g., river network and coastline, to be simplified to coarser ESM length scales on the order of 2–60 km; and (c) close alignment of the complex geography of coastline, watershed boundaries and river networks (Engwirda and Liao, 2021). The global mesh is developed to allow for more seamless coupling of the land, river and ocean components in E3SM (Fig. 2) for more consistent modeling of global surface processes. In the mid-Atlantic region, the mesh resolution is refined to ∼3 km to better resolve local coastal and watershed processes (Fig. 2b). Significant effort is also made to ensure that the cells in the high-resolution mesh match the prescribed dam locations and the orientations of the edges conforming to the flow direction along the main channel.

Figure 2 (a) The global unstructured mesh of E3SM. (b) A magnified view of the mid-Atlantic region (orange square in a). (c) A magnified view near the mouths of the Susquehanna River and Delaware River (orange rectangle in b). The red rectangles in (b) and (c) represent the downstream section of the Delaware River used for backwater propagation analysis.

The river networks and flow directions are modeled using HexWatershed (Liao et al., 2020, 2022; Liao, 2022), a watershed and flow direction model that supports both structured and unstructured meshes for river routing models. HexWatershed uses a topological relationship-based approach to define river networks in the mid-Atlantic region (Lehner et al., 2008). To generate the flow direction for the entire domain, HexWatershed uses a hybrid depression filling and breaching stream burning algorithm to remove local depressions while minimizing modifications to surface elevation and produces flow routing parameters including the flow direction map, channel slope and drainage area, which are critical for accurately representing coupled land–river–ocean processes.

2.4 Extreme event selection

The extreme events of fluvial flood (FF) and storm surge (SS) are separately selected based on their corresponding time-series observations. The selection of FF follows the strategy proposed in (Zhang et al., 2021). In the time-series of discharge data, a flood event is identified by first selecting the flood peaks using the peaks-over-threshold (POT) approach (Lang et al., 1999). The threshold is determined based on automatic threshold selection and the independence of the peak series is examined with a declustering method (Zhang et al., 2021). The start and end dates are specified using the empirical formulation:

$$\begin{array}{}\text{(5)}& {\displaystyle }{Q}_{\mathrm{S}}\le a{Q}_{\mathrm{P}};T_{\mathrm{P}}-T_{\mathrm{S}}\le b\left(\mathrm{5}+\ln \left({\displaystyle \frac{A}{{\mathrm{1.609}}^{\mathrm{2}}}}\right)\right),\text{(6)}& {\displaystyle }{Q}_{\mathrm{E}}\le a{Q}_{\mathrm{P}};T_{\mathrm{E}}-T_{\mathrm{S}}\le b\left(\mathrm{5}+\ln \left({\displaystyle \frac{A}{{\mathrm{1.609}}^{\mathrm{2}}}}\right)\right),\end{array}$$

where T_P, T_S and T_E are the peak date, start date and end date, Q_P, respectively, Q_S and Q_E are the discharge on the corresponding date, and A is the basin drainage area. The empirical parameters a and b are specified as 0.5 and 1.5, respectively.

The selection of SS is performed in three major steps. The first major step is to extract the SS component from the hourly total water level (TWL) data at the NOAA tidal gauge (i.e., 8534720 in Fig. 1): (a) the TWL time series is detrended by extracting the annual mean sea level, which removes the effects of SLR; (b) the predicted astronomical tides are then derived from the harmonic tidal analysis performed on the detrended TWL on a year-by-year basis using eight major tidal constituents; (c) SS is the non-tidal residual obtained by extracting the tides from the detrended TWL.

The second major step is to filter the extreme SS events using a peak detection algorithm (van Brakel, 2014). When a data point is m-fold standard deviations away from the moving mean, an event peak is identified. The start and end times of the corresponding event are the nearest data points that change signs. In addition to m, the other input parameters of this algorithm are lag and influence: lag is the number of observations to smooth the data, or the length of the moving window; influence represents the influence of new signals on the threshold. Here we selected the SS event by setting m= 5, lag = 30 and influence = 0. These parameters are determined to ensure that the selected SS events include all documented hurricane events in the mid-Atlantic region.

The third major step is to select the extreme SS events of interest by extracting the events with SS peaks greater than the 99.5th percentile. When there is an overlap between the FF and SS events, a compound flood event is identified, for which the duration is defined as the combined period of the two events. The applied event selection method can be more accurate than those used in continental and global applications, where the event period is simply determined using a predefined time window (e.g., ±1 d of a peak event) (Nasr et al., 2021; Ward et al., 2018; Wu et al., 2021).

2.5 Decomposition of backwater drivers

The backwater drivers are decomposed into three different levels (Fig. 3). The first level considers the river topology and the forcing that affects the river flow directly. The direct forcing in the backwater zone is considered as the upstream discharge and the TWL at the river mouth. To address the impact of topology, we compare the backwater effects in the Susquehanna River and Delaware River, which differ significantly in terms of the riverbed elevation along the downstream section and the receiving water body. As previous studies have shown the crucial impact of sea level variations on the backwater effects (Yamazaki et al., 2012), we further decompose the TWL into low-frequency surge (LFS) and tide in the second level. The predicted tide is estimated using harmonic tidal analysis and the LFS is obtained by subtracting the tide from the TWL. It should be noted that this level of decomposition must be applied to a tidal gauge since the harmonic tidal analysis requires measurable tidal effects. The third level decomposes the LFS to high discharge, SS and their compound effects. It is worth noting that the LFS differs from the SS in that the LFS is extracted from the TWL at the NOAA gauge nearest to the river mouth (e.g., 8545240 in Fig. 1) where the river discharge can have a significant influence on the LFS. In particular, the LFS is dominated by both river discharge and SS in narrow tidal rivers, such as Delaware Bay. By contrast, the river discharge influence is negligible in larger estuaries or coastal oceans, such as Chesapeake Bay. Previous studies mainly focused on the drivers of river discharge and SS (Nasr et al., 2021; Ward et al., 2018), assuming they are distinct mechanisms that have few mutual interactions. This assumption is valid when the river drains directly to a large receiving water body. However, cases with a river ending within a small estuary have not been explored. In such cases, the high discharge during FF increases the water level within the estuary and attenuates the spatial variation of water stage along the river (Luo et al., 2017), creating backwater effects. Thus, in the third level, we attempt to understand the respective role of high discharge, SS and their interactive impacts on the LFS-induced backwater effects. These drivers are separated based on the selected flood events in Sect. 2.4. Events dominated by the drivers of high discharge, SS and their compound effects are denoted as $\mathrm{LFS}\cap \mathrm{FF}\cap \stackrel{\mathrm{‾}}{\mathrm{SS}}$, $\mathrm{LFS}\cap \mathrm{SS}\cap \stackrel{\mathrm{‾}}{\mathrm{FF}}$ and LFS∩compound, respectively. The symbol ∩ means the intersection period between two events and the overline means to exclude the different drivers if there is an intersection. Each event is measured in terms of the event duration and the peak water level during the event.

Figure 3 Decomposition of the drivers of backwater effects. TWL is total water level, LFS is low-frequency surge, and SS is storm surge.

The trend analysis was performed to assess the annual trends of the backwater drivers. For each event type, we calculated the annual duration, occurrence and the maximum peak value on an annual basis. The impact of SLR is addressed by selecting the SS event without detrending the TWL data, while the other steps are the same as those in Sect. 2.4. The corresponding case is referred to as SS + SLR. The nonparametric Mann–Kendall (MK) test (Tosunoglu and Kisi, 2017) was employed to statistically assess whether there is a monotonic trend. The null hypothesis (H₀) and the alternative (H_a) are no monotonic trend is present and monotonic trend is present, respectively. For each MK test, we set the significance level at 0.1 and calculated the standard MK statistics (Z) and p value. The positive (or negative) value of Z corresponds to the increase (or decrease) trend.

2.6 Numerical experiments

MOSART simulations were performed from 1990 to 2019 with the first year excluded from analysis as the spin-up time. This period has sufficient data coverage of both runoff and the water level at NOAA gauges. We configured five simulations based on the aforementioned downstream CBCs: (a) normal depth, (b) total water level (TWL), (c) mean sea level (MSL), (d) low-frequency surge (LFS) and (e) tide. The backwater effects are quantified by comparing the TWL and MSL simulations in terms of two quantification metrics along the main channel: water depth change Δh and water volume change ΔV:

$$\begin{array}{}\text{(7)}& {\displaystyle }\mathrm{\Delta }h\left(t,i\right)=h_{\exp }\left(t,i\right)-h_{\mathrm{MSL}}(t,i),\text{(8)}& {\displaystyle }\mathrm{\Delta }V\left(t\right)={\sum }_i^N\left(h_{\exp }\left(t,i\right)-h_{\mathrm{MSL}}\left(t,i\right)\right)L\left(i\right)W\left(i\right),\end{array}$$

where the subscript exp represents TWL, LFS or tide, h is the main channel water depth, t is the model output time step and i is the grid cell index, L and W are the length and width of the main channel within the ith cell, respectively, and N is the number of cells with nonzero Δh. The two metrics measure the backwater-induced changes within the river channel that are created from the variation in the downstream water level. The model performance is assessed using the coefficient of determination (r²), RMSE and Kling–Gupta efficiency (KGE): (Gupta et al., 2009)

$$\begin{array}{}\text{(9)}& \mathrm{KGE}=\mathrm{1}-\sqrt{{\left(r-\mathrm{1}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{{\mathit{\delta }}_{\mathrm{m}}}{{\mathit{\delta }}_{\mathrm{o}}}}-\mathrm{1}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\sum X_{\mathrm{m}}}{\sum X_{\mathrm{o}}}}-\mathrm{1}\right)}^{\mathrm{2}}},\end{array}$$

where X_m and X_o represent the model simulation and the observation, respectively, δ_m and δ_o are the corresponding standard deviations, and r is their linear correlation.

2.7 Extreme events

The flood events in the Delaware River during Hurricane Irene (followed by Tropical Storm Lee) and Hurricane Sandy are used to demonstrate the selection of extreme events and the quantification of the backwater effects. Hurricane Irene, one of the most destructive tropical cyclones in US history, made its first landfall on the coast of North Carolina on 27 August 2011 as a Category 1 hurricane, followed by another landfall in southeastern New Jersey on 27 August and a third landfall in New York City. The SS peaked at 1.8 m along the coast of New Jersey and the wind speed was up to 105 km h⁻¹. The maximum rainfall is 10 in. in DRB. Tropical Storm Lee is the subsequent storm event that formed over the Gulf Coast and swept the East Coast. Lee brought 10–12 in. of precipitation to the mid-Atlantic region, resulting in mainly fluvial processes rather than coastal surges (Ye et al., 2020). The combined events caused two consecutive flow peaks in the Delaware River from 28 August to 10 September 2011 (Fig. 4). Irene is affected by the interactive SS and precipitation-induced FF and has been studied extensively as an example of compound flood (Xiao et al., 2021; Zhang et al., 2020). Hurricane Sandy, the largest Atlantic hurricane on US record, made landfall on the coast of New Jersey on 29 October 2012 with a sustained wind speed of 130 km h⁻¹. Sandy caused a maximum SS of 4 m near New York City and 1.5 m near the Delaware coast. The observed peak flow during Sandy was ∼ 800 m³ s⁻¹ at Trenton.

Figure 4 The river discharge evaluation in SRB and DRB: (a) r², (b) KGE, (c) hydrograph at USGS gauge 01463500, (d) hydrograph at USGS gauge 01578310. The triangles and the circles in (a) and (b) represent the USGS gauges in the main stem of the Susquehanna River and Delaware River, respectively.

3 Model evaluation

In this section, the performance of MOSART in simulating river discharge and water level is compared between the experiments using the downstream CBC of normal depth, TWL and MSL to demonstrate the importance of imposing an appropriate downstream CBC. The model performance in reproducing the observed water level in the downstream section of the Delaware River is significantly improved when the TWL is enforced at the boundary.

3.1 River discharge

The MOSART simulated daily discharge is compared with the USGS observations over the simulation period (1991–2019) at the gauges along the main stem of the Susquehanna River and Delaware River (Fig. 4). The coefficient of determination r² and KGE are calculated for each gauge. The MOSART simulation compares reasonably well with r² and KGE, with both being over 0.5 across all gauges (Towner et al., 2019). The model performance is generally higher in the Susquehanna River and decreases toward the upstream regions. The highest r² and KGE (≥ 0.75) are found at the downstream gauges of the two rivers, as the forcing may not capture the runoff accurately for smaller drainage areas. A closer look at the time series of discharge at these two gauges from 2011 to 2012 shows that the model can capture the hydrograph and smaller peaks of river discharge well, and there is no significant difference in model performance among the different CBC configurations because even the USGS gauges closest to the river outlets are still too far upstream to capture effects of dynamic CBCs. The model has a large bias in the Susquehanna River during Hurricane Irene (21–30 August 2011) with the extreme flood peak significantly underestimated. While the observed peak flow is about 20 000 m³ s⁻¹, the simulated flow is about 5000 m³ s⁻¹. This bias is likely caused by uncertainty in the GRFR runoff forcing (Yang et al., 2021). The global runoff forcing, despite being bias corrected, may underestimate the extremes over specific regions. This is a known challenge in global runoff generation schemes. Additionally, there also exists uncertainty in USGS discharge data estimated from rating curves during extreme events (Di Baldassarre et al., 2012), because these measurements occur less frequently and usually do not cover extreme events (Turnipseed and Sauer, 2010). Overall, the evaluation indicates that the MOSART simulated river discharge captures reasonably the spatial and temporal variability of the observed discharge in the mid-Atlantic region.

3.2 Water level

The simulated water level is compared with the observations at four NOAA tidal gauges along the downstream section of the Delaware River. The model performance is quantified in terms of r² and RMSE (Fig. 5). The TWL simulation results in the best performance, in which r² is over 0.5 among all the gauges, much higher than the other two configurations. The lowest RMSE is also obtained from the TWL simulation. By setting the TWL as the downstream CBC, the model's capability to reproduce the water level variation is greatly enhanced. The same conclusion can also be drawn in the time series comparison from 2011 to 2012 (Fig. 6). The TWL simulation accurately captures the small variations in the observed water level, which are missing in the simulations with normal depth and MSL boundary conditions. The extreme peaks are overestimated in the TWL simulation, as well as the normal depth simulation in which no data are enforced at the downstream boundary. The MSL simulation tends to produce smaller variations and lower peaks as the downstream boundary is forced by a constant water level. The overestimation in water level peaks by the TWL and normal depth simulations is most likely a result of the uncertainties in the channel topology in MOSART. Moreover, the diffusive wave equation (Eq. 1) simplifies the momentum transport by neglecting the inertia terms (local and convective). In the diffusive wave method, because the flood wave is considered as subcritical and diffusive (Trigg et al., 2009), the water level is mainly controlled by the upstream discharge. In low-lying rivers, while gravity and friction may not be the dominant forcing, the inertial force related to velocity changes in space and time dominates the flow momentum. As such, the flood wave propagation from the downstream boundary is underestimated in the backwater zone. As shown in Fig. 5, the improvement of the TWL simulation in predicting water level is reduced toward upstream. This is not unexpected in a reduced-physics river model (Hodges, 2013) because as implied in the diffusive wave equation (Eq. 1), the energy head at the downstream CBC is created by the pressure gradient as a result of the variation in water surface but it is lost rapidly upstream due to the increase in riverbed elevation.

Figure 5 The r² and RMSE values of the MOSART simulations at the four NOAA gauges.

Figure 6 Comparison of daily-averaged simulated and observed water level at NOAA gauge 8546252.

The model evaluation results illustrate that a river model on a regionally-refined global mesh can represent backwater effects at the basin scale when properly specified downstream CBC is used. Thus, the model can be used to further examine the contribution of the backwater drivers. It should be noted that the 1D river models are by no means comparable to 3D hydrodynamic models (e.g. models in Gori et al., 2020 and Zhang et al., 2020) at reproducing coastal flood events or resolving the complex flow dynamics (Neal et al., 2012). Therefore, our analysis focuses on a larger temporal scale by extracting the extreme events from a long period, which are then used to quantify the backwater drivers.

4 Flood event simulation

The event selection method (Sect. 2.4) is used to select the extreme SS, LFS and FF events during Hurricane Irene, Tropical Storm Lee and Hurricane Sandy from long-term observations at the NOAA coastal gauge (8534720), the NOAA gauge closest to the river mouth (8545240) and the USGS gauge at Trenton (01463500), respectively (Fig. 7). Gauge 8545240, despite being located at the upstream reach of Delaware Bay, is dominated by semidiurnal tides even at high-flow conditions during extreme storm events (Xiao et al., 2021). A lag time of 4 h is added to the water level data at the coastal gauge to compensate for the phase lag between this NOAA tidal gauge and that at the river mouth. As there is an overlap between the SS and FF events during Irene (Fig. 7a), a compound flood is identified over the combined period. The obvious difference is observed between the LFS and SS events, which were obtained using the same method but at different locations. At the Delaware River mouth, the LFS event can be attributed to both SS and river discharge, resulting in a duration period comparable to an FF event and much longer than an SS event. This highlights the importance of considering the effect of high river discharge on compound flooding, particularly for rivers contributing to a small receiving water body. Sandy did not induce significant fluvial processes, and the LFS event was primarily caused by SS. Thus, no compound flood is identified over this period.

Figure 7 The top panels are the three types of selected extreme events (SS, LFS, FF) overlaid on the corresponding time series of SS, LFS and discharge, represented by solid, dashed and dotted gray lines, respectively, for Hurricane Irene and Tropical Storm Lee (a) and Hurricane Sandy (b). The compound flood event is marked between the two vertical black lines. The bottom left panels are the riverbed elevation along the backwater propagation extent in the Delaware River and the bottom right panels are Δh (color shading) along the upstream distance and ΔV (black curve) derived from the MOSART simulations.

Because the TWL simulation shows a reasonable performance at reproducing both river discharge and water level during the hurricane periods (Sect. 3), the TWL simulation is used to estimate the water depth change (Δh) and water volume change (ΔV) in Eqs. (7) and (8) to assess the backwater effects. Δh shows the backwater propagation extent, which is roughly 60 km upstream from the river mouth (Fig. 7). This extent is determined by the riverbed elevation. Increasing the elevation to ∼ 5 m, Δh implies a large increase in the downstream water level. As a spatially aggregated quantity Δh, ΔV is also consistent with LFS, with ΔV following the LFS variation and peaking on the same dates. By comparing Δh and ΔV over the two hurricanes, it is not difficult to conclude that the compound flood caused by two consecutive events of Hurricane Irene and Tropical Storm Lee has a much larger impact on the backwater effects than that of Hurricane Sandy, even though their LFS peaks are at similar levels. This is probably because the duration of the SS event during Sandy is much shorter than the combined SS and FF events during Irene.

5 Backwater drivers

The interactive effects of water level and discharge inspire us to further decompose the backwater drivers. This section provides the analysis of the drivers based on the three-level decomposition introduced in Sect. 2.5. We examine the contribution of each driver to the backwater effects and their corresponding long-term trend under climate change.

5.1 Decomposition of backwater drivers

5.1.1 Discharge, TWL and topology

The first decomposition level assesses the impacts of river discharge (Q), TWL and topology by comparing ΔV between the Susquehanna River and Delaware River over the entire simulation period (Fig. 8). While ΔV is computed from the TWL simulation, Q and TWL are obtained from the paired streamflow and tidal gauges nearest to the river outlets, i.e., gauge 01578310 and 8573364 for the Susquehanna River and gauge 01463500 and 8545240 for the Delaware River. The result shows the key role of river topology and TWL in affecting backwaters. The maximum ΔV in the Susquehanna River is roughly 3 orders of magnitude smaller than that in the Delaware River. This is the result of a larger gradient in the riverbed elevation profile of the Susquehanna River that impedes the backwater propagation. Over the 5 km downstream section, the elevation increases from 0 to ∼ 20 m in the Susquehanna River but by less than 1 m in the Delaware River.

Figure 8 Scatter plot of ΔV against TWL in the Susquehanna River (a) and Delaware River (b). Colored circles represent the corresponding river discharge (Q). The TWL and Q data are obtained from gauges 8573364 and 01578310, respectively, for the Susquehanna River, and from gauges 8546252 and 01463500 for the Delaware River. The solid lines are the fitted linear regression under different discharge ranges.

Figure 9 The correlation coefficient (r) matrix of ΔV, Q and TWL in the Susquehanna River (a) and Delaware River (b).

In both rivers, ΔV is dominated by TWL, with the corresponding correlation coefficient (r) over 0.9 (Fig. 9). However, the influence of Q differs significantly between the two rivers. In the Susquehanna River, Q is negatively correlated with ΔV ($r=-\mathrm{0.11}$). The increase in Q reduces the TWL impact on ΔV. For instance, at the same TWL, a smaller Q could result in a higher ΔV (Fig. 8a). This behavior is also evident in the slopes of the fitted linear regression lines: the fitted slope for Q≥ 1000 m³ s⁻¹ is smaller than those for low discharge (Fig. 8a). This result is expected because high upstream discharge can attenuate the propagation of downstream backwaters. In the Delaware River, ΔV increases with Q, and r between ΔV and Q is 0.36. The regression slopes are similar at different discharge conditions. These contrasting results between the two rivers imply that the impact of Q on the backwater effects depends on the river's receiving water body. Because the Delaware River contributes to Delaware Bay, a much narrower estuary than Chesapeake Bay, the effect of its discharge on the water level variation of the estuary is much stronger than that of the Susquehanna River. Consistently, there is a much higher r between Q and the TWL (0.27) in the Delaware River than that (−0.03) in the Susquehanna River (Fig. 9). In addition, the channel constriction in the Delaware River might also facilitate the formation of backwaters (Castelltort et al., 2020).

Between the two draining estuaries, Chesapeake Bay behaves more like an ocean as the river impact is limited and the coastal and fluvial processes are distinct. This situation is usually taken as the general case for compound effects and has been addressed extensively using statistical models (Nasr et al., 2021; Ward et al., 2018) and large-scale river models (Ikeuchi et al., 2017). By contrast, cases like the Delaware River have rarely been documented in any global-scale studies, even though they are ubiquitous and may witness more coastal backwaters. A possible reason for why such cases were overlooked is that previous global meshes have a significant deficiency in resolving narrow estuaries properly. Thus, we focus the following analysis on the Delaware River.

5.1.2 LFS and tide

Given the key impact of TWL, the second decomposition level examines the respective role of the LFS and tide. The simulations configured with the downstream BCs of the TWL, LFS and tide (Sect. 2.6) are used to derive the maximum water depth change (Δh_max) for each grid cell:

$$\begin{array}{}\text{(10)}& \mathrm{\Delta }{h}_{\max }\left(i\right)=\max \mathit{\{}\mathrm{\Delta }h\left(t,i\right):t=\mathrm{1},\mathrm{2},\mathrm{\dots },T\mathit{\}},\end{array}$$

where T is the simulation period.

The Δh_max comparison shows the dominance of LFS over tide in increasing the maximum water depth (Fig. 10). In the along-channel profile (Fig. 10a), the value of Δh_max in the LFS simulation is close to that in the TWL simulation: ≥ 1 m within the 40 km upstream range of the outlet and then gradually reduced to 0 with a sharp increase in the riverbed elevation. By contrast, the tide simulation produces a much smaller Δh_max with the value never exceeding 0.5 m. Among the simulation cases, the highest Δh_max occurs at roughly 25 km upstream from the mouth. This along-channel profile reveals the interaction of the discharge and the upstream propagation of tide and surge momentum. It is also observed that Δh_max is slightly higher in the LFS simulation than in the TWL simulation, which is likely the result of the negative impact of low tide on TWL. The spatial variation of Δh_max is shown in Fig. 10b. The backwater effects are limited to the low-lying section of the Delaware River, i.e., downstream of Trenton. Backwater propagation occurs along the main channel as well as some small contributing tributaries, for which the extent is determined by the corresponding elevation. The propagation extent is similar between the TWL and LFS simulations and is much smaller in the tide simulation.

Figure 10 (a) The along-channel profiles of the maximum Δh obtained using the TWL, SS and tide configurations. The river outlet is at x= 0 km. (b) The corresponding spatial map of the maximum Δh over a downstream region of Trenton. This region is specified as the red rectangles in Fig. 2. The black dots represent the Delaware River outlet.

5.1.3 High discharge, SS and compound effect

The LFS impact is further decomposed into high discharge, SS and their compound effect using the LFS simulation. We compared the variation of the event accumulated ΔV with respect to the event duration and peak water level among the drivers (Fig. 11). Regardless of the drivers, ΔV is mainly determined by the event duration. Its value increases linearly with the duration with high correlations. The peak water level provides the secondary effect. Higher peaks generally increase ΔV, resulting in values above the fitted regression line.

Figure 11 Scatter plot of the accumulated ΔV against the event duration in days for the LFS events (a), the SS events (b), the FF events (c) and the compound flood events (d). The color represents the corresponding peak LFS.

Our results indicate that the influence of each driver on ΔV is more dependent on the frequency and duration of the corresponding events rather than their extremes. For example, the FF events are more influential on ΔV because they are more frequent than the SS and compound flood events. Also, the ΔV is higher in the FF and compound flood events than in the SS events because the latter lasts much shorter. A remarkable case occurred during Hurricane Irene when the highest ΔV in our study period was produced by the combined long-lived FF and compound flood events (Fig. 7). Notably, the slope and correlation between ΔV and duration is the largest in the compound events, which means that the strength of the compound events increases more rapidly with duration. In all, our driver comparison indicates that high discharge is the key driver of the backwaters in the Delaware River due to the higher frequency and longer duration of the corresponding FF events.

5.2 Trend analysis

This section presents the analysis of annual trend in the LFS and decomposed drivers, as well as their influences on the backwater effects. The results reveal the impacts of SLR and increasing frequency of SS events during 1990–2019 on exacerbating the backwater effects.

5.2.1 Trend in the backwater drivers

The annual trend in the backwater drivers shows an increasing trend of SS due to both SLR and increasing SS frequency (Fig. 12), with p values of 0.039 and 0.031, respectively, for annual duration and occurrence in the SS + SLR case (Table 1). When SLR is considered, the number of SS occurrences increases from ∼ 1 to 2–5 times per year over the study period (Fig. 12b). Accordingly, the annual duration of the SS event increases from ≤ 15 d (1990–2005) to > 25 d over multiple years (2005–2019) and reaches over 50 d in 2016 (Fig. 12a). The SS peaks are also increased by SLR but do not present any trend (Fig. 12c). We do not notice any clear trends from the FF events or the compound events (Table 1). The annual characteristics of the FF events vary significantly between wet and dry years. In the very wet year of 1996, the duration of the FF events reaches 100 d with up to six events per year. The frequency of compound flood events is low, only occurring from one to three times per year between 2003 and 2012.

Figure 12 The interannual variability of (a) the event annual duration, (b) the event occurrence and (c) the maximum peak value for different types of flood events.

Table 1 The MK statistics of long-term trends in annual duration, occurrence, peak value and ΔV.

Affected by the drivers of high discharge and SS, the LFS events show an increasing trend in both duration and frequency with p values of 0.094 and 0.042, respectively. No clear trend is found in the LFS peaks (Table 1). The LFS trend is basically consistent with that of the SS events with a few exceptions in flood years (i.e., 1996 and 1998) when the LFS events were caused by high discharge. Except for the flood years, the LFS occurrence is as low as one time per year in the 1990s and increases to two to six times per year from 2003. Correspondingly, the LFS duration increases from ≤ 5 to 20–40 d from 2004. The increased trend in LFS is most likely a result of SLR that leads to more frequent occurrences of the SS events.

Figure 13 Similar to Fig. 12 but for the interannual variability of ΔV.

5.2.2 Trend in the backwater effects

The annual trend of the backwater effects is analyzed for the different drivers (Table 1) in terms of the annually accumulated ΔV (Fig. 13). The trend of ΔV is consistent with the event duration and frequency trends of the corresponding drivers. The SS- and LFS-induced backwater effects increase but no clear trends can be observed for the FF and compound flood-induced backwater effects. The ΔV trends are significant in the SS + SLR, SS and LFS cases and insignificant in the FF and compound cases (Table 1). Our results also demonstrate the critical impact from high discharge. The resulting ΔV by FF over the flood years (e.g., 1996) can be several times higher than the ΔV caused by the other drivers. This is probably due to the long duration of the FF events. Because SLR and intensified SS increased coastal backwaters in river channels, our analyses call for better representations of the related processes in ESMs for predictive understanding of associated flood risks under climate change and effects on the water and biogeochemical cycles through land–river–ocean interactions and possible impacts on atmospheric processes. However, we caution against attributing changes based on modeling and analysis of a 30-year period to climate change because internal climate variability and other anthropogenic effects most likely also play important roles in the increasing sea level and SS frequency.

6 Discussion

Our study shows that by using the diffusive wave method, large-scale river models configured on a coastal refined mesh are capable of reproducing backwater effects in low-lying river channels with appropriate downstream boundaries. The global unstructured mesh alleviates the computational burden in ESMs by relaxing the resolution in the inland and offshore regions to ∼ 100 km while embedding regions of “high” resolution of O(1 km) near the river–ocean interface. Although the resolution is still not comparable to that used in local-scale models, the mesh is able to resolve the complex river networks near the coastline without having to merge multiple outlets into a single cell. The downstream BC is critical for connecting the coastal and fluvial processes, transferring the water head energy upstream and thus simulating the backwater effects. The widely used diffusive wave method that uses the normal depth boundary and the more simplified kinematic wave method may be only applicable in high-gradient regions, as these methods do not incorporate any downstream information.

As an important finding, this study revealed the crucial difference in flood drivers between two distinct coastal rivers (i.e., Susquehanna vs. Delaware), with the former connected to a wide ocean-like estuary and the latter connected to a narrow river-like estuary, which is usually ignored by ESMs. The difference is mainly caused by the effects of river geometry and estuary size on LFS, which is a direct driver of backwaters. For an ocean-like estuary, such as Chesapeake Bay, river discharge from its drainage basins hardly affects the water level fluctuations of the estuary. But when a coastal river drains to a narrow estuary, its LFS would be driven by not only SS but also river discharge and their compounding. Further, the backwater effects created by high discharge and SS are different. While SS generates an upstream-propagated energy head, high discharge gradually builds up the water level of the receiving water body. The increased water level would slowly move upstream, attenuating the river stage fluctuation and flood waves. High discharge that presents a higher frequency and a longer duration can occur in close successions with SS during compound flooding, e.g., Hurricane Irene and Tropical Storm Lee, creating extended backwater effects. In all, we show that in addition to the conventional flood drivers, such as riverbed elevation and sea level, ESMs need to properly represent the flood drivers for small estuaries, such as river discharge.

We demonstrate that the backwater effects are significant in the low-lying watersheds and have an increased trend over the past 30 years. In the Delaware River, the propagated backwaters account for up to 1.2 m increase in water depth and 1×10⁸ m³ increase in water volume per day during an extreme event. The effects could be several orders higher for larger river basins. This increased flood risk will otherwise be underestimated if the backwater effects are not properly represented in ESMs. Furthermore, our simulation also shows the increased influence of climate change on backwaters, with SLR and more frequent SS increasing the strength of backwaters in the mid-Atlantic region.

Notably, there are still a few limitations in the river model used in this study, which may introduce uncertainties to our simulations. First, as a large-scale river transport model, MOSART simplifies the channel cross section as rectangular and trapezoidal when the water depth is below and above the river's bankfull depth, respectively, partially due to lack of large-scale river cross-section data (Li et al., 2013). This simplification may affect the accuracy of simulated water depth for rivers with very irregular channel cross sections. The river channel width and bankfull depth estimated from empirical formulations may introduce uncertainties. Although such estimation achieves reasonable accuracy at local basins, more reliable river geometry data should be considered at least for regions wherever the data are available. While global river width datasets have been developed for rivers with width > 90 m (Allen and Pavelsky, 2018; Yamazaki et al., 2014), the river bankfull depth may also be derived from high-resolution remote sensing data. However, it remains challenging to upscale the observed river geometry to model resolution given the river is resolution free (Liao et al., 2022). Second, the global runoff forcing may underestimate the event extremes, such as the discharge peak during Hurricane Irene, affecting the reliability of backwater quantifications during the corresponding event. As we aim at generalizing the analyses to the global scale, the bias-corrected global forcing that can capture extremes is desired. Third, we demonstrate that the direct impact of tide on the backwater effects is limited compared to LFS, but the quantification assumes that LFS and tide can be separated by removing tide from the total water level, ignoring the nonlinear interaction between tide and surge. In reality, as a dynamic component in Delaware Bay, the interaction between high tide and coastal surge may stimulate a further increase in the water level (Krien et al., 2017). Last but not least, because our simulation takes the in situ observation as the boundary data, it does not account for the interaction between fluvial and coastal processes. The backwater effects are constrained by the prescribed boundary sea level. However, the river–ocean interface is a dynamic region that features complex multiscale processes. In the context of an intense storm, the upstream propagation of SS impedes river discharge, which in turn modulates the water level (Dykstra and Dzwonkowski, 2020). It remains unclear whether this mutual interaction is critical in flood modeling and how it responds to SLR (Kulp and Strauss, 2019) and enhanced tidal dynamics (Talke and Jay, 2020) due to climate change. Thus, this study provides a basis for modeling coastal induced flooding in river basins. We aim to couple MOSART with the E3SM ocean model interactively in our next step to resolve the complex interactions at the river–ocean interface and eventually couple all land-river-ocean processes within E3SM to improve predictive understanding of compound flooding through pluvial, fluvial and coastal processes (Xu et al., 2022).

7 Conclusion

This research assesses the capability of the global-scale MOSART river transport model to simulate the coastal backwater effects at the basin scale by imposing the observed water level at the downstream boundary and using a coastal refined unstructured mesh. The simulation is evaluated at two major river basins of the densely populated mid-Atlantic region. MOSART shows reasonable agreement with the observed river discharge in both rivers and captures the water level variations in the downstream section of the Delaware River, indicating the model's capability to represent the backwater effects. We performed numerical experiments and extracted extreme flood events to examine the contribution of various backwater drivers. Our analyses revealed the dependence of the backwater drivers on the river geometry and the river's receiving water body. While SS is considered a dominant forcing, the impact of high discharge can be significant in a narrow river-like estuary, such as the Delaware River, as the discharge modulates the low-frequency water level variations within the estuary. In addition, high discharge when occurring simultaneously with SS creates strong compound flooding. The extreme impact of a compound event should be mainly attributed to the extended duration of the combined coastal and fluvial processes rather than to extreme flood peaks. The trend analysis of the backwater drivers shows that the strength of backwaters in the Delaware River has been increasing in recent decades due to SLR and more frequent SS. In the future, we plan to extend the current work from the mid-Atlantic region to the global domain and refine the coastal mesh globally. A framework of two-way coupled river and ocean models will be established to understand the complex river–ocean interactions.