Comparisons of Two Types of Particle Tracking Models Including the E ﬀ ects of Vertical Velocity Shear

: In this study, two types of particle tracking models were presented to investigate the applicability in the two-dimensional solute mixing simulations. The conventional particle tracking model, denoted as PTM, was developed based on Fick’s law, which adopted the dispersion coe ﬃ cient to calculate the random displacements. The other model is the particle dispersion model (PDM), which computes the shear dispersion process by dividing into two computation procedures as the shear translation and the vertical mixing. The PTM and the PDM included the e ﬀ ects of vertical proﬁles of velocity in the computation of dispersion coe ﬃ cients and the shear translation step, respectively. The main di ﬀ erence between the two models is whether the shear dispersion process is reproduced using Fick’s law or the direct computation method. These di ﬀ erences were clearly revealed by comparing with the analytic solution of the advection-dispersion equation. The concentration curve resulting from the PTM shows the Gaussian curves, which were well-ﬁtted with the analytic solution in both initial and Taylor periods. Meanwhile, the PDM presented skewed curves in the initial period and gradually turned to the symmetric shape in the Taylor period. The inherent di ﬀ erences of the two particle tracking models were scrutinized against the two-dimensional tracer test results, which show the non-Fickian mixing properties. The comparisons of concentration–time curves reveal that the PDM reproduced a more accurate shape of the curves than the results by the PTM by demonstrating skewed concentration curves.


Introduction
The increase of urban and industrial areas become a threat to aquatic environments by increasing possibilities of water pollution. Accurate predictions of pollutant mixing in rivers can reduce damage from water pollution accidents and secure public health by drawing essential information for disaster mitigation activities. The analysis of soluble pollutant mixing introduced in rivers is mainly conducted using two-dimensional (2D) models, as the vertical mixing is completed during a relatively short period [1,2]. The 2D mixing can be described based on the shear dispersion theory, as described in Figure 1, in which a pollutant cloud is stretched by shear flows in the longitudinal and transverse directions, and simultaneously diffused through flow depth by the turbulent diffusion [3]. In the 2D mixing model, the shear dispersion theory is modeled using Taylor's study [4], in which the mass flux where h is the flow depth; z is the vertical direction; u i = u i − u i and c = c − C are the vertical fluctuations, which are subtractions of the depth-averaged values (u i , C) from the vertical distributions (u i , c), of velocity and concentration, respectively; and D ij are the dispersion coefficients in tensor notation. The relation shown in Equation (1) is only valid after the initial period, defined as below.
where t I is the initial period, h is flow depth, and ε z is the vertical turbulent diffusion coefficient. Thus, in the 2D mixing simulations, determination of dispersion coefficients is one of the crucial concerns to conduct pollutant mixing simulations.
Water 2020, 12, x FOR PEER REVIEW 2 of 18 the mass flux term derived from the depth-averaging process is assumed to be proportional to the concentration gradient based on Taylor's assumption as follows [5].  (2) where I t is the initial period, h is flow depth, and  z is the vertical turbulent diffusion coefficient. Thus, in the 2D mixing simulations, determination of dispersion coefficients is one of the crucial concerns to conduct pollutant mixing simulations. Several studies adopted empirical formulas obtained from tracer test results or theoretical approaches to determine the dispersion coefficient for solute transport simulations [6][7][8]. The empirical formulas of dispersion coefficient are presented as follows. properties of rivers [9][10][11][12][13][14]. On the other hand, the theoretical formula of L D suggested by Elder [15] has been widely used for simulations instead of using the empirical formulas [2,16]. The theoretical value of L a is derived as 5.93 under the assumption of an infinitely wide width of straight channel and the log distribution in the longitudinal velocity shear. However, the theoretical value shows large discrepancies from the results of field studies, which suggested the value of L a to be in the range of Several studies adopted empirical formulas obtained from tracer test results or theoretical approaches to determine the dispersion coefficient for solute transport simulations [6][7][8]. The empirical formulas of dispersion coefficient are presented as follows.
where D L and D T are the longitudinal and transverse dispersion coefficient, respectively; a L and a T are the empirical coefficients; and u * = ghS f is the shear velocity, in which S f is the energy slope, which can be obtained from Manning's formula. In case of D T , several studies have been reported that a T has a wide range of 0.15-3.0 according to hydraulic and geographical properties of rivers [9][10][11][12][13][14].
On the other hand, the theoretical formula of D L suggested by Elder [15] has been widely used for simulations instead of using the empirical formulas [2,16]. The theoretical value of a L is derived as 5.93 under the assumption of an infinitely wide width of straight channel and the log distribution in the longitudinal velocity shear. However, the theoretical value shows large discrepancies from the results of field studies, which suggested the value of a L to be in the range of 13-150 [17][18][19]. Instead of the aforementioned dispersion coefficient formulas, theoretical approaches reported by Fischer et al. [5] have been adopted to obtain dispersion coefficients by conducting triple integrations of vertical profiles of velocity [20]. The velocity-based method to calculate the dispersion coefficient is reported to adequately reproduce solute mixing in meandering channels [21]. Despite several studies to define the dispersion coefficients, it is still difficult to employ adequate values for two-dimensional solute mixing simulations following the Fick's law. The particle tracking method is an alternative to the grid-based solute transport models. The method is known as avoidable from the numerical oscillation and diffusion of the grid-based models [22]. Furthermore, the particle tracking approach is more suitable to simulate the transport of discontinuous properties of solute [23]. For these advantages, particle tracking simulations have been presented for solute transport in rivers [24] and estuaries [25,26]. The two-dimensional particle tracking model is derived from the mathematical analogy between the Fokker-Plank equation and the advection-dispersion equation [27]. Thus, the conventional particle tracking models adopted the dispersion coefficients into the random transport term following Fick's law. However, several field studies have presented skewed concentration curves, which violate the Fickian mixing properties because of complexities of channel geometries [28][29][30]. To demonstrate these non-Fickian mixing behaviors, the particle tracking method has been adopted for two-dimensional mixing simulations [24,31].
In this study, the performance of two types of particle tracking models was evaluated by comparing the simulation results. The first model is derived from the conventional particle tracking model following Fick's law. The second model is the particle dispersion model suggested by Park and Seo [31], in which the shear flow dispersion was reproduced using the sequential computations of shear translation and vertical diffusion instead of using the Fickian model. From the simulation results, the shapes of depth-averaged concentration curves were compared, and applicability of the two models was evaluated using Chatwin's parameter [32] and statistical properties.

