A coupled hydrodynamic (HEC-RAS 2D) and water quality model (WASP) for simulating flood-induced soil, sediment, and contaminant transport.

Increased intensity and frequency of floods raise concerns about the release and transport of contaminated soil and sediment to and from rivers and streams. To model these processes during flooding events, we developed an External Coupler in Python to link the Hydrologic Engineering Center-River Analysis System (HEC-RAS) 2D hydrodynamic model to the Water Quality Analysis Simulation Program (WASP). Accurate data transfer from a hydrodynamic model to a water quality model is critical. Our test results showed the External Coupler successfully linked HEC-RAS and WASP and addressed technical challenges in aggregating flow data and conserving mass during the flood event. We ran the coupled models for a 100-year flood event to calculate flood-induced transport of sediment-associated arsenic in Woodbridge Creek, NJ. Change in surface sediment and arsenic at the end of 48-h flood simulation ranged from a net loss of 13.5 cm to a net gain of 11.6 cm, and 16.2 to 2.9 mg/kg, respectively, per model segment, which demonstrates the capability of the coupled model for simulating sediment and contaminant transport in flood.

INTRODUCTION

Floods are one of the most common natural disasters, causing significant loss of life and property in the United States. During the 2018 water year (October 1st, 2017 to September 30th, 2018), floods caused 84 deaths and one billion USD in damages to buildings, infrastructure, and croplands in the United States (National Weather Service, 2018). In addition to direct damages, floods have caused unintended release and dispersal of contaminated soils and sediments from contaminated sites (Horney et al., 2018; Su et al., 2008). Flooding from Hurricane Katrina (2005) damaged 21 oil refinery facilities in southeastern Louisiana and caused widespread sediment deposition in flood-inundated areas (Rotkin-Ellman et al., 2010). A post-Hurricane Katrina toxic trace assessment in surface sediment samples found arsenic and vanadium contamination in flood-affected areas of New Orleans, Louisiana (Su et al., 2008). Hurricane Harvey (2017) impacted Superfund sites in Texas and Louisiana. Among these sites, flooding damaged the San Jacinto River Waste Pits Superfund site cap and released underlying dioxins stored in the soil layer. As a result, cap repairs were made to reduce the dioxin concentration to the EPA recommended clean-up level for the site (U.S. Environmental Protection Agency Region VI, 2017). The role of fluvial flooding in contaminated sediment transport is also a concern, such as repeated flooding of the Chattanooga Creek, TN, causing mobilization of polycyclic aromatic hydrocarbons in the floodplain (Vulava et al., 2017). As the trend of increasing intensity and frequency of floods continues (Di Baldassarre et al., 2010; Vogel et al., 2011), there are increasing concerns of release and transport of contaminated soil and sediment (United States Government Accountability Office, 2019).

Soils and sediments act as sinks for contaminants in the environment because their organic and inorganic components provide attractive binding sites for low solubility toxic chemicals (Wölz et al., 2010). Binding to soil and sediment particles often immobilizes the contaminants in the environment by protecting them against external forces like transport of surface water (Brinkmann et al., 2010). In flood conditions, the opportunity for detachment and transport of contaminated sediment and soil particles increases due to the significant velocity and shear stress of floodwaters. The importance of sediment particles to contaminant transport under flood conditions has led to multiple studies on the topic (Baborowski et al., 2004; Foulds et al., 2014; Macklin et al., 2003; Miller et al., 1999; Moody et al., 2000; Navrátil et al., 2008; Vulava et al., 2017). A comprehensive list can be found in Ciszewski and Grygar (2016). However, a limited number of studies have used multidimensional hydrodynamic and water quality models to physically simulate and predict sediment and contaminant transport in flood events (Ciffroy et al., 2000; Hardy et al., 2000; Hostache et al., 2014; Schulz et al., 2009).

Typical two-dimensional (lateral and longitudinal) unsteady flood models use a fine-scale computational grid size on the order of meters with time steps on the order of seconds to accurately capture the complexity of river hydraulics (Bomers et al., 2019). Water quality models typically do not require such fine temporal and spatial resolution (Ganju et al., 2016). The resolution disparity causes flood models to become computationally expensive when solving water quality equations. By linking flood and water quality models, hydrodynamic data can be aggregated and used for water quality simulation while reducing computational cost (Ganju et al., 2016).

