The HEC-RAS software (Hydraulic Engineering Center, developed by the United States Army Corps of Engineers), which yields reliable results in the hydrodynamic modeling of natural rivers, was used for one-dimensional unsteady hydrodynamic simulations in this study. The software uses a weighted four-point implicit finite difference method (Gary, 2020) to solve the Saint-Venant equations by discretization and iteration: ∂ A ∂ t + ∂ S ∂ t + ∂ Q ∂ x − q l = 0 (1) ∂ Q ∂ t + ∂ ( v Q ) ∂ x + g A ( ∂ z ∂ x + S f ) = 0 (2) where t is time, x is the longitudinal distance along the river channel, A is the cross-sectional area of the portion of the channel occupied by the flow, Q is river discharge, S is the cross-sectional area of the portion of the channel without water flow, q _l is lateral incoming flow, v is flow velocity, g is the gravitational acceleration, and S _f is friction slope.

The current study focused on the river segment impacted by the Enbridge Line 6 B Oil Discharge, which stretches from Talmadge Creek near Marshall to downstream Morrow Lake. River bathymetry was created by interpolating the field data from the United States Geological Survey (USGS) report (Reneau et al., 2015) with the bathymetry interpolation tool from Merwade et al. (2008). The bathymetry map was then merged with a digital elevation model (DEM, with a resolution of 1 m × 1 m) published by the USGS in order to produce a corrected elevation map of the study area. The constructed main channel model for the HEC-RAS was 64.27 km in length, with an average slope of 0.72‰. The model contained a total of 3893 cross sections with spacing of 1.4–30.9 m apart. Inflow was applied to nine cross sections, the locations of which are illustrated in Figure 1.

Flow data from the 12th of April to the 16th of April 2013 were used to verify the model (Reneau et al., 2015). The model was calibrated using water surface elevation and discharge data from USGS stations 04105500, SR1485, and SR0585 by adjusting Manning’s roughness coefficients. More details can be found in the supplementary data.

In the Lagrangian scheme, the displacement components of each particle in the longitudinal (x), lateral (y), and vertical (z) directions in Cartesian coordinates at time (t + Δt) (i.e., x i t + Δ t , y i t + Δ t , and z i t + Δ t ) were determined by its location at time t, advective movement, and diffusive transport during the time interval Δt (Wu et al., 2019). The Visser scheme was employed to describe the vertical movement of particles (Visser, 1997). The governing equations are as follows. x i t + Δ t = x i t + u i t Δ t + R − 1 1 2 σ R − 1 D x t Δ t (3) y i t + Δ t = y i t + v i t Δ t + R − 1 1 2 σ R − 1 D y t Δ t (4) z i t + Δ t = z i t + w i t Δ t + D z t ′ Δ t + R − 1 1 2 σ R − 1 D z 1 t Δ t (5)

Here i is the number index of a particle, Δt is the time step; u i t and v i t are streamwise and lateral flow velocity components at location ( x i t , y i t ) (m·s⁻¹), respectively; w i t is the average vertical velocity of the particle between t and t + Δt, which is calculated using the Stokes equation (m·s⁻¹); R − 1 1 is a uniform distribution between −1 and 1; σ R is the standard deviation of random number R, which is set at 1/3 (Visser, 1997) and D x t , D y t , and D z t denote the streamwise, lateral, and vertical turbulent diffusivity coefficients at time t, respectively.

In the present study, the transverse flow velocity ( v i t ) was ignored because the longitudinal drifting of the oil particles was more valuable for a fast emergency response. The longitudinal velocity ( u i t ) is assumed to obey the law of the wall (Jones and Garcia, 2018), which is written as: u i t u * = 1 κ ln ( u * z i t ν ) + 5.5 (6) where u * denotes the shear velocity of the flow (m·s⁻¹); ν denotes the dynamic viscosity of water (m²·s⁻¹); and κ denotes the von Karman constant (0.41).

The streamwise and lateral diffusions, D x t and D y t , are given by (Garcia, 2008): D x t = D y t = 0.6 H u * (7)