Particle Tracking Model (PTM)
The general particle tracking model is derived from the Itō stochastic differential equation as follow [33].
where dx i = x i (t + ∆t) − x i (t) is displacement of a particle for ∆t; A i is an arbitrary vector; B ij is an arbitrary tensor; and dW j = R √ ∆t is the Wiener process, which determines the random motion of particles using a random number, R. For a number of particles and infinitesimally small ∆t, Equation (5) can be described using the Fokker-Plank equation, as shown in Equation (6) [27].
where p is the conditional probability density function. The Fokker-Plank equation mathematically coincides with the advection-dispersion equation as follows.
Therefore, A i and B ij in Equation (5) are defined using Equation (7), and Equation (5) can be rewritten as follows under the assumption, which has the coordinate system aligned with the flow direction. x From Equation (8), displacement of a pollutant particle after ∆t is determined by both the deterministic and the random transports. The fourth-order Runge-Kutta scheme was adopted to solve the ordinary differential equation shown in Equation (8). In this study, the particle tracking model using Equation (8) is denoted as PTM. The aforementioned computation procedure was depicted  Figure 2. Pollutant cloud is transported by the depth-averaged velocity fields, and the scale of longitudinal (s) and transverse (n) direction mixing is determined by the computation of the random transport term.
From Equation (8), displacement of a pollutant particle after t is determined by both the eterministic and the random transports. The fourth-order Runge-Kutta scheme was adopted to lve the ordinary differential equation shown in Equation (8). In this study, the particle tracking odel using Equation (8) is denoted as PTM. The aforementioned computation procedure was epicted in Figure 2. Pollutant cloud is transported by the depth-averaged velocity fields, and the ale of longitudinal (s) and transverse (n) direction mixing is determined by the computation of the ndom transport term. In the 2D analysis, the deterministic displacement in Equation (8) is governed by the deptheraged velocity fields calculated from the flow analysis model, Hydro Dynamic Model-2D (HDM-), which has been validated for several open channel flows [34]. The random transport term is ntrolled by the dispersion coefficient and the random number following the Gaussian distribution 5]. Generally, the dispersion coefficients are determined using empirical formulas. However, the pirical formulas have large bound of errors according to hydraulic conditions, which cause ifficulties in determining the dispersion coefficients [36]. Thus, the dispersion coefficients were mputed using the theoretical formula reported by Fischer et al. [5], as shown in Equations (9) and 0).
here    In open channel flows,  z is known as the arabolic distribution over flow depth, as shown below [5].
here  is the von Karman constant. By incorporating Equations (9) and (10) into the advectionispersion equation, previous research shows that simulation results of solute transport can display proved accuracy compared with the results using empirical formulas for dispersion coefficients 6,21]. The vertical velocity shears ( n u , s u ) in Equations (9) and (10) were reproduced using the In the 2D analysis, the deterministic displacement in Equation (8) is governed by the depth-averaged velocity fields calculated from the flow analysis model, Hydro Dynamic Model-2D (HDM-2D), which has been validated for several open channel flows [34]. The random transport term is controlled by the dispersion coefficient and the random number following the Gaussian distribution [35]. Generally, the dispersion coefficients are determined using empirical formulas. However, the empirical formulas have large bound of errors according to hydraulic conditions, which cause difficulties in determining the dispersion coefficients [36]. Thus, the dispersion coefficients were computed using the theoretical formula reported by Fischer et al. [5], as shown in Equations (9) and (10).
where u s = u s − u s and u n = u n − u n are the velocity deviations from the depth-averaged velocity in the longitudinal and transverse direction, respectively. In open channel flows, ε z is known as the parabolic distribution over flow depth, as shown below [5].
where κ is the von Karman constant. By incorporating Equations (9) and (10) into the advectiondispersion equation, previous research shows that simulation results of solute transport can display improved accuracy compared with the results using empirical formulas for dispersion coefficients [16,21]. The vertical velocity shears (u n , u s ) in Equations (9) and (10) were reproduced using the theoretical formulas and the depth-averaged flow fields from the simulation results by the HDM-2D. The theoretical formulas for u n and u s were demonstrated using the log profile [37] and the linear profile [38] as follows.

