To evaluate the applicability of an integrated framework for total water level modeling in upstream tidal areas, we established a dedicated hydrodynamic numerical modeling domain specifically for the Potomac River. The ADCIRC coastal model was utilized for this purpose. The framework underwent calibration to accurately simulate astronomical tides and water level variations driven by riverine discharge, as described in section 2.2. Furthermore, the model's performance was assessed by comparing its results with six historical flooding events recorded at the Washington, DC NOAA station (section 2.3). This validation process ensured that the model accurately captures the simultaneous effects of multiple flood drivers. Subsequently, our focus shifted to quantifying the changes in water levels resulting from individual flood drivers, as well as their co-occurrence, using return period events (section 2.4). Finally, we introduced a ranking system that evaluates the influence of each flood driver on localized water level changes. This ranking system was established through return period event simulations (section 2.5).
2.1. Study Area
This study focuses on Washington, DC, specifically its location at the confluence of the Potomac and Anacostia River. NOAA maintains a water level recording station in Washington, DC (referred to as WASD), which has documented an increase in the frequency of minor-level flooding days from 24 in 2010 to 43 in 2020. Minor-level flooding, as defined by the National Weather Service (NWS), refers to incidents with minimal to no property damage but the potential for some public threat.
The Potomac River, the largest tributary of the Chesapeake Bay, spans approximately 166 km from Washington, DC to the Chesapeake Bay. Figure 1 displays the location of the Potomac River, which experiences tidal fluctuations due to its connection with the Chesapeake Bay downstream. However, the tidal influence diminishes near the Little Falls Pump station (LFMD), which serves as a United States Geological Survey (USGS) river discharge gauge upstream of Chain Bridge.
LFMD has a total drainage area of 29,940 square kilometers (km2) and exhibits an average daily flow of 334 cubic meters per second (m3/s). During heavy rainfall events, the flow can exceed 600 m3/s, as indicated by the NWS's Action Stage. On the eastern side of Washington, DC, the Anacostia River joins the Potomac near the Navy Yard (WNDC). The Anacostia River is shallower and narrower than the Potomac River (McDowell, 2016), and the average daily flow near Bladensburg is significantly smaller (< 4 m3/s) compared to the inflow from LFMD due to its drainage area of 315 km2 (Huanxin et al., 1997).
Moreover, several small urban streams discharge into the Potomac River between Chain Bridge and Fort Washington Park, as indicated by the blue arrows in Fig. 1. The drainage areas of these streams range from 104 to 521 km2. The mean tidal range at the Washington, DC station is approximately 0.9 meters (m), with a tidal phase lagging 5 hours (h) behind Lewisetta (LWTV) and 11.5 h behind Hampton Roads at Sewells Point (SWPV), which is situated at the mouth of the Chesapeake Bay (NOAA tides and currents).
Figure 1 Map of study area. a) Overview map indicates the location of the study area in reference to Chesapeake Bay. Locations of two NOAA recording stations (Lewisetta and Swells Point) used for model validations is shown as black diamonds. b) Computational Model (Red Polygon) domain in Potomac River. Downstream boundary condition is represented as a blue line near Lewisetta NOAA gage (LWTV). Other NOAA gages used for model validations are shown as black diamonds. NOAA tide predictions are shown as light blue squares. Locations of stream flow gages is indicated using orange triangles. c) Shows an insert focusing on Washington, DC. Purple arrows indicate location of river inflows whereas blue arrows show location of urban flows or lateral flows.
Figure 1 illustrates the entire modeling domain of the study area, including three boundary conditions: (1) major upstream river boundaries based on observed discharge data from LFMD on the west, as well as the Northeast and Northwest Anacostia River USGS streamflow gauges (NWMD and NEMD) on the east, shown as purple arrows; (2) small stream flow boundaries based on available USGS observed discharge, represented by blue arrows; and (3) a downstream boundary (stage hydrograph) based on water observations collected at the NOAA Lewisetta, VA gauge (LWTV).
Water level observations and predicted tides from NOAA stations were acquired from the Center for Operational Oceanographic Products and Services (NOAA Tides and Currents) to assess the performance of the model. The positions of the NOAA water level and tide stations are indicated in Fig. 1b, represented by black and light blue dots, respectively. A comprehensive account of each recording station utilized in this study can be found in Table A1 of the appendix.
2.2. Modeling Framework
In this study, we employed a model framework whose schematic structure is illustrated in Fig. 2. The components of the model framework were interconnected in an offline mode through boundary linkage. The coastal contribution, particularly astronomical tides and storm surges (referred to as storm tide), were provided by the downstream input location. The river inflow locations accounted for the contribution from rainfall-runoff processes. Additionally, the urban inflow was considered as lateral flows into the model domain. To incorporate local changes in water levels at the NOAA Washington, DC tide gage, we utilized local wind speed observations from the NOAA Washington, DC meteorological station and the Washington National Airport (WAS) recording station. These observations were uniformly interpolated over the entire modeling domain.
2.2.1. Advanced Circulation (ADCIRC)
The Advanced Circulation (ADCIRC) model, originally developed by Luettich et al. in 1992 (R. A. Luettich et al., 1992), is a finite-element hydrodynamic model that employs the generalized wave continuity equation (GWCE). By solving these equations on an unstructured computational grid in both space and time, the ADCIRC model can effectively simulate the dynamics of open water bodies such as oceans, lakes, and rivers. The model is capable of accurately capturing the effects of various driving forces, including astronomical tides, coastal storms, and river inflows. Its versatility has made it a popular choice for modeling historical storm surges and forecasting floods (Blain et al., 2010; Dresback et al., 2013; Funakoshi et al., 2012; J. L. Garzon et al., 2018; Hanson et al., 2013; Khalid & Ferreira, 2020; Shen et al., 2006).
In our study, we utilized the two-dimensional, depth-integrated version of ADCIRC, known as ADCIRC-2DDI, operating in the barotropic mode with a constant density assumption. This configuration allowed us to simulate the combined impacts of astronomical tides, river inflows, urban runoff, local winds, and storm surges on the overall water levels. ADCIRC is implemented using FORTRAN and is an open-source numerical model that benefits from extensive documentation available in both the published scientific literature (R. Luettich & Westerink, 2004a) and the official ADCIRC website (http://adcirc.org/).
2.2.2. Model Input Forcing
The study utilizes four primary input variables for model simulation, namely downstream water levels, upstream major river discharges, stream flows, and wind observations obtained from the WASD, which covers the entire modeling domain. The downstream boundary is determined based on observed water levels at the LWTV gauge. In cases where LWTV data is unavailable, data from the SWPV gauge is used instead, adjusted for a time lag of -5 hours. The upstream discharge boundary is derived from daily-observed flow data at the USGS Little Falls Pump Station gage (LFMD), which is interpolated to hourly flow using spline interpolation to ensure model stability. Additionally, urban runoff (defined as lateral flows) is considered as a boundary condition, but only when available for historical validation events. Prior to 2008, wind observations at the WASD station were not accessible. Therefore, wind data from the Washington National Airport meteorological station were utilized, with adjustments made by converting from a 5 m elevation to a 10 m elevation using a multiplier of 1.09, as described by Hsu et al. (Hsu et al., 1994). No wind observations were available for model validation before 1937. It is important to acknowledge that due to the interpolation of daily to hourly data and the use of proxy data when observed data is absent, there is a potential for uncertainty in the model results during the validation process.
2.3. Model Setup
The unstructured computational grid for the ADCIRC model was created using the automated mesh generator, OceanMesh2D (Roberts et al., 2019). To accurately represent the coastline in the numerical grid, the high-resolution coastline data from the Global Self-consistent, Hierarchical, High-resolution Geography Database (GSHHG) (Wessel & Smith, 1996) was manually updated within the study domain. The open ocean boundary, located near LWTV (Fig. 1), was defined as a non-periodic, time-varying elevation boundary condition. By placing this ocean boundary near LWTV, we were able to incorporate coastal forcing directly from NOAA's station, minimizing the error in tidal predictions resulting from low-resolution Global Tidal Models.
For validation storms, the river discharge boundary was added as a time-dependent flux boundary condition at LFMD. Since ADCIRC interprets flow as flux (Pandey et al., 2021), the discharge data was converted to m2/s using the width between the flow boundary nodes. The topography data in the localized model was extracted from the USGS 1/9 arc-second Digital Elevation Model (DEM), while bathymetry data was obtained from various sources including NOAA nautical charts (Austin, 2005), NOAA National Centers for Environmental Information (NCEI) (Caldwell et al., 2015), and the Coastal National Elevation Database (CoNED) (Thatcher et al., 2016) during the calibration phase to determine the most accurate bathymetry. All elevation and bathymetry datasets were in meters and vertically adjusted to the North American Vertical Datum (NAVD88) before being used in the ADCIRC model.
The land cover in the modeling setup was based on the National Land Cover dataset (NLCD) of 2019 (Dewitz, 2021). ADCIRC's f13 utility (https://adcirc.org/) was employed to generate nodal attributes such as mannings_n_at_sea_floor, primitive_weighting_in_continuity_equation, surface_canopy_coefficient, and surface_directional_effective_roughness_length. In our study, the simulations included various ADCIRC model features such as the wetting and drying algorithm, non-linear bottom friction, advection, finite amplitude terms, convective acceleration, and the time derivative of convective acceleration.
2.4. Model Calibration
To achieve the study objective, we performed four set of model simulations to determine the best ADCIRC modeling setup to simulate astronomical tides within the study area. These experiments focused on testing different configurations of four modeling input parameters namely, numerical mesh resolution, bathymetry data source, bottom friction, and eddy horizontal viscosity.
The testing range of values for calibration of bottom friction variable manning’s n in waterways (0.01 to 0.02), eddy horizontal viscosity (0 to 40) and grid resolution for bathymetric representation in channels (500 m to < 50m) were based on test values provided by previous literature (Bacopoulos & Hagen, 2017; Bakhtyar et al., 2020; Bastidas et al., 2016; J. Garzon & Ferreira, 2016; Kerr et al., 2013; Mied et al., 2006; Thomas et al., 2021). For bathymetry testing, we used three publicly available data sources for Potomac River (Modified NOAA nautical charts, NCEI, CONED) along with a modified mosaic of available bathymetry datasets (Merged DEM). Table 1 lists parameter-testing values used for each calibration test. Figure A1 in appendix provides the comparison of model bathymetries. The downstream ocean boundary was set as NOAA predicted astronomical tides from LWTV starting at January 1, 2020 and ending at January 30, 2020, and each calibration test simulation lasted for 30 days.
Table 1
Calibration Parameters and associated average percentage error for each test in Experiment 1
Parameters | Identifier | | Description | Nodal information |
Numerical Model Resolution | low_noCh | | Low resolution with no defined channels | NC: 6253, MiG: 450, MaG: 550 |
low_Ch | | Low resolution with defined channels | NC 8068, MiG: 350, MaG: 450 |
high_noCh | | High resolution with no defined channels | NC: 133101, MiG: 60, MaG: 225 |
| high_Ch | | high-resolution numerical mesh (high_Ch) | NC: 276647, MiG:15, MaG: 200 |
Bathymetry Datasets | B1 | | Modified NOAA nautical charts (NOAA) | Same as high_Ch |
B2 | | National Centers for Environmental Information (NCEI) | Same as high_Ch |
B3 | | Coastal National Elevation Database (CoNED) | Same as high_Ch |
B4 | | Merged DEM | Same as high_Ch |
Bottom Friction (Manning’s n) | M1 | | 0.02 | Same as high_Ch |
M2 | | 0.015 | Same as high_Ch |
M3 | | 0.013 | Same as high_Ch |
M4 | | 0.01 | Same as high_Ch |
Horizontal Eddy Viscosity (m2/s) | E1 | | 0.5 | Same as high_Ch |
E2 | | 1 | Same as high_Ch |
E3 | | 5 | Same as high_Ch |
E4 | | 10 | Same as high_Ch |
* NC: Total ADCIRC node count; * MiG: ADCIRC minimum node spacing (meters); * MaG: ADCIRC maximum node spacing (meters) |
2.5. Model Validation
For model validation, we employed the optimal mesh resolution, bathymetry sources, bottom friction, and eddy horizontal viscosity parameters established during the model calibration phase. We conducted simulations for six historical storm events using available data on river discharge, ocean water levels, and meteorological conditions as boundary inputs to assess the model's accuracy in replicating historical total water levels at the NOAA Washington, DC tide gauge, often referred to as WASD.
These case studies encompass two historical events for each flooding source category: River, Coastal, and Compound. A "River" flood event is characterized by a predominant influence of the upstream discharge from the USGS Little Falls Pump Station gauge (LFMD) on the observed water level increase at WASD. In contrast, a "Coastal" flood event results from a pronounced storm tide signal recorded at LWTV, which propagates upstream to affect WASD. Coastal floods typically involve relatively low discharge from the LFMD gauge (below 4000 m3/s).
Lastly, a "Compound" flood occurs when both high river discharge and elevated coastal water levels coincide, leading to a "Minor" flooding alert issued by the NWS (defined as 0.85 m water level above NAVD88 at WASD). These model validation events were simulated for approximately 23 days, with the simulations commencing at least 10 days prior to the observed peak water levels at WASD to allow for model initialization and enhanced stability.
Table 2 provides a detailed breakdown of event types, observed maximum water levels, wind conditions, discharge data, and event timing. Additional information regarding the selected validation events can be found below.
Table 2
Validation Events for evaluating recreation of historical water levels
Events | Type | Year | Dates | Avail. Data | Flooding Level | Observed Maximum |
| | | | | | Water level (m) | Discharge (cfs) | Wind speeds (m/s) |
Great flood of 1936 | River | 1936 | 03/14 − 03/24 | Q, H* | Major | 2.79 | 14500 | - |
Hurricane Agnes | River | 1996 | 01/14 − 01/25 | Q, H*, W+ | Major | 2.22 | 9458 | 14.9 |
April 2011 | Compound | 2011 | 04/08 − 04/27 | Q, H, W | Moderate | 1.35 | 4585 | 10.5 |
May 2014 | Compound | 2014 | 05/08 − 05/27 | Q, H, W | Moderate | 1.33 | 4073 | 10.4 |
Hurricane Isabel | Coastal | 2003 | 09/15 − 09/27 | Q, H, W+ | Major | 2.7 | 4348 | 20.1 |
TS Ernesto | Coastal | 2006 | 08/28 − 09/05 | Q, H, W+ | Major | 1.61 | 513 | 14.9 |
Q = Discharge at LFMD, H = Water levels at LWTV, H* =Water levels at SWPV, W = Wind speeds at WASD, W+ = Wind speeds at DCA |
In the early 1900s, riverine floods were common in the Washington, DC area and led to significant historical floods, including the Chesapeake-Potomac Hurricane of 1933, the Great Potomac Flood of 1936, the Flood of 1942, and Hurricane Agnes in 1972 (National Capital Planning Commission, 2008). The USGS Little Falls Pump Station gauge has continuously recorded daily river discharge since the 1930s. The Great Flood of 1936, resulting from snowmelt and heavy rainfall, saw LFMD discharge exceed 14,500 m3/s, causing severe flooding in Washington, DC, with observed water levels over 2.79 m above NAVD88 (Winters, 2018). The Blizzard of 1996, characterized by heavy snow followed by rain, led to flash floods, with LFMD recording a peak discharge of 9,260 m3/s and resulting in a maximum water level of 2.04 m above NAVD88 at the WASD gauge station.
Coastal flooding in Washington, DC is primarily caused by strong storm tide signals from the NOAA Lewisetta station in the Chesapeake Bay. Two notable coastal events, Hurricane Isabel in 2003 and Tropical Storm Ernesto in 2006, caused flooding at WASD due to recorded water levels exceeding 1 m at LWTV. Coastal floods at WASD are generally smaller than riverine floods, except for Hurricane Isabel, which led to "major" flooding at WASD with levels reaching 2.7 m above NAVD88 (Montgomery et al., 2006). Tropical Storm Ernesto generated a record storm tide of 1.5 m above NAVD88 but resulted in moderate flooding at WASD, reaching 1.61 m above NAVD88 (Knabb et al., 2006).
Compound floods in Washington, DC result from simultaneous high river discharge and coastal water levels. Major compound flood events occurred in 1937 and 1972, but they were not simulated due to data unavailability. Two distinct compound events in April 2011 and May 2014 were used for validation, combining a large storm tide from LWTV with high river discharge (> 4,500 m3/s) from LFMD, resulting in water levels at WASD rising to 1.35 and 1.3 m above NAVD88, classified as moderate flooding by the NWS. Wind speeds during both events exceeded 10 m/s, with gusts reaching up to 15 m/s at WASD.
2.6. Model Forcing Sensitivity
To assess the impact of various flood drivers on the total water level in the Potomac River, we conducted a sensitivity analysis that encompassed the following scenarios: 1) downstream boundary conditions (storm tide); 2) upstream boundary conditions (major river discharges); 3) lateral boundary conditions (urban runoff); 4) combined upstream and lateral boundary conditions; and 5) local-scale wind forcing. Each scenario involved unsteady ADCIRC simulations spanning 15 days, which included a 3-day spin-up period. These simulations were compared to a baseline scenario representing calm day conditions, driven by tidal influences alone.
For the downstream boundary conditions, we applied storm tide values based on NOAA's calculated 25-, 50-, and 100-year return periods at LWTV, resulting in surge values of 1.1 m, 1.2 m, and 1.3 m above the NAVD88 datum, respectively. To isolate the influence of downstream boundary conditions, we excluded stream flow and local wind forcing from these simulations.
To quantify the impact of upstream major river flows, we conducted simulations with discharge inputs based on stream flows corresponding to 25-, 50-, and 100-year return periods at LFMD and Anacostia river gages, sourced from USGS StreamStats analysis (Ries III et al., 2017).
For analyzing the effects of urban runoff, we introduced discharge inputs as lateral flows at various urban stream locations along the Potomac River, with return period flows of 25-, 50-, and 100-year, also based on USGS StreamStats analysis.
We further performed simulations with combined flows from major rivers and urban runoff to evaluate the cumulative impact on water levels when both factors contributed to Potomac River flow. In these combined forcing simulations, no time lag was considered to represent worst-case compound flooding conditions.
Lastly, to investigate the influence of local wind changes on water level variations at WASD, we conducted tests with eight wind directions (N, NE, E, SE, S, SW, W, and NW) and wind magnitudes ranging from 5 to 35 m/s, with a peak duration of 12 hours. These wind magnitudes were based on the wind records at WASD and the WAS meteorological station, taking into account predominant wind directions and magnitudes exceeding 35 m/s.
Table 3
Flow characteristic for major rivers and urban runoffs
| Full name | Code | Drainage Area (km2) | Max measured flow (m3/s) | 25 year return period | 50 year return period | 100 year return period |
Major Rivers | Little Falls at Potomac River | RF1 | 29940 | 13705 | 8680 | 10752 | 12908 |
Bladensburg at Anacostia River | RF2 | 239 | 349 | 233 | 317 | 580 |
Urban Runoffs | Rock Creek | LF1 | 521 | 77 | 330 | 417 | 518 |
Four Mile Run | LF2 | 122 | 36 | 308 | 390 | 483 |
Oxon Run | LF3 | 98 | - | 121 | 160 | 209 |
Cameron Run | LF4 | 228 | 116 | 349 | 426 | 512 |
Broad Creek | LF5 | 173 | 35 | 166 | 218 | 283 |
Piscataway Creek | LF6 | 420 | 127 | 231 | 300 | 384 |
Little Hunting Creek | LF7 | 65 | - | 23 | 31 | 39 |
The input conditions employed in these scenarios are theoretical and represent simplified representations of actual situations, spanning from typical daily weather to severe weather extremes. These scenarios serve the purpose of enabling a more precise assessment of the influence of each flood driver on total water levels. They underscore the significance of either incorporating or omitting specific factors to attain accuracy in total water level modeling in the Potomac River, particularly at the WASD tide gauge.
2.7. Ranking Flood Drivers
Following the aforementioned experiments, we conducted a detailed analysis to isolate the percentage change in water levels above normal daily tides at the WASD gauge station. This allowed us to quantify the individual contributions of each flood driver. In this process, flood drivers that exhibited the highest percentage change in a given return period experiment were assigned a higher rank in comparison to the other flood drivers. This ranking system provides valuable insights into the relative impact of different drivers on total water levels at WASD.
2.8. Model Evaluation
For a quantitative assessment of the modeled water levels, we employ various statistical metrics, which encompass bias, mean absolute error (MAE), and percentage change. Bias assesses the model's propensity to either overestimate or underestimate peak water levels. MAE represents the average of all absolute errors computed during the storm's duration. Percentage change is employed to gauge the variation in the modeled peak water level in comparison to the maximum observed water level at WASD. These evaluation metrics are defined as follows:
$$Bias= {x}_{mod} -{x}_{obs}$$
$$MAE= \frac{1}{n}{\sum }_{i=1}^{n}| {x}_{mod} -{x}_{obs}|$$
$$Percent Change= \frac{{x}_{mod} -{x}_{obs}}{{x}_{obs}}* 100$$
Where \({x}_{mod}\)and \({x}_{obs}\) are the modeled and observed water level at a given recording station at a given time step. The \(n\) is the number of time steps in hours for error calculation.