Articles | Volume 26, issue 21
Hydrol. Earth Syst. Sci., 26, 5473–5491, 2022
Hydrol. Earth Syst. Sci., 26, 5473–5491, 2022
Research article
03 Nov 2022
Research article | 03 Nov 2022

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

Investigating coastal backwater effects and flooding in the coastal zone using a global river transport model on an unstructured mesh
Dongyu Feng1, Zeli Tan1, Darren Engwirda2, Chang Liao1, Donghui Xu1, Gautam Bisht1, Tian Zhou1, Hong-Yi Li3, and L. Ruby Leung1 Dongyu Feng et al.
  • 1Atmospheric Sciences and Global Change Division, Pacific Northwest National Laboratory, Richland, WA, 99354, USA
  • 2T-3 Fluid Dynamics and Solid Mechanics Group, Los Alamos National Laboratory, Los Alamos, NM, 87545, USA
  • 3Department of Civil and Environmental Engineering, University of Houston, TX, 77204, USA

Correspondence: Zeli Tan (


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 km2 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 km2 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 km2 to Delaware Bay and contributes 58 % of freshwater to the estuary (Whitney and Garvine, 2006). Delaware Bay has 2030 km2 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 1An 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)

(1) h x - S 0 + S f = 0 ,

where h is the water depth in the channel, S0 is the riverbed slope and Sf 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:

(2) v = S f S f n - 1 R 2 3 S f 1 2 ,

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 (S0) 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, Sf 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 (Atotal) using empirical formulations (Bent and Waite, 2013):


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 (Sf) is equal to the riverbed slope (S0) 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:


where TP, TS and TE are the peak date, start date and end date, QP, respectively, QS and QE 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 LFSFFSS, LFSSSFF 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 3Decomposition 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 (H0) and the alternative (Ha) 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:


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 (r2), RMSE and Kling–Gupta efficiency (KGE): (Gupta et al., 2009)

(9) KGE = 1 - r - 1 2 + δ m δ o - 1 2 + X m X o - 1 2 ,

where Xm and Xo 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−1. 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−1. 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 m3 s−1 at Trenton.