of 18
where C h is the Chezy coefficient and r c is the radius of curvature. Park et al. [21] reported that the velocity formulas shown in Equations (12) and (13) provide adequate computation results of the spatially varying dispersion coefficients in curved bends.

Particle Dispersion Model (PDM)
Fischer [39] modeled the shear dispersion in rivers into two sequential steps: the shear transport and the turbulent diffusion. Based on the study by Fischer [39], Park and Seo [31] developed the particle dispersion model (PDM), which computes the two-dimensional solute mixing using step-by-step computations of shear transport and vertical mixing. In the shear translation step, a particle position (x i ) after ∆t is determined using Equation (14).
where z 1 is the initial vertical position of a pollutant particle, ε h = αhu * is the horizontal diffusion coefficient, and α is an empirical constant according to the hydraulic conditions. According to Equation (14), deterministic displacements of particles in longitudinal and transverse directions are determined by the vertical profiles of velocity. After finishing the shear translation, the vertical position of pollutant particles belonging to the computational mesh are repositioned by the vertical mixing algorithm to finish the shear dispersion for ∆t. The vertical mixing algorithm is given below.
where n p is the number of particles in a computational mesh; a is an integer from 1 to L; L is the number of artificially divided vertical layer; β = ∆t/t m when ∆t is shorter than t m , otherwise β = 1; and t m is the completion time for vertical mixing. The value of β represents the degree of progress in complete vertical mixing. Thus, β% of entire particles within a computation mesh is evenly distributed across the artificially divided vertical layers using Equation (15), and the other particles are randomly rearranged by the turbulent diffusion shown in Equation (16). From the study of Fischer et al. [5] and Rutherford [10], t m can be defined as follows.
From the aforementioned step-by-step computations, the PDM computes the two-dimensional solute transport, as described in Figure 3. According to this figure, solute particles are transported by vertical velocity shear and the random fluctuations in the s-n plane on each vertical layer. The velocity profiles were generated using Equations (12) and (13), which are used for computing the dispersion coefficients of the PTM. The transported particles were rearranged in a vertical layer by the vertical mixing algorithm using Equations (15) and (16). Figure 3b shows the particle distributions in the s-z plane after the vertical mixing is completed.
Water 2020, 12, x FOR PEER REVIEW 6 of 18 Figure 3. Descriptions of the two-dimensional particle mixing procedures using the particle dispersion model (PDM).

