Turbulent Mixing of Floating Pollutants at the Surface of the River

Based on the concept of molecular diffusion, Taylor (1921) first introduced Lagrangian statistics to analyze the scattering of particles attributable to turbulent motion. Following Taylor’s study, after a sufficient period of time ($t\gg T_I$), the variance of the ensemble-averaged concentration distribution for clouds dispersing in a stationary homogeneous field of turbulence grows linearly with time aswhere $⟨x_i^2⟩$ = ensemble mean size of the cloud, in which $i$ denotes the streamwise and spanwise coordinates; $u_i^{\prime }=u_i-{\overline{u}}_i$ = turbulent fluctuation that becomes the velocity of the particle, with zero mean (${\overline{u}}_i=0$); $⟨{u^{\prime }}_i^2⟩=⟨u_i^{\prime }(0)u_i^{\prime }(0)⟩$ = square of the turbulence intensity; $T_I$ = Lagrangian integral timescale, and this timescale is the time that particles remember their own velocity in the homogeneous and stationary turbulence field; $T_I$ = integral of the Lagrangian autocorrelation function ($R_i$) with given timeand $R_i$ is defined as

When $t$ is significantly large, $R_i$ decreases from 1 to 0 because the motions of the particle gradually lose their correlation with time $t$.

In this study, following the moment method (Fischer et al. 1979), Eq. (1) was rearranged to calculate the diffusion coefficient aswhere $K_i$ = diffusion coefficient in streamwise and spanwise directions, and ${\sigma }_i^2$ = variance of the particle cloud. Eq. (4) shows that $K_i$ can be obtained from the time evolution of variance with the time series of particle distributions (Rutherford 1994). Richardson and Stommel (1948) reported that the relative diffusion between two particles was measured by the ensemble mean of separating distance. Further, using the fact that the ensemble mean is related to the size of the particle cloud, Rowinski et al. (2005) calculated the average size of the particle cloud from the position data of particles. Thus, in this study, from the discrete particle distribution, variance in the $i$-direction was calculated using the following equation (Rowinski et al. 2005):where ${\overline{x}}_i=1/N\sum _{k=1}^N{x}_i^k$ = centroid of particle distribution at time $t$; $\overline{()}$ indicates the average value; and $N$ = number of particles. When particles introduced at the water surface are translated and scattered by turbulent motion, ${\sigma }_i$ in Eq. (5) is the average displacement of the scattered particles, and this value represents the size of the particle cloud. In this study, a number of GPS-equipped floaters were deployed onto the surface of the river, and ${\sigma }_i^2$ was then calculated from the time series data of the floaters’ displacements.

Mixing of floating pollutants in turbulent flows can be modeled using either an Eulerian or Lagrangian description in a fluid continuum. In the Eulerian description, pollutant mixing is represented in a fixed control volume at a specific time $t$ and is described using the advection-diffusion equation. The advection-diffusion equation for the turbulent flow in the Eulerian description can be derived from Eq. (6) by incorporating Fick’s law into the mass conservation equation (Fischer et al. 1979)where $C$ = time-averaged concentration; $t$ = time; $x_i$ = Cartesian coordinates in three dimensions; ${\overline{u}}_i$ = time-averaged velocity components; and $K_{ij}$ = turbulent diffusion coefficient. Since grid-based models based on Eq. (6) suffer usually from artificial oscillation and/or numerical dispersion, particle tracking methods in which pollutant particles are traced using the Lagrangian approach have been developed to predict the spatial mixing of particles and the temporal concentration field (Lonin 2000). The particle tracking model was initially developed in the 1950s to investigate the contaminant transport in groundwater flows (Scheidegger 1955). The Lagrangian description is more suitable than the Eulerian description for tracing the trajectories of discrete pollutant particles. The Lagrangian approach is a grid-free numerical method, so it can avoid artificial oscillations and numerical dispersion and can save computational cost. Therefore, a number of particle tracking models have been developed to simulate the contaminant transport in turbulent flows of inland and coastal waters. In conjunction with the 3D hydrodynamic model, particle tracking models such as the Lagrangian transport model (LTRANS) (Schlag and North 2012), the oil transport model (OILTRANS) (Berry et al. 2012), and the LPT model of EFDC (Dunsbergen and Stalling 1993) were developed to simulate the mixing of drifting particles, sediment, microorganisms, and oil droplets, primarily in coastal waters. Even though the Lagrangian approach is an effective method to simulate the mixing of floating pollutants, the particle tracking method has some limitations; a sufficient number of particles needs to be introduced for the random fluctuations, and some numerical difficulties are found in the results in a highly distorted grid (Salamon et al. 2006).