We coupled the Hydrologic Engineering Center-River Analysis System (HEC-RAS) 2D and the Water Quality Analysis Simulation Program (WASP) to simulate flood-induced sediment and contaminant transport at contaminated sites. HEC-RAS is a commonly used hydrodynamic model that simulates flood hydrodynamic routing in one- and two-dimensions (Brunner, 2016). The model is an approved tool for studying flood hazards designated by the U.S. Federal Emergency Management Agency (FEMA) (U.S. Federal Emergency Management Agency, 2020) and the U.S. Army Corps of Engineers (Brunner, 2016). HEC-RAS has been coupled with Qual2k for water quality studies (Fan et al., 2009) and Modular Three-Dimensional Finite-Difference Ground Water Flow Model (MODFLOW) to simulate surface water and ground water interaction (Rodriguez et al., 2008). Several programs such as HEC-RAS API controller, PyRAS (Peña-Castellanos, 2015), and Pyflood (Vimal, 2015) also have been developed to externally run and manipulate the model’s input and output, details of these programs can be found in Goodell (2014a), Goodell (2014b), and Dysarz (2018).

WASP is a fate and transport water quality framework developed by the U.S. Environmental Protection Agency (Ambrose & Wool, 2017) and is used worldwide (Knightes et al., 2019). WASP has been used in several sediment and contaminant studies (Carroll et al., 2000; Caruso, 2004; Franceschini & Tsai, 2010; Schultz, 2001; Suk & Fikslin, 2011). WASP does not simulate multidimensional hydrodynamics, but its flexible structure allows for importing hydrodynamic data from an external model. Although WASP accepts hydrodynamic data from several models (EFDC, DYNHYD, RIVMOD, CE-QUAL-RIV1, and SWMM) (Ambrose & Wool, 2017), and other linkages have been developed (Defne et al., 2017; Rifai et al., 2013; Wool et al., 2003), a coupler between HEC-RAS 2D and WASP framework has not been previously available. A coupled HEC-RAS 2D and WASP provides an opportunity to study fate and transport of sediment and contaminant in rivers and floodplains, which was not possible via the previously developed WASP linkages.

The objective of this study was to couple HEC-RAS 2D and WASP to simulate flood-induced fate and transport of sediment and contaminants and demonstrate its application at a contaminated site. We hypothesize that coupling these two models will successfully simulate sediment and contaminant transport in dynamic flood conditions. The modeling results can provide information for contaminated site managers and associated communities for building resiliency. This research presents the HEC-RAS 2D to WASP coupling process and an example application for Woodbridge Creek in New Jersey, US. The linked HEC-RAS 2D and WASP were used to simulate fate and transport of sediment and contaminant (arsenic) during the 100-year flood event.

METHODS

HEC-RAS 2D description

HEC-RAS (v5.0.7) simulates steady and unsteady flow 1D and 2D. The 2D module simulates flow hydraulics in the river channel and over the floodplain where the flow is mostly two-dimensional (Blanch, 2017). The model uses computational cells to represent terrain for calculating flow and velocity along a horizontal plane, from a computational cell into adjacent cell. HEC-RAS 2D calculates flow rate for a cell boundary using hydraulic properties of the cell and water surface elevation of adjacent cells (Brunner, 2016). The model solves the 2D shallow water equations (SWE) or diffusion wave equations (DWE) with an implicit finite volume solution algorithm. SWE is derived using continuity and momentum equations, and DWE is an approximation of SWE obtained by neglecting the inertial terms of the momentum equations (Moya Quiroga et al., 2016). Compared to DWE, SWE is more momentum conservative and it can model turbulence and Coriolis effect. Both methods, however, are expected to result in accurate simulations, except SWE should be used in situations when there are extreme changes in flood wave and velocity (Brunner, 2016). Within the HEC-RAS, DWE is the default method since it improves model stability by simplifying calculations (Blanch, 2017). Since there are no drastic changes in flood wave and velocity during our flood simulation, we used DWE for HEC-RAS calculations.

HEC-RAS 2D uses spatially varying Manning’s roughness coefficients to calculate the flow velocity. HEC-RAS v5.0.7, used here, does not support sediment scour and deposition in two dimensions. However, the recently released HEC-RAS v6.0 (Beta) supports 2D sediment scour and deposition simulation and simulates 1D transport of some water quality parameters (e.g. dissolved oxygen and phosphorus). For this application, we linked HEC-RAS v5.0.7 to WASP and simulated fate and transport of sediment and associated toxicants in WASP.

WASP description

WASP is a widely used mass balance fate and transport modeling framework for simulating concentrations of contaminants in sediment and surface water (Ambrose & Wool, 2017; Arlos et al., 2014; Bouchard et al., 2017; Carroll et al., 2000; Wool et al., 2003; Yang et al., 2012). The framework applies a differential mass balance equation for simulating advective and dispersive transport of a constituent in 1D, 2D, or 3D. WASP 8.4 includes the Advanced Toxicant module, which uses hydrodynamic data to calculate fate and transport of sediment and chemicals in surface waters and sediment layers. To simulate dissolved and suspended materials in 2D the model uses hydrodynamic data to solve the 2D advection-dispersion equation (Wool et al., 2020): where C is the concentration (M/L3); t is time; U and U are velocities (L/T) in the x and y direction; E and E are longitudinal and lateral diffusion coefficients (L2/T); S is direct and diffuse loading rate (M/L3T); S is boundary loading rate (M/L3T); and S is total kinetic transformation rate (M/L3T). Positive terms indicate sources and negative terms indicate sinks.