Model Validations
The simulation results of the two particle tracking models were validated by comparing with the analytic solution of the advection-dispersion equation. By assuming uniform flow in a straight channel with a rectangular cross section, the concentration of instantaneously introduced pollutant can be calculated using Equation (18).
where M is the total mass of pollutant (kg) and   , which was derived by Elder [15], who assumed the log-profile for the longitudinal velocity shear. In the case of the transverse dispersion coefficient, hu  was adopted following the study by Fischer et al. [5]. In this straight channel, the flow velocity is 1 m/s and the flow depth is 0.3 m. Particle tracking simulations were conducted under the aforementioned conditions by inputting 1 kg mass of pollutant using both the PTM and the PDM. The number of particles for particle tracking simulations was 30,000, and the concentration fields were calculated by counting the number of particles in a computational mesh as follows. Figure 3. Descriptions of the two-dimensional particle mixing procedures using the particle dispersion model (PDM).

Model Validations
The simulation results of the two particle tracking models were validated by comparing with the analytic solution of the advection-dispersion equation. By assuming uniform flow in a straight channel with a rectangular cross section, the concentration of instantaneously introduced pollutant can be calculated using Equation (18).
where M is the total mass of pollutant (kg) and (s 0 , n 0 ) is the injection point. The longitudinal dispersion coefficient was inputted as D L = 5.93hu * , which was derived by Elder [15], who assumed the log-profile for the longitudinal velocity shear. In the case of the transverse dispersion coefficient, D T = 0.15hu * was adopted following the study by Fischer et al. [5]. In this straight channel, the flow velocity is 1 m/s and the flow depth is 0.3 m. Particle tracking simulations were conducted under the aforementioned conditions by inputting 1 kg mass of pollutant using both the PTM and the PDM. The number of particles for particle tracking simulations was 30,000, and the concentration fields were calculated by counting the number of particles in a computational mesh as follows.
C(s, n, t) = mn p (s, n, t) h∆A (19) where m = M/n p is the mass of a particle (kg), n p is the number of particles in a computational grid, and ∆A = ∆s∆n is a grid size (m 2 ). The concentration variations according to distance along the center of the channel are compared in Figure 4, where the range of the initial period was computed using Equation (2). The results by the PTM presented that the longitudinal dispersion was successfully reproduced by adopting Equations (9) and (12), instead of using D L = 5.93hu * . Thus, the concentration curves of the PTM were well-fitted to the analytic solution in both the initial and the Taylor periods. On the other hand, the results by the PDM show skewed concentration curves in the initial period, and the results were gradually fitted to the analytic solution after t = 30 s, as reported by Park and Seo [31]. These results show the distinct differences between the two models, where the PTM derives the results by assuming a complete balance between the shear advection and the vertical diffusion, whereas the PDM directly calculates the process of achieving the balance.
Water 2020, 12, x FOR PEER REVIEW 7 of 18 fitted to the analytic solution after t = 30 s, as reported by Park and Seo [31]. These results show the distinct differences between the two models, where the PTM derives the results by assuming a complete balance between the shear advection and the vertical diffusion, whereas the PDM directly calculates the process of achieving the balance. The shape of the temporal variations of concentration curves of the PTM and the PDM were estimated using Chatwin's parameter as follows [32,40,41].
where max t is the time when the maximum concentration ( max C ) occurs. The parameter shown in Equation (20)   The shape of the temporal variations of concentration curves of the PTM and the PDM were estimated using Chatwin's parameter as follows [32,40,41].
where t max is the time when the maximum concentration (C max ) occurs. The parameter shown in Equation (20)  After that, the deviations were gradually narrowed within the Taylor period. These results indicate that the concentration curves present tail and asymmetric curves in the initial period and they are gradually turned to symmetric curves even though there are no storage zones [31,39].

Temporal Variations of Concentration
The performance of two types of particle tracking model was evaluated by comparing with the tracer tests conducted in the meandering channel, showing the non-Fickian mixing properties [42]. The meandering channel has two bends with / 2.5 c r W  and the bends are connected with a straight part. The channel has a rectangular cross section of 1 m width and 15 m long. In this channel, the two-dimensional tracer tests were conducted using a salt solution. Here, 100,000 ppm of salt solution was instantaneously introduced as a vertical line source and the temporal changes of concentration were measured at 5 points on each of the 12 sections using conductivity sensors. Table  1 shows a summary of the experiment and the simulation conditions for flow and solute mixing analysis. The contour plot shown in Figure 6 is the velocity magnitude (U) distributions for Case C306 resulting from the HDM-2D. The velocity distributions show the increase of velocity in the inner bank of the channel bends, and the results were adequately reproduced flow in the channel bends compared with the measurements [42]. The particle tracking simulations were conducted using the results from the HDM-2D, and the vertical profiles of velocities were reproduced using Equations (12) and (13), known to adequately represent the shear flow in the meandering channel [21,31].