In the Lagrangian models, the transport of pollutant particles in the water body is described by a combination of the deterministic translation and random motion (Itō and Nisio 1964). If the position of particles, $x_i$, has a conditional probability density function, then the Fokker-Planck equation was defined as (Gardiner 1985)where $p=p[x_i,t|x_i(t_0),t_0]$ = conditional probability density function; $A_i$ = arbitrary vector to transport particles with deterministic forces; and $B_{ij}$ = arbitrary tensor that drives particles to random motion. Comparing Eqs. (6) and (7), which are mathematically similar, the following equation can be derived:Eq. (8) includes the drift term, which covers the indeterministic transport of the simple random walk model from the Wiener process (Monti and Leuzzi 2010).

In this study, as a particle tracking model, the LPT module of EFDC was applied. The LPT model uses Eq. (8) as the governing equation in Lagrangian coordinates (Dunsbergen and Stalling 1993). With the isotropic turbulent diffusion in the horizontal direction, Eq. (8) can be decomposed intowhere $u$, $v$, and $w$ = velocity components in $x$-, $y$-, and $z$-directions, respectively; $K_H$ = horizontal turbulent diffusion coefficient; $K_V$ = vertical turbulent diffusion coefficient; and $Z$ = random number that follows the uniform distribution with a mean quantity of 0 and a range from $-1$ to 1. In this study, for the analysis of pollutant mixing at the surface of natural rivers, vertical displacements were neglected, and only two horizontal transports were considered. In this study, for the flow simulation, a quasi-3D hydrodynamic model of EFDC was applied, in which a hydrostatic assumption was used in the vertical direction.

The field measurements were conducted in the middle reach of the Nakdong River in Korea as shown in Fig. 1. In this river, to control the water discharge and the water level of the river, eight movable weirs were recently constructed in 2012. In this study, two series of experiments were conducted: Series 1 at Site 1, located between the Gangjeong-Goryeong weir (GGW) and the Dalsung weir (DW) in 2012 and 2013, and Series 2 at Site 2, located between the Gumi weir (GW) and the Chilgok weir (CW) in 2013. During the experiments, the water levels at Sites 1 and 2 were maintained as fixed values of 14.0 and 25.5 m, respectively. In these reaches of the river, the flow rate was manipulated by operating the gates of the movable weir to maintain the water level. At Site 1, the Geumho River merges with the Nakdong River downstream of GGW. Two large industrial complexes are located on the left side of the Geumho River and the right side of the Nakdong River. At Site 2, the Gam Creek and the Han Creek merge with the Nakdong River, and industrial complexes have been constructed near the tributaries. Downstream of the confluence in Sites 1 and 2, many water intake facilities operate to draw water from the river for municipal and agricultural purposes. Since water pollution accidents by oil spill are highly possible in this reach, the mixing behavior of floating pollutant needs to be studied based on tracer experiments and numerical simulations.

At Site 1, three tests, Cases GF11, GF12, and GF13, were conducted in the upper, middle, and lower portions of the reach, respectively, as shown in Fig. 1. In Case GF11, the floaters were released at the downstream of the Samunjin Bridge near the Geumho River confluence. The test for Case GF12 was conducted near the Dasan industrial complex, and that for Case GF13 was conducted at the downstream of the Seongsan Bridge, which is located upstream of the DW. At Site 2, the experiments were conducted near the confluence of tributaries in which water pollution accidents had occurred previously. Case GF21 was located near the confluence of the Gam Creek, which is 1.5 km downstream of the GW. At 13 km downstream from the Case GF21 site, and downstream of the confluence of Han Creek, the floaters were deployed on the left bank for Case GF22 and on the right bank for Case GF23.

In this study, floaters consisting of a plastic bulb in which the GPS sensor was inserted were used as the tracer for Lagrangian particle tracking experiments. As shown in Fig. 2, the GPS floater has a spherical shape to minimize the disturbance around the floater and the pseudodiffusion induced by the inertial force and wind stress. The diameter of the floater ($d$) was 10 cm in order to embrace the 7.2-cm-long GPS in which a self-logger was embedded to store time-series position data without additional equipment. The submerged depth of the floater was set as 4.5 cm with a weight of 145 g, in order for it to drift by the surface currents and turbulent motions. The GPS was situated at the center of the floater to receive stable data. The accuracy measurements of the GPS are usually defined with a dilution of precision (DOP) (Langley 1999), and the GPS in this study has 0.95 horizontal DOP (HDOP) and 2.11 vertical DOP (VDOP). Han (2016) maintained that this specially fabricated floater, named the GPS floater, could adequately represent the behavior of the floating pollutant; the laboratory experimental results of Han (2016) showed that the difference between the water velocity and transport velocity of a GPS floater was in range of 0.84–4.23%. This floater is neutrally buoyant, and its size is small enough to cover the spatial scale of the smallest turbulence at the surface of the river. Furthermore, unlike the image analysis required with the particle image velocimetry (PIV) and particle tracking velocimetry (PTV) methods (Suara et al. 2015), the data of particle position versus time were easily acquired without any further processes by tracking the position of the floater using the embedded GPS when drifting in river flows. In this study, the position data of each GPS floater along with time data were stored every 3 s. Before deploying GPS floaters, preliminary experiments were conducted with a dummy floater, which had the same size and weight as the GPS floater, to find a consistent trajectory in each trial. In the main experiments, 25–35 GPS floaters were deployed to minimize the interference between floaters that causes false spreading on the water surface. The logged data from the GPS floaters were stored with the format of World Geodetic System 84 (WGS84), a coincident coordinate system that is necessary to calculate the number of particles in a particular computational grid. The computational grid was defined in transverse Mercator (TM) coordinates, and the position data from the GPS floaters were transformed to TM coordinates (National Geographic Information Institute 2005). From the coincident coordinates system, the horizontal diffusion of the GPS floaters and LPT simulation results were able to be compared, and the floating pollutant mixing was analyzed.

