On the application of the depth-averaged random walk method to solute transport simulations

Pollutant transport in waters has been extensively studied since it is relevant to a wide range of environmental problems and may cause significant economic losses (Rajar 1997). Among the natural waters, many can be regarded to be shallow, such as estuaries, lakes, reservoirs and nearshore waters. The horizontal scale of these flows is much larger than their depth. Therefore, it can be assumed that the solute is well-mixed vertically over the water column, allowing a depth-integrated approach to model the solute transport process. This paper focuses on the solute transport in shallow waters.

With the development of computing techniques, many modern numerical methods have been established to solve such problems (Lin & Falconer 1997; Gupta et al. 2004; Begnudelli & Sanders 2006; Benkhaldoun et al. 2007). Previously, most research has been based on the Eulerian approach to solve the advection-diffusion equation using finite-difference, finite-volume or finite-element method (Mingham et al. 2001; Liang et al. 2006; Yang et al. 2018). However, these grid-based methods have proven to be problematic in addressing steep concentration gradients and tend to produce numerical diffusion (Wu & Liang 2018). By contrast, Lagrangian methods use a group of discrete particles to represent the pollutant cloud, and these particles are tracked independently in flows. The Lagrangian approach can usually yield a good estimation regardless of the concentration gradients (Zhang & Chen 2007; Wu et al. 2018), and the computation can be easily parallelised. Nowadays, thanks to the rapid progress in electronic technology, the high computational cost of the Lagrangian approach has become increasingly affordable. Its merits include inherently conserving mass and being free from artificial diffusion.

The computational cost of the Lagrangian approach depends on two factors: the number of particles to represent the pollutant cloud and the size of the computational time step. The choice of these two factors is crucial to the efficiency of random walk modelling. Traditionally, it is necessary to manually tune the number of particles and the size of time steps. Such a trial-and-error procedure is time-consuming and relies on the modeller's experience. This paper first uses an ideal case study to investigate the sensitivity of numerical solutions to the time step and particle number in depth-averaged random walk modelling. Then, the model is used to solve a solute oscillation problem along a hypothetical tidal estuary. The results of a TVD–MacCormack model are used as a reference. Finally, the paper describes the application of the random walk scheme to simulate the wind-driven chaotic mixing in a shallow circular lake. The results are qualitatively compared with those in previous research, and the importance of diffusion is highlighted.

The source term q_s in the previous equation is neglected in the new form, as all the following case studies do not include any source or sink. According to the stochastic process theory, S can considered as a probability density function. U and V in Equations (8) and (9) represent the modified advective velocity components. The more detailed explanations of this scheme can be found in Wu & Liang (2018). This depth-averaged random walk model is then implemented by performing the advection and diffusion transport processes in each time step, as explained in the following sections.

Particles follow a ‘random walk’ trajectory with the flow inside the computational domain. The boundary treatment is designed to prevent them from crossing solid boundaries. In this study, the fully reflective boundary condition is applied to particles that penetrate the solid boundaries after implementing Equations (18) and (19). When the time step is not too large, such a treatment then reflects the particles back into the computational domain.

In this case study, the total amount of solute material M = 233.06 kg is released suddenly at the initial location (x₀, y₀). As shown in Figure 1, the water depth is set to be h = 1 m, and the Chézy coefficient is 40 m^1/2/s over the whole study area. The streamwise dispersion and transverse diffusion coefficients are typical values of 13.0 and 1.2, respectively, in straight open channel flows (Falconer 1991). Two flow conditions are considered in this section. The first one is a uniform flow with u = 1 m/s along the x-axis (θ = 0). The solute material is initially located at (x₀, y₀) = (0, 400 m). The second one is a uniform flow aligned diagonal direction (θ = 45°). The velocity is set to be u = v = 1/m/s. For computations in both scenarios, the particle numbers P is 2.33 × 10⁶ and the time step is 1 s.

The development of concentration contours for the flows along the x-axis direction and diagonal direction are presented in Figure 1(a) and 1(b), respectively. It is notable that the major axis of the elliptical patches is along the flow direction. The cloud of the solute experiences rapid elongation along its path. The reason is that streamwise dispersion is over ten times larger than the transverse diffusion. The results for the same test case obtained by a grid-based method, TVD–MacCormack, can be found in Liang et al. (2010).

The advantage of the random walk model includes high accuracy and small numerical diffusion. However, the disadvantage that is associated with this Lagrangian approach is its high computational cost. As seen in the scheme description, the amount of the computation depends on two factors: the size of the computational time step and the number of particles to present the pollutant cloud. Therefore, the choices of these two factors are crucial to the efficiency of the random walk model. The following part of this section is to discuss the influence of the two parameters on the random walk model when applied to the instantaneous release problem in uniform flows.