WASP sediment transport module allows simulation of up to ten solids classes with discrete particle size and density in descriptive and mechanistic approaches (Knightes et al., 2019). In the mechanistic approach used in this study, the water shear stress and particle size mainly control erosion, deposition, and resuspension of the sediment, while settling is a function of particle size and density.

Depending on the magnitude of shear stress generated by the flood and the soil and sediment particle size, erosion can occur in cohesive and non-cohesive forms (Mitchener & Torfs, 1996). WASP calculates erosion and resuspension velocities and fluxes for cohesive and non-cohesive erosion and for mixtures of cohesive and non-cohesive sediments. For cohesive sediment, the model calculates erosion flux from the excess shear stress power law equation (Lick et al., 1994); where Ecoh is the cohesive erosion flux (g/m2 s); Fcoh is the fraction of the surface bed that is cohesive; M is the shear stress multiplier (g/m2 s); n is the shear stress exponent; and τ* is the excess shear stress (N/s): where τ and τce are water shear stress and critical shear stress (N/m2) for cohesive sediment erosion. For non-cohesive sediment, WASP calculates erosion velocity and flux for each particle class separately. For details of WASP sediment flux and sediment transport calculations see Ambrose and Wool (2017).

WASP simulates sorption using equilibrium or kinetic models. In the equilibrium model applied in this study, the interaction between chemicals and solids is assumed to be fast relative to other environmental processes such as oxidation and photo-transformation (Ambrose et al., 2017). At equilibrium, concentrations of a chemical in dissolved and solid phases are determined by the partition coefficient (Equation (5)) to solve transport equation for chemicals (Equation (1)): where C and C are the concentrations of the chemical in dissolved and solid phases (mg/L); and K is the partition coefficient (L/kg).

HEC-RAS 2D to WASP external coupler

WASP requires water volume, depth, velocity, and average interfacial flow at each time step to calculate sediment and contaminant transport in 2D. The model can obtain the 2D hydrodynamic data through a linkage file generated by the Hydrolink Application Program Interface (API). The Hydrolink API is included in the WASP package; it allows users to convert arrays of hydrodynamic information into binary input for the model (Ambrose & Wool, 2017). To accurately extract the hydrodynamic output from the HEC-RAS 2D flood simulation and use it as an input for WASP, we developed an External Coupler using Python (v.3.7) to link the two models via the Hydrolink API. The External Coupler connects HEC-RAS 2D to WASP in an offline and one-way direction by: (a) reading hydrodynamic outputs from HEC-RAS 2D, (b) aggregating those outputs into larger cell sizes suitable for WASP, and (c) converting the aggregated data into a hydrodynamic input file for WASP using the Hydrolink API.

Differences in WASP and HEC-RAS 2D spatial resolution

WASP uses explicit finite difference methods to solve mass balance and fate and transport equations (Ambrose & Wool, 2017). Providing HEC-RAS 2D hydrodynamic flow data at fine spatial and temporal time scales may cause numerical errors for WASP. In addition, this setup would result in excessively long runtimes because WASP uses a single processor to solve fate and transport equations (Defne et al., 2017). To create a computationally robust model, the HEC-RAS 2D computational grid cells should be aggregated for use in WASP.

HEC-RAS 2D reports velocity and flow components for cell faces. To preserve the fidelity of these data for WASP the aggregated cells (segments) must align with the underlying HEC-RAS 2D cells at their boundaries. We developed an automated aggregation routine (included in the External Coupler) to read HEC-RAS 2D cells and WASP segments in ESRI shapefile (ESRI, 1998) format and align the models’ grids by spatially intersecting HEC-RAS 2D cells’ centroid and segment polygons. That is, WASP segment boundaries are modified based on the HEC-RAS 2D cells within the segment. The program creates a connectivity table to catalogue mapping between WASP segments and HEC-RAS 2D cells to transfer data between the two models. The External Coupler also generates a flow direction map illustrating the neighboring information for WASP segments to simulate the 2D flow transport within the grid.

Spatially aggregating HEC-RAS 2D flood hydrodynamic data for WASP

We calculated WASP segment volume in two steps. First, segment volume at time = 0 was calculated from the sum of HEC-RAS 2D cells’ volumes estimated from the water surface elevation output. For this purpose, a polynomial interpolation method was implemented within the coupler using the SciPy Python library (v.1.2.1). This method uses water surface elevation and estimates the volume for the HEC-RAS 2D cells from their elevation-volume table. The time series of volume for a WASP segment was calculated from the net inflow to the segment using Equation (6): where V (m3) is the volume for the WASP segment; t is time; Qin and Qout are inflow and outflow (m3/s) to the segment.