In this study, both the GPS floater experiments and LPT simulation results produced discrete particle position information. Even though these discrete particle positions are useful to demonstrate the floating pollutant mixing, this information is not suitable to anticipate the peak values of the pollutant concentration and distribution of pollutant clouds. Therefore, in this study, the particle position was converted to the concentration field using the concentration conversion technique given below. Since the particle density is defined as the number of particles in a unit volume (Dimou and Adams 1993; Israelsson et al. 2006), the concentration value was calculated at a point by using Eq. (10) (Suh 2006)where $C(x,y,t)$ = concentration field; $m$ = mass contained in each particle; $n_p(x,y,t)$ = number of particles in a particular computational grid; $H$ = water depth; and $\mathrm{\Delta }x$ and $\mathrm{\Delta }y$ = grid size in the $x$-and $y$-directions, respectively.

The river velocity data were collected using the acoustic Doppler current profiler (ADCP) at Section 1 of each site shown in Figs. 3 and 4, and the cross-sectional average values of velocity and depth from the ADCP measurements along with wind speed data were listed in Table 1. The secondary current profiles and velocity magnitude contours from the ADCP transects at Section 1 of each case were plotted in Fig. 5 using the VMT software (Parsons et al. 2013). From the velocity measurements in Fig. 5, the secondary current intensity (SCI) was calculated using Eq. (11) (Seo et al. 2006) and listed in Table 1where $u_n^{\prime }=u_n-{\overline{u}}_n$, $u_n$ = spanwise velocity at vertical measurement points; ${\overline{u}}_n$ = depth-averaged transverse velocity; $U_p$ = cross-sectional averaged streamwise velocity; and $L$ = number of lateral measurement points. As shown in Table 1, even though the SCI of Series GF1 was larger than those of Series GF2, it is difficult to directly relate the value of the SCI with the turbulent diffusion by surface currents as shown in Table 2 because the secondary currents that were induced by the channel curvature and irregularities of the channel boundary contributed to an increase of the transverse dispersion rather than the turbulent diffusion (Rutherford 1994; Sharma and Ahmad 2014).

The surface current can be evaluated by the surface velocity and wind stress. In this study, the surface velocity was taken from the top cell of the ADCP measurements. The cross-sectional mean values of surface velocity, $U_S$, and the magnitude of spanwise component, $U_{ST}$, were given in Table 1. As shown in this table, as with the SCI, the $U_{ST}$ of Series GF1 were larger than those of Series GF2. Also, in this table, wind speed data were given to consider the wind stress effect on the floaters, in which the local wind speed and direction were measured with an anemovane, which has 3% accuracy. This data showed that the wind speed was lower than $2.0\mathrm{m}/\mathrm{s}$; thus, in this study, the force induced by wind on the floater was considered to be insignificant.

In Figs. 3 and 4, the trajectories of all GPS floaters for each case are depicted as continuous path lines. These lines can be considered as streak lines because GPS floaters were simultaneously introduced at the same injection point as that in the test reach. However, these lines do not represent the streamlines because the particles had drifted as a result of both the mean motion and the unsteady turbulent fluctuations of the river water. In Case GF11, the GPS floaters released downstream of the Samunjin Bridge were transported downstream by about 716 m for 85 min; the floaters then reached the left bank of the Nakdong River, where they became trapped, as shown in Fig. 3(a). This figure shows that for Case GF11, the scattering of the path lines in the spanwise direction was relatively large compared to those of the other cases. This was considered to be a result of the stronger lateral mixing induced by the velocity fluctuation compared to the advection caused by the mean motion at the water surface. The GPS floaters were transported down the river by about 414 m for 20 min in Case GF12 and by about 500 m for 26 min in Case GF13.