A great deal of research on time steps has been undertaken for the Euler scheme. The size of is restricted by the Courant–Friedrichs–Lewy condition. Usually, the smaller the time step is, the more stable the simulation will be. However, the behaviour of the random walk model is very different. Figure 2 presents the longitudinal concentration profile at 600 s after the solute release. Regardless of the time step of 600 s or 0.01 s, the same concentration profile is predicted by the model. Variations of the peak concentration are not affected by the size of time steps as well, as seen in Figure 3.

Therefore, the results will obey the distribution of no matter what size of the time step is used. It can be concluded that the change of the time step does not affect the accuracy of the computation for this case of uniform flow with constant water depth.

Figure 4 shows a qualitative illustration of the solute transport process predicted by the random walk model using different numbers of particles. In general, the contour turns out to be more reasonable with the particle number increases. Larger particle numbers significantly improve the visualisation of the elliptical cloud for solute distribution.

To get a quantitative analysis, the predicted results are compared with analytical solutions. Figure 5 shows variations of peak concentration with different particle numbers used in the simulation. The size of sampling bins is set to be a circle with a radius of 1 m. It is apparent that small particle numbers tend to produce numerical oscillations. On the contrary, by introducing a larger number of particles in the model, the prediction further approaches the analytical solution. For example, when 2.33 × 10⁷ particles are introduced in the model, the peak concentration is identical to the analytical solution. The relative errors with different particle numbers and different sampling bin sizes are compared in Figure 6(a). In the legend, r represents the radius of the bin. It is clear that, with the increase of particle numbers involved in the simulation, the error decreases significantly. The error for insufficiently small particle numbers can jump beyond 18%, which is unacceptable in most simulations. When the number of particles in each sampling bin approaches 200 to 300, the relative error is reduced to less than 5%. This means that the particle should be guaranteed to be higher than hundreds per bin to give a reasonable estimation of the solute concentration. In addition, it is worth noting that the error for both x-direction and diagonal flows is always very similar, as shown in Figure 6(b). This means the random-walk model is not affected by the flow direction, unlike most of the mesh-based models.

This section considers the situation where the solute is transported forward and backward in a hypothetical tidal estuary. As shown in Figure 7, the boundary condition for the left end is regarded as a sinusoidal tide, while the right side of this estuary is considered as a vertical wall. The total length of the estuary is 13,800 m, and the bed level changes from −5 m at the end of the seaward side to 0 m at the right end. In this case, the tidal flow has an average water level of 0 m and amplitude of 2 m, rising from the average water level at the beginning of the simulation. Before the pollutant transport simulation, the predictions of the flow field were obtained by using TVD–MacCormack scheme to solve the SWEs (shallow water equations). Detailed information can be found in Liang et al. (2010). The discretised velocity field is then reconstructed over the whole domain using linear interpolation.

The initial concentration is set to be 100 units at the cell located at x = 10 km, while the concentration at the other positions is set to zero. An illustration of the particles' position and their distribution at different times is given by Figure 8. It is well presented that particles oscillate with tidal flows and gradually disperse in the longitudinal direction. The instantaneously released solute moves with the rising tide to the closed wall under the influence of tidal currents, while the receding tide carries the particles back to the open boundary. At the same time, the spread of the solute grows because of the effect of dispersion.

The aforementioned sampling algorithm is required to convert the scatter of particles, as illustrated in Figure 8, to the concentration profile, as illustrated in Figure 9. This sampling algorithm is crucial in the interpretation of the results of the random walk method and in the evaluation of the concentrations. Traditionally, this concentration profile, i.e., histogram, uses sampling bins of fixed length. However, this fixed length may not suit all the positions. In this work, the length of the bin is chosen automatically by including 200 particles within each bin. This dynamic determination of the bin size is found to be able to avoid spurious fluctuations and achieve good comparisons with analytical and previous results. A series of smooth concentration profiles along the estuary at both high and low water levels are shown in Figure 9. When the solute cloud moves close to the landward position, the tidal flow generally decelerates because of the wall boundary. It is notable that the peak concentration at 27 hours is even higher than that at 9 hours, as the water depth is smaller near the right-side wall which amplifies the value of the concentration to conserve mass of the solute.