WASP is not designed to simulate dry segments. Since WASP is designed to calculate concentration by dividing mass by volume, if the segment volume equals zero a fatal error occurs. This is a significant challenge, as rapid wetting/drying processes are possible throughout the duration of a flood. To address this issue, the coupler assigns a minimum water volume for dry segments to allow WASP to operate. This volume should be selected with the objective of maintaining the smallest volume possible (representative of a dry segment) while conserving mass for WASP segments.

WASP uses water velocity to calculate shear stress, which determines the erosion and deposition velocities and fluxes for the segment. To calculate the velocity for WASP input, the vector of HEC-RAS 2D face velocities are weighted by their hydraulic areas and averaged to compute the cell velocity (Equation (7)). The calculated HEC-RAS 2D cell velocities (derived from cell-face velocities) are then weighted by the cell area and averaged to calculate an aggregate velocity for each WASP segment (Equation (8)); where u and u are cell and segment velocities (m/s); A is hydraulic area of the cell face; A is area of cell (m2). A similar approach, weighted average by the cell area, is used to calculate water depth for WASP segments derived from the HEC-RAS 2D water depth of each cell.

WASP requires interfacial flow and its direction in each time step to calculate advective mass transport between segments. The velocity and hydraulic area at each cell face are used to calculate the interfacial flow, Q (m3/s), at time t for HEC-RAS 2D cells: The interfacial flow for the WASP segments was computed by summing the HEC-RAS 2D interfacial flows at the segment boundaries. Since WASP requires integrated or average interfacial flow across the time step, a linear interpolation (averaging the flows at output time t and t − 1) method was used to integrate the calculated flow.

Generating the hydrodynamic linkage

To generate an input file for WASP, the coupler organizes the aggregated flood hydrodynamic data into an input file for the Hydrolink API. The following data are created within the input file: (a) model set up information, such as number of segments, number of flow directions, time step, and time stamps of the flood simulation; (b) the flow direction, including the location of upstream and downstream boundaries for the model; (c) HEC-RAS 2D aggregated hydrodynamic outputs; (d) initial value volume for the surface sediment layer. By calling the Hydrolink API, the coupler-generated input file will be converted to a binary input file for WASP.

Case study

We applied the coupled HEC-RAS 2D and WASP complete modeling approach in the Woodbridge watershed, New Jersey, USA. The approach for simulating flood-induced sediment and contaminant transport consists of five steps: (a) data collection, (b) flood discharge simulation using the Hydrologic Engineering Center—Hydrologic Modeling System (HEC-HMS) watershed model, (c) flood hydrodynamic simulation using HEC-RAS 2D, (d) linking HEC-RAS to WASP using the External Coupler, and (e) simulating sediment and contaminant transport using WASP.

Study area

The Woodbridge watershed is in Perth Amboy, NJ and has a drainage area of 44.5 km2 (Figure 1(a)). Part of the watershed (5.3 km2) is a designated FEMA Special Flood Hazard Area (SFHA) because of its vulnerability to tidal and fluvial flooding (U.S. Federal Emergency Management Agency, 2014). Woodbridge Creek and its tributary, Spa Spring Creek, empty into the Arthur Kill tidal strait (Figure 1(b)).

FIGURE 1

Map of the study area. (a) Location of the watershed in New Jersey, US; (b) elevation of the Woodbridge Creek watershed based on the Digital Elevation Model (Office for Coastal Management, 2020), model domain, and FEMA SFHA (U.S. Federal Emergency Management Agency, 2019); (c) land use/land cover (New Jersey Department of Environmental Protection, 2015) and model domain. Background image source—ESRI, GEBCO, NOAA, National Geographic, DeLorme, HERE, Geonames.org, and other contributors

Soils and sediments near Woodbridge Creek’s mouth are contaminated with toxic chemicals (TRC Environmental Corporation, 2016) from past industrial activities, primarily oil and gas production and storage. The location of soil and sediment contamination in the SFHA suggests potential for mobilization, transport, and dispersal of the contaminants during flood conditions.