For D z t , the parabolic-constant profile (Rijn, 1984) was employed and expressed as: D z t = { 0.25 κ u * H z H ≥ 0.5 z 1 ( 1 − z 1 H ) κ u * z H < 0.5 (8) where H is the water depth (m). In the lower half of the depth, D z t is estimated for location z ₁, which is given by z 1 = z i t + D z t ′ Δ t / 2 .

A reflective boundary condition was used when the oil particles exceeded the left, right, and bottom boundaries of the channel. When the particles surfaced, they were assumed to be absorbed into the oil slick. The resuspension probability, p, was calculated to determine whether the particles would reenter the water column (Nordam et al., 2019). p is given by p = 1 − e − ∆ t ⋅ λ (9) where λ is the lifetime of the oil at the surface (Nordam et al., 2019).

In particular, a random depth is assigned if the particle is resuspended; otherwise, the mass loss of the oil droplet due to evaporation should be calculated during the surfaced duration.

Based on a series of documents released by the National Transportation Safety Board, in the present model, the leakage entered the upstream of Talmadge Creek at the cross section located approximately 2.5 km away from the confluence of Talmadge Creek and the Kalamazoo River, and the start time was 18:00 on 25 July 2010. Oil particles were released at a rate of 720 per hour in the calculation domain for 10 h. Based on the oil droplet size distribution, which was reported in the literature based on the breakup of oil in water (Li et al., 2009; Johansen et al., 2015; Zeinstra-Helfrich et al., 2016), nine oil droplet sizes, that is, diameters of 0.5 mm, 2 mm, 4 mm, 6 mm, 8 mm, 10 mm, 12 mm, 16 mm, and 20 mm, were considered. The density of oil at different temperatures is given by (Berry et al., 2012): ρ o i l = ρ r e f [ 1.0 − 8.0 × 10 − 4 ( T − T r e f ) ] (10)

The density of oil after evaporation was expressed by (Berry et al., 2012), ρ o i l = ρ r e f ( 1.0 + 0.018 F e v ) (11) where ρ _oil is the density of the spilled oil (kg·m⁻³); ρ _ref is the reference density at a reference temperature T _ref , which is approximately 920 kg m⁻³ at 20°C; T is the temperature (°C), equal to 23°C; F _ev is the mass loss in percentage terms (%); and t is the duration of evaporation (min).

Based on the measured data obtained by Waterman (Waterman, 2015), a fitting function in logarithmic form (R ² = 0.985) that describes the mass loss due to evaporation was introduced as follows: F e v = 1.961 ⁡ ln ( t ) (12)

The hydraulic data were calculated using the HEC-RAS software using the flow data between July 25th and 31st, 2010 (Reneau et al., 2015). The particle-tracing algorithm was performed using MATLAB software. The computational time step is 2 s.

Two scenarios, one with containment booms and the other with no measures, were conducted. The containment booms were positioned similarly to those in actual situations and extended 0.05 m deep underwater. Figure 2.

When no containment measure was taken, particles with different sizes exhibited similar drifting trajectories in the simulated domain. Figure 3 illustrates the maps of the particle positions when the particle diameter equals 1 mm. The black arrows in the figure indicate the leading edges of the particle swarm at 1 h, 6 h, 12 h, 24 h, 72 h, and 120 h. The particles arrived at Battle Creek at approximately 19.5 h, which roughly matched the time when a resident reported the appearance of an oil slick near Battle Creek, according to documents released by the National Transportation Safety Board.

Figure 4 shows the statistics for the 10th, 25th, 50th, 75th, and 90th percentiles of the distance drifted downstream of the particle swarm, as well as the temporal variation in the average distance. Overall, the polluted zone continued to expand. The 90th and 10th percentiles represent the positions of the front and rear edges of the polluted zone, respectively. The length of the polluted zone in the downstream direction can then be calculated as the difference between the 90th and 10th percentiles. This was calculated as 12.06 km, 22.81 km, 29.61 km, and 34.96 km, at 12 h, 24 h, 36 h, and 48 h, respectively, after the oil entered the water. The length of the polluted zone decreased after 48 h, because the front edge arrived at the downstream boundary of the simulated area.

Plotting the means of the particle swarm’s longitudinal position (X) and vertical position (Z) at different times provides insight into how the center of the particle swarm varies spatially and temporally. As shown in Figure 5, the dimensionless height Z/H was used as well as Z/H = 1.0, which represents the water surface. In general, as the particle size increased from 0.5 mm to 10 mm, the vertical velocity, which is influenced by the buoyancy force, increased. Advective movement driven by buoyancy dominated the vertical position of the particles; thus, more particles (larger than 4 mm) appeared close to the water surface. Small particles were gradually dispersed in water and forced by diffusive movement. After a trip of 10 km downstream, the center of the 0.5 mm particle cloud stabilized at approximately Z/H = 0.6.

After the containment booms were added in the same positions as the actual containment booms that were installed in real response operations, particles with different diameters changed their movement paths. The monitoring sections were set in the model at downstream distances of 10 km, 30 km, and 50 km, as shown in Figure 3. As shown in Figure 6A, the number of particles passing through each section was determined. The containment booms implanted before the monitoring section at a distance of 10 km exhibited a slight interception effect on particles smaller than 10 mm. Approximately 26.12% of the 16 mm particles and 65.93% of the 20 mm particles were contained in the upstream booms. At a distance of 30 km, the containment rates for the 8 mm, 10 mm, and 20 mm particles increased to 7.57%, 19.72%, and 91.23%, respectively. At a distance of 50 km, the corresponding containment rates further increased to 27.54%, 52.55%, and 95.11%, respectively. Unlike large oil droplets, which migrated in the upper part of the water column, small oil droplets were hardly blocked by the implanted booms as these particles were suspended over depth, as revealed in the previous section. The containment effectiveness was enhanced by increasing the number of containment booms. However, their relationship is complex because it is also affected by the actual particle properties and hydraulic conditions.

Statistical analysis was performed for the time when the particles arrived at each monitoring cross section. In the model with no containment measure implanted, the mean time for particles to reach the sections 10 km, 30 km, and 50 km downstream was 11.80 ± 0.41, 35.15 ± 1.34, and 55.01 ± 1.84 h, respectively. The average arrival times in the presence of containment booms are plotted in Figure 6B. We found that the diffusion of oil droplets greater than 6 mm was slowed by the containment booms.