The numerical solution of the TVD–MacCormack model is used as a reference for this case study. For such a grid-based method, solutions are sensitive to the grid size. The finer the grid is, the more accurate the simulation will be. The numerical diffusion for the Eulerian model is most obvious at the early stage of the simulation because of the relatively sharp concentration gradient. Take Figure 10(a) as an example, which presents the concentration profile at 3 hours. Several different cell sizes, from 3 m to 270 m, are used in the mesh-based simulations. When the mesh size is increased to 270 m, the concentration distribution is far too flattened because of large numerical diffusion. The peak concentration for the grid size of 3 m is nearly three times higher than that for the grid size of 90 m. However, very small grid size leads to very large computational cost. On the contrary, the random walk model does not involve numerical diffusion. A narrower distribution and a higher peak concentration are predicted by the present model, and it is similar to the prediction of the mesh-based method with very small grid size.

These physical parameters are set to be the same as those in Kranenburg (1992) for comparison purposes.

For the purpose of comparison with previous studies, the diffusive process is first ignored so as to test the behaviour of the advection part of the random walk model. These particles are tracked for 500 periods. For small μ, as shown in Figure 12(a), the superposed particles with the same phase form elliptical patterns that occupy the left and right halves of the lake. The area of the chaotic region increases with the increase of μ, while the sizes of the noticeable elliptical patterns decrease. When μ is beyond 0.70, the regular patterns are no longer obvious, and only chaotic distribution is noticed (for μ = 0.84). These findings are consistent with previous results of Kranenburg (1992).

Then, both advection and diffusion processes are considered as tracer particles are also spread by turbulence diffusion and shear dispersion in reality. Figure 13 depicts the impact of streamwise dispersion (the transverse diffusion is set to 0 in this case study) on this chaotic mixing phenomenon. For equal to 0.001, it is still clear that the circular lake is divided into two halves and patterns of the elliptic movement are still obvious. With the increase of the streamwise dispersion coefficient, this pattern becomes more and more blurred until the KAM (Kolmogorov–Arnold–Moser) (Kranenburg 1992) curves no longer exist (⁠ > 0.1). It is worth mentioning that the growth of the streamwise dispersion will not change the mean position of points, but it will remove the large local concentration gradients by increasing the random motion.

For μ = 0.28 and 0.42, the entire calculation domain is chaotic when the dispersion is as small as 0.1, as shown in Figures 14 and 15. This means that the increase in the duration of the periodic storm will emphasise the influence of the diffusion coefficient, making it easier to achieve a chaotic state. Figure 16 shows the advection and diffusion behaviour of a line of 10,000 particles after 32 periods of the storm event. The line is initially positioned along the x-axis. The whorl and tendril structures coexist when μ is 0.28. With the increase of the dispersion parameter, these fine structures both disappear and particles are spead anywhere in the lake, although it can still be seen that the particle distributions in the left and right halves are symmetrical.

The traditional random walk model has been extended to solve the depth-integrated advection-diffusion equation. First, this model is verified by solving an instantaneous release problem in a uniform flow. The analytical solution is used as a reference. The results reveal several advantages of this model, including high accuracy, stability and simplicity. Extensive parametric studies have been carried out to investigate the sensitivity of the predictions to the computational parameters. It has been found that simulations are independent of the size of the time step. The particle number significantly influences the performance of the random walk model. Too few particles degrade the visual inspection and make it impossible to quantitatively examine the solute distribution. In uniform-flow applications, a relatively large time step will reduce the computation expense without compromising accuracy. Particle numbers should be chosen to ensure at least 200–300 particles in each sampling bin. Subsequently, investigations are carried out regarding the oscillation of a pollutant cloud in a tidal estuary. The sampling method is optimised to convert the particle distributions into concentration profiles. The random walk simulations display high accuracy, which can only be achieved by the mesh-based simulations with extremely fine resolutions. The mesh-based methods are shown to be highly sensitive to grid resolution. Finally, the model is used to simulate the chaotic mixing process in a circular lake. The results for the case of pure advection are consistent with the findings reported by previous researchers. Because of the presence of turbulent diffusion and bottom friction in any real-world lake, the advection is always accompanied by transverse diffusion and longitudinal dispersion. This study shows that diffusion plays an important role in the mixing. The well-structured chaotic mixing pattern produced by pure advection cannot be observed in reality.

In summary, this study demonstrates that the random-walk model is highly stable and free of artificial diffusion in solving the solute transport problems in aquatic environments. In particular, this study proposes a generic sampling technique to convert the scatter of discrete particles to the distribution of solute concentration. In order to obtain a good estimate of the solute concentration, this study proposes a guideline for estimating the total number of particles needed in the simulation.

We are grateful for financial support provided by the Royal Academy of Engineering UK-China Urban Flooding Research Impact Programme (No. UUFRIP\100051), the 111 Project (No. B17015) and the China Scholarship Council (No. 201708060090).