We confined the modeling domain to 4.0 km upstream of the Arthur Kill along Woodbridge Creek. This corresponds to locations where the contaminated soils and sediments are found in the floodplains and river channel (Figure 1(c). The modeling domain is 5.8 km2 and characterized by relatively flat terrain (Figure 1(b), with a mean slope of 4.0%. Industrial and residential land uses are dominant, covering 70% of the model domain (Figure 1(c)). Average annual precipitation is 1025 mm with greatest rainfall in spring and summer months (Middlesex County Hazard Mitigation Plan, 2015). Tropical storms also affect the watershed. The 100-year, 24-h storm corresponds to a precipitation depth of 255 mm (National Oceanic and Atmospheric Administration, 2017). Flashy stormwater flows in the watershed result in high flood elevations, increased velocities, and wider floodplains in Woodbridge Creek and Spa Spring Creek (Township of Woodbridge Floodplain Management Plan, 2015).

HEC-HMS flood discharge simulation

Time-varying discharge at the upstream boundaries was simulated using HEC-HMS. Using the SSURGO soil dataset (USDA NRCS, 2019) and a 2015 land cover map (New Jersey Department of Environmental Protection, 2015) as inputs, we simulated the 100-year flood discharge for Woodbridge Creek and Spa Spring Creek. Flood discharge was simulated using the Soil Conservation Service Runoff Curve Number Method (USDA NRCS, 1986) for a 24-h storm with a 100-year return period (National Oceanic and Atmospheric Administration, 2017). HEC-HMS simulated flood discharge was calibrated against the FEMA’s steady 100-year flood peak (U.S. Federal Emergency Management Agency, 2014) by adjusting runoff curve numbers within 10% range of default value.

HEC-RAS 2D set up for flood simulation

HEC-RAS 2D was set up using a LiDAR Digital Elevation Model (DEM) with 1 m spatial resolution (Office for Coastal Management, 2020). The DEM was modified by superimposing river channels and inserting aboveground structures (e.g., buildings and storage tanks) that were manually digitized using 60 cm resolution ESRI World Imagery in ArcGIS 10.7 (ESRI, 2019) to accurately capture surface topography and flow paths. A computational grid composed of 34,806 rectangular and irregular (greater than four faces/sides) cells of a nominal size between 5 and 30 m was generated to capture variability of the underlying terrain for the 2D flood simulation.

Manning’s roughness coefficients for the study area were defined based on land cover (Manandhar, 2010) derived from the 2015 New Jersey Department of Environmental Protection land cover map (New Jersey Department of Environmental Protection, 2015). Woodbridge and Spa Spring Creeks were established as two upstream boundaries for our modeling domain and their associated 100-year flood was input from HEC-HMS. The stage for the model’s downstream boundary was calculated from the Manning equation based on the simulated flows in the outlet.

We ran HEC-RAS 2D using Diffusion Wave equations with a 1-s computational interval to simulate the 100-year flood for a 48-h duration. The HEC-RAS 2D hydrodynamic simulation output was exported in 1-min output intervals, as an input to WASP. Due to the lack of 100-year floodplain inundation observations and water surface elevation measurements for Woodbridge Creek, the HEC-RAS 2D was calibrated against the SFHA 100-year floodplain extent and water surface elevation in 102 locations ((Figure 1(c)) by locally adjusting value of Manning’s roughness coefficients for river channel and floodplain.

WASP segments

According to the SFHA, 61% of the HEC-RAS 2D cells would be inundated in a 500-year flood. The WASP simulation domain was therefore confined to the 21,098 HEC-RAS cells that were inside the 500-year SFHA. This provided an upper bound on the model domain, with the potential to expand simulations to the 500-year flood. Using the External Coupler program 21,098 HEC-RAS 2D cells data were aggregated into 812 WASP computational segments ranging from 900 to 21,600 m2 (mean 2500 m2). The WASP segments were aligned to HEC-RAS 2D cell boundaries to accurately transfer data between the two model grids (Figure 2). WASP segments were set up using surface water and surface benthic layer. A surface benthic layer with 15 cm depth was generated for WASP segments to simulate transport of sediment and contaminant between water column and surface sediment.

FIGURE 2

(a and b) HEC-RAS 2D computational grid cells (21,098) converted to WASP computational segments (812). Model domain is the FEMA 500-year floodplain SFHA

WASP setup for sediment and contaminant simulation

WASP must be parameterized with initial and boundary conditions after completing the hydrodynamic linkage. Clay, silt, and organic matter that can potentially adsorb contaminants (Ji, 2017) and fine sand served as solids classes for the WASP simulation of flood-induced sediment and arsenic transport. The initial volumetric concentrations of these solids in the bed (mg/L), within a depth of 15 cm, were defined for the model based on their percentage, soil bulk density in the SSURGO dataset (USDA NRCS, 2019), and 14 in situ surface sediment samples that were collected in Woodbridge Creek and Spa Spring Creek in 2013 and 2014 (Figure 1(c)).

Due to a lack of suspended sediment concentration measurements, the initial concentration of solids in the surface water and their loads at the upstream model boundaries were estimated from Landsat 8 imagery. For this purpose, 14 Landsat 8 images for the time period from 2013 to 2014 were processed using Acolite software (Vanhellemont, 2019) to calculate the total particulate matter (TPM) in the surface water. The maximum calculated TPM was selected as a representative for flood condition and it was partitioned into individual solids classes using the Stokes’ settling velocity for the solids (Lerman et al., 1974). This assumes the fraction of a solid in TPM is proportional to the solid’s settling velocity, calculated by dividing the solid’s settling velocity by settling velocity of the smallest particle (clay). For example, silt concentration was 25 times smaller than clay due to its larger settling velocity, 7.1 m/day versus 0.28 m/day.

We performed a one-at-a-time sensitivity analysis for WASP sediment transport by individually adjusting values of eight model parameters: critical shear stress (τce) for cohesive and non-cohesive erosion, shear stress multiplier (M) and exponent (n) for cohesive and non-cohesive resuspension, and lower and upper critical shear stress for solids deposition. The sensitivity analysis showed that the critical shear stress for cohesive erosion is the most sensitive parameter. The critical shear stress value for cohesive erosion was set to 0.03 N/m2 following Berenbrock and Tranmer (2008).

Arsenic measurements (555 samples) in soil and sediment were used to initiate initial concentrations in WASP segments (Figure 1(c)). Inverse Distance Weighting (IDW) in ArcGIS 10.7 was used to spatially interpolate arsenic samples measured within a soil and sediment depth of 15 cm. For segments beyond the IDW boundary, arsenic initial concentration was set to zero. Initial and boundary concentrations of arsenic in surface water were set to zero due to lack of water quality measurements. The partition coefficient of arsenic to soil and sediment particles was set to 10,000 L/kg (Allison & Allison, 2005). WASP was run with a 1-min computational interval to calculate the sediment and contaminant transport in the 100-year flood for a 48-h duration.

We performed an uncertainty analysis for sediment and contaminant simulations by testing a combination of high and low bounds of estimated suspended sediment (28 and 78 mg/L). Arsenic concentration at the upstream boundary was set to 0.00135 mg/L in the water column, a median value for 36 water samples measured in neighboring watershed (USGS-01395000). The results were compared with the sediment and contaminant simulation to address the sensitivity of suspended sediment and arsenic concentration in WASP.

RESULTS

Flood simulation

The calibrated HEC-HMS 100-year, 24-h flood discharge of 64.8 m3/s and 23.5 m3/s reproduced FEMA’s 100-year steady flood peak discharge (U.S. Federal Emergency Management Agency, 2014) for Woodbridge Creek and Spa Spring Creek respectively (Figure 3), reported at their junction (Figure 1(c)).

FIGURE 3

HEC-HMS calibrated 100-year flood discharge for Woodbridge (blue solid) and Spa Spring (orange solid) compared to their FEMA’s corresponding steady flood peak discharges, dashed blue and orange lines

The extent of the 100-year flood simulated in HEC-RAS 2D matched the 100-year floodplain from FEMA SFHA (Figure 4(a)) with a mean absolute error of 0.07 m in water surface elevation, minimum and maximum of 0.007 and 0.17 m, in 102 locations (Figure 4(b)). Deviation in water surface elevation and area was minimal. The differences are likely due to the HEC-RAS 2D unsteady flood simulation compared with the FEMA 1D steady flood simulation. In 2D, the flood wave is represented as a function of time and influenced by topography and structures such as buildings, better capturing the flood wave and pattern horizontally (Cook & Merwade, 2009). The simulated 100-year flood raised the Woodbridge Creek stage from 0.5 to 4.5 m above mean sea-level, which caused partial inundation of industrial and residential areas with water depths less than 2.6 m and 0.5 m, respectively (Figure 4(b)).

FIGURE 4

(a) HEC-RAS 2D simulated 100-year floodplain extent versus the FEMA floodplain extent; (b) HEC-RAS 2D simulated maximum 100-year water depth. Basemap image source—ESRI, DigitalGlobe, GeoEye, i-cubed, USDA FSA, USGS, AEX, Getmapping, Aerogrid, IGN, IGP, swisstopo, and the GIS User Community

Coupled HEC-RAS and WASP 2D sediment transport simulation in 100-year flood

Using the HEC-RAS 2D 100-year flood hydrodynamic simulation data via the coupler, WASP calculated flood-induced sediment erosion and transport for 812 model segments. Time series of flow velocity, shear stress, surface sediment volume (in the top 15 cm of the surface sediment layer), and suspended sediment concentration for two-channel segments are shown in Figure 5, in which one segment experiences erosion and the other deposition. The WASP-calculated shear stress exceeded the minimum critical shear stress for cohesive erosion (Figure 5(c)–(d)) when the floodwater velocity was ≥0.06 m/s (Figure 5(a)–(b)). By exceeding the critical shear stress for cohesive erosion, the flood scoured the sediment layer and resuspended the eroded solids for both segments, thereby decreasing the volume of the surface sediment layer when erosion was greater than deposition (Figure 5(e)). The concentration of suspended sediment reached its maximum at the flood peak, coinciding with the greatest volume of sediment eroding at peak velocity and shear stress and load from upstream (Figure 5(g)–(h)). Note, the WASP sediment transport model is parameterized based on available data, but further sensitivity analysis and calibration are required to address high concentrations of total suspended sediment in surface water.

FIGURE 5

(a, b) floodwater velocity; (c, d) flood shear stress, (e, f) surface sediment depth; and (g, h) suspended sediment concentration simulated for two WASP segments within the Woodbridge Creek channel. Left-panel figures (a, c, e, g) correspond to the erosional (upstream) segment, and right-panel figures (b, d, f, h) correspond to the depositional (downstream) segment

The concentration of suspended sediment in the surface water at the end of the 48-h flood simulation ranged from 0–2000 mg/L, of which about 80% were clay particles (Figure 6(a)). Concentration of suspended sediment increases from the river channel toward the floodplain and is greatest in model boundary segments. High concentrations of suspended sediment (>200 mg/L) are caused by clay particles trapped in dry segments at the end of the 48-h flood in HEC-RAS 2D output; the coupler forces dry segments to maintain a minimum water volume (20 m3, a value that we identified as optimal volume for dry segments based on the Woodbridge Creek model grid) in WASP. Surface sediment at the end of 48-h flood simulation ranged from a net loss of 13.5 cm to a net gain of 11.6 cm per model segment (Figure 6(b)), with an average of 3.7 cm and 0.14 cm for net loss and net gain, respectively. Changes to the surface sediment depth occurred in channels and on the floodplain. Although most segments experienced net deposition at the conclusion of the 48-h simulation, erosion was the dominant process near the upstream model boundary, where velocity and shear stress were greatest. A few segments at the model’s downstream boundary edge exhibit high erosion and deposition because these segments overlap a small strip of the Arthur Kill strait, which was included in the HEC-RAS 2D model domain.

FIGURE 6

(a) WASP simulated suspended sediment concentration in surface water at time = 48 h for Woodbridge Creek and Spa Spring Creek. Dark and light brown colors present high and low concentrations of suspended sediment, respectively; (b) change in WASP simulated sediment layer depth in 100-year flood, negative values indicate loss of sediment (erosion) and positive values indicate gain (deposition). Blue and red colors present erosion and deposition, respectively. Darker colors correspond to greater changes in sediment layer depth. Background image source—ESRI, GEBCO, NOAA, National Geographic, DeLorme, HERE, Geonames.org, and other contributors

Arsenic concentration in surface sediment at the end of the 48-h simulation is presented and compared to the contaminant initial concentration in soil and sediment in Figure 7(a)–(b). Predicted changes in arsenic concentration in surface sediment after the 100-year flood mostly ranged from −1.0 to 1.0 mg/kg (light blue and red colors in Figure 7(c)), except for the few segments at the model’s downstream boundary that experienced larger changes due to greater erosion and deposition (Figure 7(c)). Deposition of eroded sediment from the Spa Spring Creek and Woodbridge Creek channels increased arsenic concentration in middle and north parts of the study area. Expansion of the extent of arsenic contamination in the north section of the model domain occurred after the 48-h simulation where deposited sediments increased arsenic concentration above zero initial concentration (Figure 7(c)). Most of the model domain showed decreases in arsenic concentration due to deposition of previously suspended sediments with lower concentration.

FIGURE 7

(a) Initial concentration of arsenic in soil and sediment interpolated by Inverse Distance Weighting method; (b) arsenic concentration in surface sediment at the end of 48-h flood simulation; (c) difference in arsenic concentration between (a) and (b). Background image source—ESRI, GEBCO, NOAA, National Geographic, DeLorme, HERE, Geonames.org, and other contributors

There was minimal change in surface sediment (Figures 6(b) and 8(a)) and contaminant concentration (Figures 7(c) and 8(b)) with boundary conditions of 0.00135 mg/L arsenic concentration in combination with 28 mg/L and 78 mg/L of suspended sediment. The exception was a small area in the north part of the modeling domain where the concentration of arsenic in soil and sediment increased from zero initial concentration to 0.1 mg/kg due to deposition of contaminated sediments (Figures 7(a) and 8(b)).

FIGURE 8

(a and b) change in surface sediment depth and arsenic concentration at the end of 48-h flood simulated using lower bound of suspended sediment (28 mg/L) and 0.0013 mg/L arsenic for the upstream boundaries. The change in arsenic using higher bound of suspended sediment (78 mg/L) looks the same as (b) when divided into 7 bins shown in the figure. Background image source—ESRI, EBCO, NOAA, National Geographic, DeLorme, HERE, Geonames.org, and other contributors

DISCUSSION

Accurate data transfer from a hydrodynamic model to a water quality model is critical. Temporal and spatial upscaling when moving from the hydrodynamic model to the water quality model can result in a loss of hydrodynamic information (Defne et al., 2017). Temporal aggregation of the HEC-RAS 2D output in 1-min intervals had minimal impact on WASP results. This is because the calculated volumes for the WASP segments were equal when the HEC-RAS 2D instantaneous interfacial flow and integrated interfacial flow were used. However, the WASP mass balance and sediment transport results appear to be sensitive to the size of segments. The WASP segments typically had smaller velocities compared to the underlying HEC-RAS 2D cells because they contain dry cells with zero velocity. Since WASP uses velocity for calculating shear stress, the underestimation can result in lower than expected erosion and resuspension. This issue was addressed by adjusting the critical shear stress for erosion and resuspension (see Berenbrock and Tranmer (2008)) within WASP.

Results of sediment transport simulation in the 100-year flood for Woodbridge Creek showed widespread deposition in river channels and the floodplain. Sediment deposition was dominant because of the flat terrain and prominent floodplains in the study area, which reduces flood velocity and shear stress (Figure 5(a)–(d)). The maximum depth of deposited sediment in 100-year flood was 11.6 cm. This is within the range of 0 to 30 cm reported by Hostache et al. (2014) for simulating transport of fine sediments in floods with similar discharge. The depth of deposited sediment in river channels, however, appears to be small (<1 cm, Figure 6(b)) in the context of the 100-year flood magnitude (Figure 3). This is because suspended particles in the Woodbridge watershed were mainly composed of clay and silt, which were transported out of the watershed by the flood due to their low settling velocity.

Soil and sediment erosion and deposition are a primary mechanism for the transport and distribution of contaminants in floods (Hostache et al., 2014; Schulz et al., 2009). By analyzing changes in arsenic concentrations under flood conditions for Woodbridge Creek (Figure 7(c)), the results showed that flood deposited sediment causes redistribution of arsenic in the river channel and floodplain. Arsenic concentration was predominantly reduced (Figure 7(c)) because the contaminant concentration sorbed to suspended sediments was lower than surface sediment. Since most of the suspended sediment was sourced from eroding segments near the Woodbridge Creek and Spa Spring Creek boundaries and upstream boundary loads with low or zero arsenic concentration, the deposited sediment reduced arsenic concentration in the river channel and floodplain due to the mixing process. Similar results have been reported in studying distribution of heavy metals in floods (Bradley & Lewin, 1982; Ciszewski, 2001; Schulz et al., 2009).

In 2D simulations of flooding, some computational cells on the floodplain progressively wet and dry over the course of the flood. WASP is not designed to simulate wetting and drying segments and will not operate with segment volumes equal to zero. In order to address this limitation, we included a minimum volume parameter, set to 20 m3, for dry segments in this study. This minimum volume ensured WASP mass balance continuity, where the internal mass check variable remains close to one. The minimum volume should be defined by the user when the External Coupler developed here is used to link HEC-RAS 2D and WASP. This gives the user flexibility in determining the best minimum volume for their specific study. In the future, this limitation can be resolved by encoding the capability of simulating dry segments in WASP.

Calibration and validation of water quality simulations under flood conditions is a long-standing difficulty due to insufficient samples, particularly at peak discharge. With limited or no calibration and validation, a robust parameter estimation and uncertainty analysis to develop uncertainty bounds around simulation results are necessary to communicate confidence in the model output.

CONCLUSIONS

Floods play an important role in mobilization of soil and sediment and transport of contaminants. Quantifying fate and transport of sediment and contaminants under flood conditions is a challenging task because it typically requires a composite hydrodynamic and water quality simulation. The development of a HEC-RAS 2D and WASP External Coupler provides the capability to adequately model these situations. The External Coupler operates in an offline, one-way manner to accurately extract hydrodynamic data (interfacial flow, water velocity, and water volume) from HEC-RAS 2D as an input to initiate sediment and contaminant transport simulations via WASP v8.4. The research in Woodbridge Creek demonstrated the success of coupling the models to simulate flood-induced fate and transport of sediment and contaminants, with sediment-arsenic redistribution as an example.

Results of the case study of the external coupler demonstrated that shear stress, initial contaminant concentrations in surface sediment (and their spatial distribution), and suspended sediment impact distribution of sediment-associated contaminants during and after a flood event. The coupled HEC-RAS 2D and WASP models developed in this study can be used to simulate sediment, toxicants, and water quality parameters under flood conditions. In addition to fluvial flooding, we could assess climate future impacts and examine the coupled model framework under tidal, storm surge, and sea-level rise scenarios in coastal watersheds. This scenario research would add another level of complexity to the coupled models by combining fluvial flooding with dynamic downstream boundary conditions. The coupled HEC-RAS 2D and WASP model is a reliable tool to predict fate and transport of sediment and contaminants in floods, with applications for contaminated site managers and decision-makers to address flood resilience. In addition, it creates a link between two widely used models and their communities of practice.