Figure 4The river discharge evaluation in SRB and DRB: (a) r2, (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 r2 and KGE are calculated for each gauge. The MOSART simulation compares reasonably well with r2 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 r2 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 m3 s−1, the simulated flow is about 5000 m3 s−1. 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 r2 and RMSE (Fig. 5). The TWL simulation results in the best performance, in which r2 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 5The r2 and RMSE values of the MOSART simulations at the four NOAA gauges.


Figure 6Comparison 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 7The 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 8Scatter 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 9The 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=-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 m3 s−1 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 (Δhmax) for each grid cell:

(10) Δ h max i = max { Δ h t , i : t = 1 , 2 , , T } ,

where T is the simulation period.

The Δhmax 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 Δhmax 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 Δhmax with the value never exceeding 0.5 m. Among the simulation cases, the highest Δhmax 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 Δhmax 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 Δhmax 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 11Scatter 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 12The 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 1The MK statistics of long-term trends in annual duration, occurrence, peak value and ΔV.

Download Print Version | Download XLSX

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 13Similar 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×108 m3 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.

Code and data availability

The water level observation data is available via (NOAA, 2022). The river discharge data is available via (U.S. Geological Survey, 2016). The source code developed in this study is available on Zenodo (, Feng, 2022).

Author contributions

DF and ZT designed the methodology and the numerical experiments. DE, CL and DX prepared the unstructured global mesh, MOSART river networks and model parameters. DF carried out the analysis and the result visualization. DF and ZT wrote the initial draft of the manuscript. All authors contributed to the discussion and review of the results and to the editing of the paper.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


All model simulations of this research were performed using resources available through Research Computing at Pacific Northwest National Laboratory (PNNL). PNNL is operated for DOE by Battelle Memorial Institute, United States, under contract DE-AC05-76RL01830. We thank two anonymous reviewers for their constructive comments.

Financial support

This research has been supported by the Earth System Model Development program areas of the U.S. Department of Energy, Office of Science, Office of Biological and Environmental Research as part of the multi-program, collaborative Integrated Coastal Modeling (ICoM) project (grant no. KP1703110/75415).

Review statement

This paper was edited by Hubert H. G. Savenije and reviewed by two anonymous referees.


Alfieri, L., Feyen, L., and Di Baldassarre, G.: Increasing flood risk under climate change: a pan-European assessment of the benefits of four adaptation strategies, Climatic Change, 136, 507–521,, 2016. 

Allen, G. H. and Pavelsky, T. M.: Global extent of rivers and streams, Science, 361, 585–588,, 2018. 

Bakhtyar, R., Maitaria, K., Velissariou, P., Trimble, B., Mashriqui, H., Moghimi, S., Abdolali, A., Van der Westhuysen, A., Ma, Z., and Clark, E.: A new 1D/2D coupled modeling approach for a riverine-estuarine system under storm events: Application to Delaware River Basin, J. Geophys. Res.-Oceans, 125, e2019JC015822,, 2020. 

Bates, P. D., Quinn, N., Sampson, C., Smith, A., Wing, O., Sosa, J., Savage, J., Olcese, G., Neal, J., and Schumann, G.: Combined modeling of US fluvial, pluvial, and coastal flood hazard under current and future climates, Water Resour. Res., 57, e2020WR028673,, 2021. 

Bent, G. C. and Waite, A. M.: Equations for estimating bankfull channel geometry and discharge for streams in Massachusetts, U.S. Geological Survey Scientific Investigations Report 2013–5155,, 2013. 

Bermúdez, M., Farfán, J., Willems, P., and Cea, L.: Assessing the effects of climate change on compound flooding in coastal river areas, Water Resour. Res., 57, e2020WR029321,, 2021. 

Bilskie, M. and Hagen, S.: Defining flood zone transitions in low-gradient coastal regions, Geophys. Res. Lett., 45, 2761–2770,, 2018. 

Caldwell, P. M., Mametjanov, A., Tang, Q., Van Roekel, L. P., Golaz, J. C., Lin, W., Bader, D. C., Keen, N. D., Feng, Y., and Jacob, R.: The DOE E3SM coupled model version 1: Description and results at high resolution, J. Adv. Model. Earth Sy., 11, 4095–4146,, 2019. 

Camelo, J., Mayo, T. L., and Gutmann, E. D.: Projected Climate Change Impacts on Hurricane Storm Surge Inundation in the Coastal United States, Frontiers in Built Environment, 207, 588049,, 2020. 

Castelltort, F. X., Bladé, E., Balasch, J., and Ribé, M.: The backwater effect as a tool to assess formative long-term flood regimes, Quaternary Int., 538, 29–43,, 2020. 

Chegini, T., de Almeida Coelho, G., Ratcliff, J., Ferreira, C. M., Mandli, K., Burke, P., and Li, H. Y.: A Novel Framework for Parametric Analysis of Coastal Transition Zone Modeling, JAWRA J. Am. Water Resour. As., 58, 86–103,, 2022. 

Chen, W.-B., Liu, W.-C., and Wu, C.-Y.: Coupling of a one-dimensional river routing model and a three-dimensional ocean model to predict overbank flows in a complex river–ocean system, Appl. Math. Model., 37, 6163–6176,, 2013. 

Chow, V., Maidment, D., and Mays, L.: Applied Hydrology, MacGraw-Hill, Inc., New York, ISBN 0-07-100174-3, 1988 

Di Baldassarre, G., Laio, F., and Montanari, A.: Effect of observation errors on the uncertainty of design floods, Phys. Chem. Earth Pt. A/B/C, 42, 85–90,, 2012. 

Dykstra, S. and Dzwonkowski, B.: The propagation of fluvial flood waves through a backwater-estuarine environment, Water Resour. Res., 56, e2019WR025743,, 2020. 

Dykstra, S. and Dzwonkowski, B.: The Role of Intensifying Precipitation on Coastal River Flooding and Compound River-Storm Surge Events, Northeast Gulf of Mexico, Water Resour. Res., 57, e2020WR029363,, 2021. 

Eilander, D., Couasnon, A., Ikeuchi, H., Muis, S., Yamazaki, D., Winsemius, H. C., and Ward, P. J.: The effect of surge on riverine flood hazard and impact in deltas globally, Environ. Res. Lett., 15, 104007,, 2020. 

Engwirda, D.: JIGSAW-GEO (1.0): locally orthogonal staggered unstructured grid generation for general circulation modelling on the sphere, Geosci. Model Dev., 10, 2117–2140,, 2017. 

Engwirda, D. and Ivers, D.: Off-centre Steiner points for Delaunay-refinement on curved surfaces, Comput. Aided Design, 72, 157–171,, 2016. 

Engwirda, D. and Liao, C.: “Unified” Laguerre-Power Meshes for Coupled Earth System Modelling, 29th International Meshing Roundtable (IMR), Virtual Conference, Zenodo,, 2021. 

ESRI: “Topographic” [basemap], Scale Not Given, “World Topographic Map”, (last access: 12 October 2022), 2012. 

Feng, D.: E3SM MOSART data ocean mode, Zenodo [code],, 2022. 

Golaz, J. C., Caldwell, P. M., Van Roekel, L. P., Petersen, M. R., Tang, Q., Wolfe, J. D., Abeshu, G., Anantharaj, V., Asay-Davis, X. S., and Bader, D. C.: The DOE E3SM coupled model version 1: Overview and evaluation at standard resolution, J. Adv. Model. Earth Sy., 11, 2089–2129,, 2019. 

Gori, A., Lin, N., and Smith, J.: Assessing compound flooding from landfalling tropical cyclones on the North Carolina coast, Water Resour. Res., 56, e2019WR026788,, 2020. 

Gupta, H. V., Kling, H., Yilmaz, K. K., and Martinez, G. F.: Decomposition of the mean squared error and NSE performance criteria: Implications for improving hydrological modelling, J. Hydrol., 377, 80–91,, 2009. 

Hellmers, S. and Fröhle, P.: Computation of backwater effects in surface waters of lowland catchments including control structures – an efficient and re-usable method implemented in the hydrological open-source model Kalypso-NA (4.0), Geosci. Model Dev., 15, 1061–1077,, 2022. 

Hendry, A., Haigh, I. D., Nicholls, R. J., Winter, H., Neal, R., Wahl, T., Joly-Laugel, A., and Darby, S. E.: Assessing the characteristics and drivers of compound flooding events around the UK coast, Hydrol. Earth Syst. Sci., 23, 3117–3139,, 2019. 

Hinkel, J., Lincke, D., Vafeidis, A. T., Perrette, M., Nicholls, R. J., Tol, R. S., Marzeion, B., Fettweis, X., Ionescu, C., and Levermann, A.: Coastal flood damage and adaptation costs under 21st century sea-level rise, P. Natl. Acad. Sci. USA, 111, 3292–3297,, 2014. 

Hirabayashi, Y., Mahendran, R., Koirala, S., Konoshima, L., Yamazaki, D., Watanabe, S., Kim, H., and Kanae, S.: Global flood risk under climate change, Nat. Clim. Change, 3, 816–821,, 2013. 

Hodges, B. R.: Challenges in continental river dynamics, Environ. Modell. Softw., 50, 16–20,, 2013. 

Hwang, S., Son, S., Lee, C., and Yoon, H.-D.: Quantitative assessment of inundation risks from physical contributors associated with future storm surges: a case study of Typhoon Maemi (2003), Nat. Hazards, 104, 1389–1411,, 2020. 

Ikeuchi, H., Hirabayashi, Y., Yamazaki, D., Muis, S., Ward, P. J., Winsemius, H. C., Verlaan, M., and Kanae, S.: Compound simulation of fluvial floods and storm surges in a global coupled river-coast flood model: Model development and its application to 2007 Cyclone Sidr in Bangladesh, J. Adv. Model. Earth Sy., 9, 1847–1862,, 2017. 

Krien, Y., Testut, L., Islam, A., Bertin, X., Durand, F., Mayet, C., Tazkia, A., Becker, M., Calmant, S., and Papa, F.: Towards improved storm surge models in the northern Bay of Bengal, Cont. Shelf Res., 135, 58–73,, 2017. 

Kulp, S. A. and Strauss, B. H.: New elevation data triple estimates of global vulnerability to sea-level rise and coastal flooding, Nat. Commun., 10, 1–12,, 2019. 

Lamb, M. P., Nittrouer, J. A., Mohrig, D., and Shaw, J.: Backwater and river plume controls on scour upstream of river mouths: Implications for fluvio-deltaic morphodynamics, J. Geophys. Res.-Earth, 117, F01002,, 2012. 

Lang, M., Ouarda, T. B., and Bobée, B.: Towards operational guidelines for over-threshold modeling, J. Hydrol., 225, 103–117,, 1999. 

Leathers, D. J., Malin, M. L., Kluver, D. B., Henderson, G. R., and Bogart, T. A.: Hydroclimatic variability across the Susquehanna River Basin, USA, since the 17th century, Int. J. Climatol., 28, 1615–1626,, 2008. 

Lehner, B., Verdin, K., and Jarvis, A.: New global hydrography derived from spaceborne elevation data, EOS T. Am. Geophys. Un., 89, 93–94,, 2008. 

Leung, L. R., Bader, D. C., Taylor, M. A., and McCoy, R. B.: An introduction to the E3SM special collection: Goals, science drivers, development, and analysis, J. Adv. Model. Earth Sy., 12, e2019MS001821,, 2020. 

Li, H., Wigmosta, M. S., Wu, H., Huang, M., Ke, Y., Coleman, A. M., and Leung, L. R.: A physically based runoff routing model for land surface and earth system models, J. Hydrometeorology, 14, 808–828,, 2013. 

Li, H.-Y., Leung, L. R., Getirana, A., Huang, M., Wu, H., Xu, Y., Guo, J., and Voisin, N.: Evaluating global streamflow simulations by a physically based routing model coupled with the community land model, J. Hydrometeorol., 16, 948–971,, 2015. 

Li, H.-Y., Tan, Z., Ma, H., Zhu, Z., Abeshu, G. W., Zhu, S., Cohen, S., Zhou, T., Xu, D., and Leung, L. R.: A new large-scale suspended sediment model and its application over the United States, Hydrol. Earth Syst. Sci., 26, 665–688,, 2022. 

Liao, C.: HexWatershed: a mesh independent flow direction model for hydrologic models (0.1.1), Zenodo [code],, 2022. 

Liao, C., Tesfa, T., Duan, Z., and Leung, L. R.: Watershed delineation on a hexagonal mesh grid, Environ. Modell. Softw., 128, 104702,, 2020. 

Liao, C., Zhou, T., Xu, D., Barnes, R., Bisht, G., Li, H.-Y., Tan, Z., Tesfa, T., Duan, Z., and Engwirda, D.: Advances in hexagon mesh-based flow direction modeling, Adv. Water Resour., 160, 104099,, 2022. 

Luo, X., Li, H.-Y., Leung, L. R., Tesfa, T. K., Getirana, A., Papa, F., and Hess, L. L.: Modeling surface water dynamics in the Amazon Basin using MOSART-Inundation v1.0: impacts of geomorphological parameters and river flow representation, Geosci. Model Dev., 10, 1233–1259,, 2017. 

Mao, Y., Zhou, T., Leung, L. R., Tesfa, T. K., Li, H. Y., Wang, K., Tan, Z., and Getirana, A.: Flood inundation generation mechanisms and their changes in 1953–2004 in global major river basins, J. Geophys. Res.-Atmos., 124, 11672–11692,, 2019. 

Marsooli, R., Lin, N., Emanuel, K., and Feng, K.: Climate change exacerbates hurricane flood hazards along US Atlantic and Gulf Coasts in spatially varying patterns, Nat. Commun., 10, 1–9,, 2019. 

Mikhailov, V. and Gorin, S.: New definitions, regionalization, and typification of river mouth areas and estuaries as their parts, Water Resources, 39, 247–260,, 2012. 

Miller, B.: Hurricane ida forces mississippi river to reverse flow, (last access: 12 October 2022), 2021. 

Muñoz, D., Moftakhari, H., and Moradkhani, H.: Compound effects of flood drivers and wetland elevation correction on coastal flood hazard assessment, Water Resour. Res., 56, e2020WR027544,, 2020. 

Nasr, A. A., Wahl, T., Rashid, M. M., Camus, P., and Haigh, I. D.: Assessing the dependence structure between oceanographic, fluvial, and pluvial flooding drivers along the United States coastline, Hydrol. Earth Syst. Sci., 25, 6203–6222,, 2021. 

Neal, J., Villanueva, I., Wright, N., Willis, T., Fewtrell, T., and Bates, P.: How much physical complexity is needed to model flood inundation?, Hydrol. Process., 26, 2264–2282,, 2012. 

NOAA: NOAA tides & currents,, last access: 28 October 2022. 

Orton, P., Lin, N., Gornitz, V., Colle, B., Booth, J., Feng, K., Buchanan, M., Oppenheimer, M., and Patrick, L.: New York City panel on climate change 2019 report chapter 4: coastal flooding, Ann. NY Acad. Sci., 1439, 95–114,, 2019. 

Paiva, R. C., Collischonn, W., and Buarque, D. C.: Validation of a full hydrodynamic model for large-scale hydrologic modelling in the Amazon, Hydrol. Process., 27, 333–346,, 2013. 

Rahmstorf, S.: Rising hazard of storm-surge flooding, P. Natl. Acad. Sci. USA, 114, 11806–11808,, 2017. 

Rayson, M. D., Gross, E. S., and Fringer, O. B.: Modeling the tidal and sub-tidal hydrodynamics in a shallow, micro-tidal estuary, Ocean Model., 89, 29–44,, 2015. 

Santiago-Collazo, F. L., Bilskie, M. V., and Hagen, S. C.: A comprehensive review of compound inundation models in low-gradient coastal watersheds, Environ. Modell. Softw., 119, 166–181,, 2019. 

Sharp, J. H. (Ed.): Excerpts from the Delaware Estuary: Research as Background for Estuarine Management and Development : a Report to the Delaware River and Bay Authority, University of Delaware, Sea Grant College Program,, 1984. 

Sun, N., Wigmosta, M. S., Judi, D., Yang, Z., Xiao, Z., and Wang, T.: Climatological analysis of tropical cyclone impacts on hydrological extremes in the Mid-Atlantic region of the United States, Environ. Res. Lett., 16, 124009,, 2021. 

Syvitski, J. P. and Saito, Y.: Morphodynamics of deltas under the influence of humans, Global Planet. Change, 57, 261–282,, 2007. 

Talke, S. A. and Jay, D. A.: Changing tides: The role of natural and anthropogenic factors, Annu. Rev. Mar. Sci, 12, 121–151,, 2020. 

Tellman, B., Sullivan, J., Kuhn, C., Kettner, A., Doyle, C., Brakenridge, G., Erickson, T., and Slayback, D.: Satellite imaging reveals increased proportion of population exposed to floods, Nature, 596, 80–86,, 2021. 

Tosunoglu, F. and Kisi, O.: Trend analysis of maximum hydrologic drought variables using Mann–Kendall and Şen's innovative trend method, River Res. Appl., 33, 597–610,, 2017. 

Towner, J., Cloke, H. L., Zsoter, E., Flamig, Z., Hoch, J. M., Bazo, J., Coughlan de Perez, E., and Stephens, E. M.: Assessing the performance of global hydrological models for capturing peak river flows in the Amazon basin, Hydrol. Earth Syst. Sci., 23, 3057–3080,, 2019. 

Trigg, M. A., Wilson, M. D., Bates, P. D., Horritt, M. S., Alsdorf, D. E., Forsberg, B. R., and Vega, M. C.: Amazon flood wave hydraulics, J. Hydrol., 374, 92–105,, 2009. 

Trigg, M., Birch, C., Neal, J., Bates, P., Smith, A., Sampson, C., Yamazaki, D., Hirabayashi, Y., Pappenberger, F., and Dutra, E.: The credibility challenge for global fluvial flood risk analysis, Environ. Res. Lett., 11, 094014,, 2016. 

Turnipseed, D. P. and Sauer, V. B.: Discharge measurements at gaging stations, US Geological Survey 2328-7055,, 2010. 

U.S. Geological Survey: National Water Information System data available on the World Wide Web (USGS Water Data for the Nation), (last access: 28 October 2022), 2016. 

van Brakel, J. P. G.: Robust peak detection algorithm using z-scores, (last access: 12 October 2022), 2014. 

Valle-Levinson, A.: Observations of barotropic and baroclinic exchanges in the lower Chesapeake Bay, Cont. Shelf Res., 15, 1631–1647,, 1995. 

Ward, P. J., Couasnon, A., Eilander, D., Haigh, I. D., Hendry, A., Muis, S., Veldkamp, T. I., Winsemius, H. C., and Wahl, T.: Dependence between high sea-level and high river discharge increases flood hazard in global deltas and estuaries, Environ. Res. Lett., 13, 084012,, 2018. 

Ward, N. D., Megonigal, J. P., Bond-Lamberty, B., Bailey, V. L., Butman, D., Canuel, E. A., Diefenderfer, H., Ganju, N. K., Goñi, M. A., and Graham, E. B.: Representing the function and sensitivity of coastal interfaces in Earth system models, Nat. Commun., 11, 1–14,, 2020. 

Whitney, M. M. and Garvine, R. W.: Simulating the Delaware Bay buoyant outflow: Comparison with observations, J. Phys. Oceanogr., 36, 3–21,, 2006. 

Wu, W., Westra, S., and Leonard, M.: Estimating the probability of compound floods in estuarine regions, Hydrol. Earth Syst. Sci., 25, 2821–2841,, 2021. 

Xiao, Z., Yang, Z., Wang, T., Sun, N., Wigmosta, M., and Judi, D.: Characterizing the non-linear interactions between tide, storm surge, and river flow in the delaware bay estuary, United States, Front. Mar. Sci., 8, 715557,, 2021. 

Xu, D., Bisht, G., Zhou, T., Leung, L. R., and Pan, M.: Development of Land-River Two-Way Hydrologic Coupling for Floodplain Inundation in the Energy Exascale Earth System Model, J. Adv. Model. Earth Sy., 14, e2021MS002772,, 2022. 

Yamazaki, D., Lee, H., Alsdorf, D. E., Dutra, E., Kim, H., Kanae, S., and Oki, T.: Analysis of the water level dynamics simulated by a global river model: A case study in the Amazon River, Water Resour. Res., 48, W09508,, 2012. 

Yamazaki, D., O'Loughlin, F., Trigg, M. A., Miller, Z. F., Pavelsky, T. M., and Bates, P. D.: Development of the global width database for large rivers, Water Resour. Res., 50, 3467–3480,, 2014. 

Yang, Y., Pan, M., Lin, P., Beck, H. E., Zeng, Z., Yamazaki, D., David, C. H., Lu, H., Yang, K., and Hong, Y.: Global Reach-Level 3-Hourly River Flood Reanalysis (1980–2019), B. Am. Meteorol. Soc., 102, E2086–E2105,, 2021. 

Ye, F., Zhang, Y. J., Yu, H., Sun, W., Moghimi, S., Myers, E., Nunez, K., Zhang, R., Wang, H. V., and Roland, A.: Simulating storm surge and compound flooding events with a creek-to-ocean model: Importance of baroclinic effects, Ocean Model., 145, 101526,, 2020. 

Zhang, Q., Zhang, L., She, D., Wang, S., Wang, G., and Zeng, S.: Automatic procedure for selecting flood events and identifying flood characteristics from daily streamflow data, Environ. Modell. Softw., 145, 105180,, 2021.  

Zhang, Y. J., Ye, F., Yu, H., Sun, W., Moghimi, S., Myers, E., Nunez, K., Zhang, R., Wang, H., and Roland, A.: Simulating compound flooding events in a hurricane, Ocean Dynam., 70, 621–640,, 2020. 

Short summary
Sea level rise, storm surge and river discharge can cause coastal backwater effects in downstream sections of rivers, creating critical flood risks. This study simulates the backwater effects using a large-scale river model on a coastal-refined computational mesh. By decomposing the backwater drivers, we revealed their relative importance and long-term variations. Our analysis highlights the increasing strength of backwater effects due to sea level rise and more frequent storm surge.