Temporal Variations of Concentration
The performance of two types of particle tracking model was evaluated by comparing with the tracer tests conducted in the meandering channel, showing the non-Fickian mixing properties [42]. The meandering channel has two bends with r c /W = 2.5 and the bends are connected with a straight part. The channel has a rectangular cross section of 1 m width and 15 m long. In this channel, the two-dimensional tracer tests were conducted using a salt solution. Here, 100,000 ppm of salt solution was instantaneously introduced as a vertical line source and the temporal changes of concentration were measured at 5 points on each of the 12 sections using conductivity sensors. Table 1 shows a summary of the experiment and the simulation conditions for flow and solute mixing analysis. The contour plot shown in Figure 6 is the velocity magnitude (U) distributions for Case C306 resulting from the HDM-2D. The velocity distributions show the increase of velocity in the inner bank of the channel bends, and the results were adequately reproduced flow in the channel bends compared with the measurements [42]. The particle tracking simulations were conducted using the results from the HDM-2D, and the vertical profiles of velocities were reproduced using Equations (12) and (13), known to adequately represent the shear flow in the meandering channel [21,31].   Figure 7 shows the particle tracking simulation results of Case C306 using the PTM and the PDM. The particle clouds of the PTM are more widely distributed in the longitudinal direction than those of the PDM. While the particle cloud crosses the first bend, the front of the particle cloud resulting from the PDM is transported to the inner bank of the bend, whereas the particles of the PTM are dispersed to the center of channel. It can also be seen in the second bend that the particle cloud of the PDM has a relatively shorter front compared with those of the PTM.  Figure 7 shows the particle tracking simulation results of Case C306 using the PTM and the PDM. The particle clouds of the PTM are more widely distributed in the longitudinal direction than those of the PDM. While the particle cloud crosses the first bend, the front of the particle cloud resulting from the PDM is transported to the inner bank of the bend, whereas the particles of the PTM are dispersed to the center of channel. It can also be seen in the second bend that the particle cloud of the PDM has a relatively shorter front compared with those of the PTM.
The temporal variations of concentration within the channel bends were plotted at y/W = 0.45, as shown in Figure 8. The value of concentration is normalized using the injection concentration (C 0 ). The main features of measured concentration curves are that the concentration curves have a steep, rising limb and long tail. These skewed concentration curves can be observed in the initial period [32] or the Taylor period, including a recirculation zone [43]. As the curvature of the meandering channel is not sharply curved enough to cause the recirculating flow [44], it can be said that the skewed concentration curves are not attributable to the storage effects. From Equation (2) and the hydraulic conditions shown in Table 1, the initial period was presumed to maintain for 144.3 s and 263.9 s in Case C306 and C406, respectively. Thus, the skewed concentration curves from the measurements were induced by the disproportion between the shear advection and the vertical diffusion. Even though the PTM adopts spatially varying dispersion coefficients to reproduce the mixing by the shear flow, the model shows the inherent limitations of the Fickian model. The results of the PTM presented nearly symmetric concentration distributions, of which the rising part shows large discrepancies from the measurements. On the other hand, the PDM properly reproduced the skewed concentration curves, which rapidly increase to the peak concentration and gradually decrease. These results seem to suggest that the tracer mixing is more adequately reproduced by the PDM than the PTM through direct computations of the shear dispersion process.  Figure 7 shows the particle tracking simulation results of Case C306 using the PTM and the PDM. The particle clouds of the PTM are more widely distributed in the longitudinal direction than those of the PDM. While the particle cloud crosses the first bend, the front of the particle cloud resulting from the PDM is transported to the inner bank of the bend, whereas the particles of the PTM are dispersed to the center of channel. It can also be seen in the second bend that the particle cloud of the PDM has a relatively shorter front compared with those of the PTM. The temporal variations of concentration within the channel bends were plotted at y/W = 0.45, as shown in Figure 8. The value of concentration is normalized using the injection concentration ( 0 The main features of measured concentration curves are that the concentration curves have a steep, rising limb and long tail. These skewed concentration curves can be observed in the initial period [32] or the Taylor period, including a recirculation zone [43]. As the curvature of the meandering channel is not sharply curved enough to cause the recirculating flow [44], it can be said that the skewed concentration curves are not attributable to the storage effects. From Equation (2) and the hydraulic conditions shown in Table 1, the initial period was presumed to maintain for 144.3 s and 263.9 s in Case C306 and C406, respectively. Thus, the skewed concentration curves from the measurements were induced by the disproportion between the shear advection and the vertical diffusion. Even though the PTM adopts spatially varying dispersion coefficients to reproduce the mixing by the shear flow, the model shows the inherent limitations of the Fickian model. The results of the PTM presented nearly symmetric concentration distributions, of which the rising part shows large discrepancies from the measurements. On the other hand, the PDM properly reproduced the skewed concentration curves, which rapidly increase to the peak concentration and gradually decrease. These results seem to suggest that the tracer mixing is more adequately reproduced by the PDM than the PTM through direct computations of the shear dispersion process.

Statistical Properties of Concentration Curves in Channel Bends
The differences in the simulation results by the PTM and the PDM were scrutinized using Chatwin's parameter and the statistical analysis. In Figure 9, the values of Chatwin's parameter were plotted using the concentration curves shown in Figure 8. In Case C306, both models show

Statistical Properties of Concentration Curves in Channel Bends
The differences in the simulation results by the PTM and the PDM were scrutinized using Chatwin's parameter and the statistical analysis. In Figure 9, the values of Chatwin's parameter were plotted using the concentration curves shown in Figure 8. In Case C306, both models show discrepancies from the measurements on S4, but the results of the PDM properly followed the variations in the measurements on S9. Especially, the results of the PDM were adequately fitted to the measurements of Case C406. Overall, the PTM appears to have a moderate change in the values of Chatwin's parameter compared with the measurements. These tendencies of the results from the PTM were also maintained in Case C406. The time variation of Chatwin's parameter indicates the temporal change of the concentration curves, in which the rising and falling parts of the concentration curves represent positive and negative values of Chatwin's parameter, respectively. Thus, the slope of Chatwin's parameter shown in Figure 9 is listed in Table 2 by dividing into rising and falling parts. In general, the slope of the rising part is steeper than the falling part, as can be seen in Figure 8, and the values shown in Table 2 appear to support such variations. The lower values of the slope in the falling part than the rising part present that a concentration curve has a long tail. In the results by the PDM, these results can also be found, and the values of slope were overall similar to the measurements, except on S9 in Case C306. The results by the PTM maintained an almost constant slope in both the rising and falling parts of the concentration curves. Thus, the rising parts show large errors compared with the measurements.  The time variation of Chatwin's parameter indicates the temporal change of the concentration curves, in which the rising and falling parts of the concentration curves represent positive and negative values of Chatwin's parameter, respectively. Thus, the slope of Chatwin's parameter shown in Figure 9 is listed in Table 2 by dividing into rising and falling parts. In general, the slope of the rising part is steeper than the falling part, as can be seen in Figure 8, and the values shown in Table 2 appear to support such variations. The lower values of the slope in the falling part than the rising part present that a concentration curve has a long tail. In the results by the PDM, these results can also be found, and the values of slope were overall similar to the measurements, except on S9 in Case C306. The results by the PTM maintained an almost constant slope in both the rising and falling parts of the concentration curves. Thus, the rising parts show large errors compared with the measurements. The statistical properties of the concentration curves shown in Figure 8 are compared in Table 3. The shape of the concentration curves can be defined using the peak concentration (C p ), time to the peak (t p ), time to centroid (µ t ), variance (σ 2 t ), and skewness (ξ t ). The values of C p /C 0 , t p , and µ t show similar errors in both models compared with the measurements. The obvious differences between the two models can be found in the values of σ 2 t and ξ t . The PTM tends to overestimate the variance on the first bend of the channel, whereas it provides similar results in the second bend. In the case of the skewness, the PTM maintains low skewness values in both channel bends. As can be found in Table 2, the PTM calculates the concentration curves close to the Gaussian curve, and the skewed curve in the initial period is difficult to reproduce, even though the model adopts the dispersion coefficients including the effects of vertical velocity shear. Thus, the errors in the statistical characteristics of the PTM were reduced in the second bend, where the non-Fickian mixing properties within the initial period decreased. In case of the PDM, the statistical values indicate that the model properly reproduced the shape of the concentration curves, although the skewness in the second bend was slightly overestimated.

Model Applications in the Hongcheon River
The performance of the two particle tracking models within the Taylor period was evaluated by comparing them with the 2D tracer test results conducted in the Hongcheon River by Park and Seo [31]. The Hongcheon River is a tributary of the North Han River located in the Kangwon province of South Korea. The study site shows a mildly curved bend, as shown in Figure 10, and the tracer tests were executed under the hydraulic conditions given in Table 4. For the experimental conditions, the initial period would be maintained for 22.9 min according to Equation (2), and the Taylor period is estimated to occur 64 m from the tracer injection point when the averaged flow velocity (U) is considered. Thus, the results indicate that Sections 1 and 2 in Figure 10 are included in the Taylor period. Water 2020, 12, x FOR PEER REVIEW 13 of 18 Figure 10. Outlines of the study site and the velocity distributions in the Hongcheon River. The particle tracking simulations were conducted following the conditions listed in Table 4. The simulation results were plotted in Figure 11 with the temporal changes of concentration in Section 1 and Section 2. As reported by Park and Seo [31], the tracer cloud has an asymmetric shape owing to sidewall friction and irregular channel geometries even in the Taylor period. Thus, the concentrationtime curves show a notably stretched tail of the falling part in Section 2, and these results were adequately reproduced using both models. These results revealed that the PTM is able to create the asymmetric concentration distribution that does not appear in the initial period. It can be inferred that the spatially varying dispersion coefficients computed by the vertical velocity profiles by Equations (9) and (10) are applicable to exhibit solute mixing in rivers including dead zones. The main difference between the two simulation results is that the rising part of the concentration curves by the PDM increases sharply compared with the results by the PTM.  The particle tracking simulations were conducted following the conditions listed in Table 4. The simulation results were plotted in Figure 11 with the temporal changes of concentration in Sections 1 and 2. As reported by Park and Seo [31], the tracer cloud has an asymmetric shape owing to sidewall friction and irregular channel geometries even in the Taylor period. Thus, the concentration-time curves show a notably stretched tail of the falling part in Section 2, and these results were adequately reproduced using both models. These results revealed that the PTM is able to create the asymmetric concentration distribution that does not appear in the initial period. It can be inferred that the spatially varying dispersion coefficients computed by the vertical velocity profiles by Equations (9) and (10) are applicable to exhibit solute mixing in rivers including dead zones. The main difference between the two simulation results is that the rising part of the concentration curves by the PDM increases sharply compared with the results by the PTM.
The statistical properties of the concentration-time curves were compared in Table 5. Both models show similar values in C p , t p , and µ t compared with the measurements. These results explain that both simulation results present appropriate solute mixing in the intermediate mixing region, where the influence of the initial period is negligible. The featured differences can be found in the value of ξ t , where the percentage error of the PTM was in the range of 14.9-22.6% and the PDM was in the range of 1.9-6.1%. From the results, the PDM is reaffirmed to be predominant for the non-Fickian mixing simulations compared with the PTM. Water 2020, 12, x FOR PEER REVIEW 14 of 18 Figure 11. Comparisons of temporal concentration variations in the Hongcheon River.
The statistical properties of the concentration-time curves were compared in Table 5. Both models show similar values in p C , p t , and  t compared with the measurements. These results explain that both simulation results present appropriate solute mixing in the intermediate mixing region, where the influence of the initial period is negligible. The featured differences can be found in the value of  t , where the percentage error of the PTM was in the range of 14.9-22.6% and the PDM was in the range of 1.9-6.1%. From the results, the PDM is reaffirmed to be predominant for the non-Fickian mixing simulations compared with the PTM.

Discussion
The performance of two particle tracking models was evaluated by comparing the analytic solution and the 2D tracer test results. The results consistently show that the PDM is more adaptable to reproduce non-Fickian mixing in the initial period. On the other hand, the results by the PTM were inappropriate to demonstrate the skewed concentration curves during the initial mixing process. These results revealed the inherent limitations of the Fickian model, which was derived under the assumption of equilibrium between the shear and diffusion processes [4]. Because the initially

Discussion
The performance of two particle tracking models was evaluated by comparing the analytic solution and the 2D tracer test results. The results consistently show that the PDM is more adaptable to reproduce non-Fickian mixing in the initial period. On the other hand, the results by the PTM were inappropriate to demonstrate the skewed concentration curves during the initial mixing process. These results revealed the inherent limitations of the Fickian model, which was derived under the assumption of equilibrium between the shear and diffusion processes [4]. Because the initially skewed concentration distributions originated from the velocity shear [11,39], the PDM included the effects of vertical velocity shear into the computation algorithm, as shown in Figure 3. In the case of the PTM, the model is based on the Fickian model and the shear effects were incorporated into the dispersion coefficients. Thus, the PTM shows discrepancies against the measurements in the initial period, even though the dispersion coefficients computed using the vertical velocity profiles were reported to improve the accuracy of the simulation results compared with empirical formulas [6,16,21].
The non-Fickian mixing behavior is also known to occur by breaking the equilibrium state in the Taylor period because of the complex channel geometries [28,30,31,40]. Thus, to exhibit skewed concentration distributions by the channel irregularities, previous studies adopted dead zone models [29,45] or statistical models [46] for one-dimensional mixing problems. However, these approaches have not been expanded to 2D mixing problems. The PTM computes the 2D shear dispersion using the spatially varying dispersion coefficients. From the application of the dispersion coefficients that change according to the velocity structures [21], the PTM reproduced the long tail of concentration curves induced by the channel geometries in the application results of the Hongcheon River. The results, however, were insufficient to reproduce the steep rising part of the concentration curves. An alternative approach has been reported, which employed vertical velocity formulas for the computation of the deterministic transport of each particle, instead of the Fickian approach [24,31]. The aforementioned direct computation method for shear flow dispersion was employed for the PDM, and the model successfully presented the skewed concentration curves due to the irregular channel geometries.

Conclusions
In this study, two particle tracking models incorporating the effects of vertical velocity shear were compared using the spatial and temporal variations of concentration curves. The PTM is developed based on Fick's law, in which the dispersion coefficients are employed to compute the longitudinal and transverse dispersion. The dispersion coefficients in the PTM were calculated using the triple integration of velocity deviations suggested by the study of Fischer et al. [5]. The other model denoted as PDM applied step-by-step computation algorithms to reproduce the shear flow dispersion instead of adopting Fick's law. The model also adopted the vertical profiles of velocity to compute the shear translation of particles placed on each vertical layer. The main difference of these two models is that the PTM is applicable only within the Taylor period, where Fick's law is valid, whereas the PDM can be used in both the initial and Taylor periods. Thus, to test the applicability of the two models, the models were applied to the test cases of the ideal straight channel and the meandering channel and the shape of the concentration curves was compared.
In the ideal straight channel case, the simulation results by the two models were validated by comparing with the analytic solution of the advection-dispersion equation. The concentration-distance curves of the PTM successfully reproduced the analytic solutions in the entire test domain. However, the PDM shows large discrepancies in the initial period compared with the analytic solution, and the differences were gradually reduced in the Taylor period. The analytic solution and the PTM produced symmetric concentration curves even in the initial period, whereas the results by the PDM show a skewed shape. These differences were distinctively found from the temporal variations of Chatwin's parameter, in which the PTM follows the linear line and the PDM shows deviations from the linear line in the falling part of the concentration curve. The deviations in the results of the PDM were gradually reduced to become a linear line, which indicates the symmetric concentration curve. Thus, the PDM properly reproduced the procedures of the balance between the shear advection and the vertical diffusion, which are not observed in the Fickian model.
Solute mixing simulations were conducted in the meandering channel, which shows asymmetric concentration curves due to the non-Fickian mixing features. The performance of the two models were investigated by comparing Chatwin's parameter and the statistical properties. The results by the PTM show the symmetric distribution, in which the slope of Chatwin's parameter was nearly constant and the skewness was maintained at a low value. On the other hand, the PDM presented a long tail of the concentration curves in the channel bends, and the skewness decreased in the second bend compared with the first bend, as observed from the measurements. The performance of the PDM was reaffirmed from the application results in the Hongcheon River. Even though the PTM adequately reproduced the skewed concentration curves because of the channel irregularities, the PDM more accurately demonstrated the asymmetric variations of concentration. These results reveal that the direct computation algorithm of the shear dispersion provides more accurate results than the model following Fick's law.
Author Contributions: Conceptualization, I.P.; methodology, I.P. and D.S.R.; software, I.P.; validation, H.S.; writing-original draft preparation, I.P. and J.S.; writing-review and editing, I.P. and J.S.; project administration, D.S.R.; funding acquisition, D.S.R. All authors have read and agreed to the published version of the manuscript.
Funding: This work is funded by the Ministry of Land, Infrastructure, and Transport, grant number 20DPIW-C153746-02.