Compared to the experimental results of Site 1, in the tests of Site 2 the floaters were transported a longer distance downstream because of the higher velocity at this reach, as shown in Fig. 4. In Case GF21, the floaters were transported about 1.6 km for 79 min, and landed on the left bank since surface currents at the curved bend were directed outward because of the secondary currents, as shown in Fig. 5. In Cases GF22 and GF23, the floaters deployed downstream of the Sanho Bridge were transported by about 768 m for 66 min in GF22 and by about 556 m for 53 min in GF23.

The measured data of the floater scatterings are plotted in Figs. 6 and 7 along with simulation results by EFDC. In this study, the floater distributions were plotted as a scatter diagram at the specific times. These figures show that the floaters were transported by the surface current while being simultaneously scattered by the turbulent fluctuations, which can be represented by the velocity deviation on the water surface in both the streamwise and spanwise directions. In this study, the advection of the particles was represented by the translation of the centroid of the cloud, and the turbulent diffusion, as mentioned in the previous section, was calculated by the scattering of each individual particle from the centroid in both the streamwise and spanwise directions.

These figures clearly demonstrate that particle scattering in the flow direction was considerably greater than the spanwise scattering. This variation in the scattering is considered to be a result of the stronger turbulent motion in the streamwise flows of the natural rivers, in which the streamwise current is also much higher than the spanwise flow. Thus, as mentioned in the previous section, the turbulent diffusions at the surface of the rivers in the streamwise direction and in the spanwise direction should be determined separately. However, in the Lagrangian particle tracking models (e.g., EFDC, OILTRANS, and LTRANS) in which there are no directional differences, only the single turbulent diffusion coefficient in the horizontal direction has usually been considered.

The numerical simulation of the particle tracking using the LPT model was conducted, and the results were compared with the experiments in Site 1. Simulation conditions were the same as the field experiments as given in Table 1. The vertical layer in the LPT model was divided into 10 layers to obtain the surface velocity, and the horizontal grid size was $7\times 10\mathrm{m}$ on average for reproducing the flow feature on the water surface. The same number of particles was introduced at the injection point where the GPS floaters were released. In the LPT simulation, the released particles were only transported at the surface layer to enable the use of the surface velocity in determining the trajectory of the floating pollutants, as mentioned in the previous section.

For the validation of the numerical model, the flow simulation results of EFDC were compared with the ADCP measurements at Section 1 as shown in Fig. 10. From the results, the EFDC simulation results show a similar surface velocity distribution in the transverse direction as that in the field measurements. The observed values of the diffusion coefficients in Table 2 were applied to the LPT simulation, and the horizontal mixing of floating particles was then compared with the experimental results. In analyzing the field measurements, directional diffusivity was considered to calculate the diffusion coefficients. However, in the LPT model in EFDC, isotropic mixing is assumed; thus, in this study, $K_H$, which was the geometric mean value of $K_L$ and $K_T$, was used in the LPT simulation.

In Figs. 6 and 7, the particle scatterings in the simulation results were plotted on the $x-y$ plane, along with the observed positions of particles from the field experiments. Floating particles from the simulation results show trajectories similar to those of the experimental results, even though the particle trajectories during early periods had shown some discrepancies. The percent error of particle displacements in Case GF11 was 28.2% after 30 min, which decreased to 10.7% after 50 min. The error also decreased from 17.4 to 14.7% in Case GF22. Thus, as aforementioned, the displacements of the floaters in the initial mixing phase were restricted because of the interaction between the floaters and the external forces induced by the limited size and weight of the floaters. In some cases, especially Cases GF12 and GF21, some floaters in the observed data deviated from the simulated cloud. It is considered that this deviation was caused by local eddies that have different length scales from those simulated by the EFDC model.

Particle distributions were converted to the concentration field using Eq. (10). Fig. 11 shows a comparison of the observed and simulated concentration contours. As aforementioned, the location of peak concentration by the LPT model differs somewhat from the observed results, especially for Case GF11. For this reason, it is suggested that the numerical model with $K_L$ and $K_T$ be applied to the analysis of the turbulent mixing in the river system, where the current and the turbulence are stronger in the streamwise direction than in the spanwise direction.

In the case of pollutant spills, decision makers need detailed information about the peak concentration and arrival time of the contaminant cloud in order to set up a disaster prevention plan for a prompt response to the water pollution accident. For this purpose, time variations of concentration curves were extracted from the concentration fields shown in Fig. 11. In addition, the $C\text{-}t\text{-}y$ curves at Section 2 for GF11 and GF13 were plotted in Fig. 12 for simulation and experimental results. In Fig. 12, $C_0$ is the injected particle concentration and $C_p$ is the peak concentration at each section. In Case GF11, the pollutant cloud was transported along the left bank; according to the decreased velocity distribution, the detention time of the particle cloud at Section 2 was about 10 min. For the pollutant cloud in Case GF13, the peak concentration occurred at $y/W=0.3$ after 15 min, and the detention time was about 3 min.