Development of a Backward–Forward Stochastic Particle Tracking Model for Identification of Probable Sedimentation Sources in Open Channel Flow

: As reservoirs subject to sedimentation, the dam gradually loses its ability to store water. The identification of the sources of deposited sediments is an effective and efficient means of tack-ling sedimentation problems. A state-of-the-art Lagrangian stochastic particle tracking model with backward–forward tracking methods is applied to identify the probable source regions of deposited sediments. An influence function is introduced into the models to represent the influence of a par-ticular upstream area on the sediment deposition area. One can then verify if a specific area might be a probable source by cross-checking the values of influence functions calculated backward and forward, respectively. In these models, the probable sources of the deposited sediments are considered to be in a grid instead of at a point for derivation of the values of influence functions. The sediment concentrations in upstream regions must be known a priori to determine the influence functions. In addition, the accuracy of the different types of diffusivity at the water surface is discussed in the study. According to the results of the case study of source identification, the regions with higher sediment concentrations computed by only backward simulations do not necessarily imply a higher likelihood of sources. It is also shown that from the ensemble results when the ensemble mean of the concentration is higher, the ensemble standard deviation of the concentration is also increased.


Introduction
Sediment that is trapped behind dams reduces reservoir capacity. In some high-profile cases reservoirs have been found to be filled with sediment, not only worsening their functions and/or rendering a weakened dam infrastructure, but also posing a safety hazard [1][2][3][4]. Despite more than six decades of research, sedimentation remains probably the most serious technical problem faced by the dam industry, and sedimentation in reservoirs is now of primary global importance [5]. To study reservoir sedimentation, adequate knowledge of the sediment transport process in rivers is needed [6].
According to the Taiwan Water Resources Agency [7], most riverbeds and banks in Taiwan are prone to erosion, resulting in high sediment concentrations in most rivers. When sediment particles in a river are transported to a reservoir, the flow velocity is reduced and thus a larger volume of sediment particles is deposited at the bottom of the reservoir, which is a process referred to as reservoir sedimentation. The Shihmen reservoir, located in Taoyuan, is one of the main reservoirs in northern Taiwan. It provides domestic as well as agricultural water supply and hydroelectricity for more than three million people in northern Taiwan. In addition, the Shihmen reservoir also plays the role as a detention basin for the Taipei Basin. It has an effective storage capacity of nearly 309 million cubic meters, of which the sediment has occupied 31.6% since 2009 [7]. Sedimentation is one of the most important engineering problems that threaten reservoirs around the world. Traditionally, dredging and sluicing are passive methods, which cannot completely resolve the sedimentation problem since the sediment particles are deposited continuously during the removal process [8].
To manage reservoirs in a more sustainable way, identifying the source of the sediment particles is one step toward dealing with the problem [4].
Sediment-related issues can be prevented in advance when the probable sources of these deposited sediments are identified and prevented from flowing into the reservoirs. The turbulent characteristics of open channel flows are important in understanding particle transport phenomena, which consist of highly disordered and chaotic fluid motions over a wide range of length scales and frequencies [9][10][11][12][13][14][15][16][17]. This study examines the probable source regions of deposited sediments using a Backward-Forward Stochastic Diffusion Particle Tracking Model (Backward-Forward SD-PTM), a method of Lagrangian dynamics [18]. Not only the particle concentrations but also the particle trajectories can be simulated.
The stochastic diffusion particle tracking model (SD-PTM) is intrinsically irreversible owing to its stochastic nature. Some investigations have implied that the results of backward tracking can provide valuable clues about the probable sources of deposited sediments [19]. The probable sources are determined not only by the backward tracking process but also by the forward tracking process. Therefore, although most partial differential equations cannot be solved reversibly such as by finding the initial conditions from the results, the places of probable sources (herein considered as initial conditions) can still be identified, thus providing valuable information for solving sedimentation problems [19].
In the SD-PTM, a Brownian motion term describes the randomness of particle motion. The probability distribution of particle positions or particle concentrations can thus be computed by performing sufficient simulation runs of a single particle, or by a single simulation run with many independent particles. Notably, the distribution of the concentrations may vary among simulations. Consequently, an attempt is made herein to simulate such distribution multiple times to acquire a more accurate concentration distribution.
The objectives of this study are threefold. First, we aim to identify the probable source of deposited particles using the backward SD-PTM, and the forward SD-PTM, respectively. Next, we intend to develop a complete backward-forward particle tracking model. We then aim at enhancing the accuracy issue of estimating diffusivity coefficients when a particle has reached the water surface. Lastly, we improve the accuracy of the concentration distribution of sediment particles by computing the ensemble statistics of concentration distributions.
In Section 2, literature about the SD-PTM, backward-forward SD-PTM will be reviewed. In Section 3, some case studies applying backward SD-PTM or forward SD-PTM will be discussed. The influence function, a quantitative link between each source-receptor pair, is introduced in Section 4. Then, the model established in Section 4 will be applied in Section 5 to illustrate how to identify probable sources with the explicit and the implicit methods, respectively. In Section 6, more scenarios of the spatial sediment concentration distributions are concerned in a case study applying the backward-forward SD-PTM to identify probable sources. In Section 7, some conclusions and recommendations will be provided.

Stochastic Particle Tracking Model
Recently, in the field of sediment transport, particle-based approaches for estimating mass concentrations or concentrations of soil particles in a tidal estuary have been developed [20][21][22][23][24]. Man and Tsai [25] later presented an SD-PTM to simulate the trajectories of sediment particles in open channel flows. An SDE is a differential equation that includes at least one stochastic process term. For example, the stochastic process term can be presented as a derivative of Brownian motion or a Wiener process as follows: In Equation (1), , ( , ), ( , ), represent a stochastic process, deterministic drift term, the size of the fluctuation of the stochastic term and a standard Wiener process, respectively. According to Man et al. [25], ( , ) = ( , ) = + . Moreover, the size of the uncertain movement of sediment particles caused by turbulence is assumed to be proportional to the diffusion coefficient . ( , ) = ( , ) can then be obtained. Numerically, can be simulated as a normally distributed random variable (0, √ ) with mean zero and variance , where is the time step for numerical simulation. Therefore, the governing equation in the numerical form of the SD-PTM can be written as follows: The SD-PTM, composed of both a mean drift term and a random term, is capable of capturing any randomly selected scenarios of particle movement. Equation (2) can be coupled with a re-suspension process to simulate the movement of suspended sediment particles.
Oh and Tsai [26] added a jump term to the SD-PTM to simulate perturbation by randomly occurring extreme events including sediment deposition and re-suspension. Subsequently, Tsai et al. [27] introduced several particle tracking models to describe particle movement under various flow conditions, i.e., the stochastic diffusion process, stochastic jump process, and stochastic jump-diffusion process. Tsai et al. [24] concentrated on the jump term in the stochastic particle tracking model by studying the influence on sediment movement of the frequency and change in the magnitude of extreme events.
Oh et al. [28] suggested using the spatial distribution of particles to estimate sediment concentration and found that a high concentration is associated with a high ensemble standard deviation. Tsai et al. [29] coupled the random arrival process of particles with the SD-PTM to simulate the trajectories of sediment particles. Oh et al. [30] examined the longitudinal distance between pairs of sediment particles using the SD-PTM model. They found that a smaller particle and a higher release point leads to a longer average travel distance to the particle deposition point.

Backward-Forward Stochastic Particle Tracking Model
The usage of backward simulations is to reduce the CPU time required to run a Monte Carlo simulation for estimating particle concentrations. Spivakovskaya et al. [31] applied a forward-reverse particle tracking system that the Monte Carlo estimator for the particle concentration can also be based on. The results show that the CPU time is reduced by at least an order of magnitude compared with the classical Monte Carlo method.
Lin et al. [32] introduced the Stochastic Time-Inverted Lagrangian Transport (STILT) model of atmospheric transport, which can be used to map the upstream influence (the probable source region) of the particle deposition point (the receptor). Lin et al. [32] were concerned with preventing inconsistencies in the backward and forward simulations of the time evolution of a particle ensemble using the STILT model. That is to say, they removed apparently irreversible particle trajectories, as can be done with the help of the "influence function," which was first proposed by Uliasz and Pielke [33]. The influence function is a quantitative indicator of a source's influence on a receptor at which a tracer particle is observed: The influence function ( , | , ) denotes the influence of an upstream source at at time on a downstream receptor at at time . The delta function ( ( ) − ) represents whether the particle appears in the source.
is the total number of particles that are released at a source. Seibert and Frank [34] also applied the linear source-receptor relationship in the backward Lagrangian particle dispersion model (LPDM) for atmospheric particles. They found that the LPDM, combined with backward simulations, is an effective tool for making point measurements, which can be handled without artificial numerical diffusion. Subsequently, Batchelder [18] employed Forward-in-Time-Trajectory/Backward-in-Time-Trajectory (FITT/BITT) modeling to identify sources of ocean organisms, such as plankton. Their research showed that the BITT approach combined with the FITT approach is more efficient in locating sources of particles than the simulations using only the FITT approach. Isobe et al. [19] later developed a method that involved backward-and forward-random walk processes for identifying the sources of particles that drift on an ocean surface. They argued that although the random walk process is irreversible, it provides valuable information when the duration of the simulation is relatively short. They found some potential probable source regions by backward-in-time simulations. The authors presented a case study of East China Sea, and statistically significant sources are well specified close to the true source since 58% to 90% of source candidates are rejected experimentally.
Most current research on backward particle tracking models is in the field of atmospheric science or ocean science. However, few investigations of backward sediment transport in an open channel can be found. Hence, the behaviors of sediments in stochastic backward particle tracking processes are of interest. In this study, some backward tracking approaches that involve the SD-PTM are used to provide some hints to the identities of probable sources of sedimentation.

Case Studies of Backward or Forward SD-PTM
In this section, the flow conditions and particle parameters are adopted from the experimental data by Kaftori et al. [35], which studied the behavior of particles in the wall region of a turbulent boundary layer. The SD-PTM is applied to simulate sediment transport in open channel flow. A vast number of random scenarios of particle trajectories are modeled using SD-PTM. Then, by performing Monte Carlo simulations, the ensemble means and ensemble variances of particle positions can be analyzed.

Model Setup
The simulation is conducted on two dimensions, which are the stream-wise direction (x) and the gravitational direction (z). In the stream-wise direction, the drift term dominates the particle velocity. Therefore, the diffusivity and turbulence terms can be ignored. In the stream-wise direction, the drift flow velocity is assumed to follow the logarithm profile, consistent with Equation (4). * = 1 30 The mean drift flow velocity in the x-direction is a function of z. * is the shear velocity and is the roughness height. When the above assumptions are taken into consideration, the numerical governing Equation (2) for forward tracking in the x-direction can be written as follows.
The numerical governing Equation (2) for forward tracking in the z-direction can be written as follows. The turbulence diffusivity in the z-direction is a function of position , and is the settling velocity.
When a particle is transported to the bed, re-suspension may occur by the mechanisms that were proposed by Wu and Lin [36], who suggested that a log-normal distribution of instantaneous velocities outperforms a normal distribution of in approaching particles on the bed when considering the pickup probability threshold. According to Wu et al. [36], the threshold is shown below.

Model Parameters
In this study, the formula proposed by Rouse [37] or by Absi et al. [38] is adopted in calculating turbulence diffusivity. Rouse [37] proposed the following vertical kinematic viscosity of a fluid: , * , , represent turbulence diffusivity in the z-direction, the shear velocity, the z-coordinate of particle position, and the flow depth respectively. If Rouse's diffusivity is adopted, then the numerical governing equation in the z-direction for forward tracking can be written as follows: Absi et al. [38] proposed the following formula for the diffusivity of suspended particles.
where ( ) is particle Stokes number, is eddy viscosity, and ( ) is turbulent kinetic energy. Turbulence intensities in the fluid phase and the solid phase are denoted as ′ and ′ , respectively. The ratio of ′ to ′ is assumed to be 1 − ⁄ . The formulas for ( ), ( ) and ( ) are listed below.

Case Studies of Backward Tracking with Rouse's Diffusivity Coefficients
In this section, the diffusivity formula adopted is Rouse's diffusivity formula. Thus, the governing equations in the x-and z-directions are derived from Equations (5) and (9) respectively. The particle parameters, initial conditions, channel settings, and simulation parameters are presented in Table 1. Based on SD-PTM and the assumption that the timescale is small enough to enable the limit of the difference quotient of the governing equation to be properly estimated, the backward SD-PTM is developed by converting the operators of the numerical Lagrangian equations. This approach can be thought of as applying a negative velocity to a particle to obtain its previous location. In the explicit method, replacing the on the right-hand sides of SDEs, Equations (5) and (9), with and vice versa on the left-hand side, yields the following equations:

Implicit Method
The numerical Equations (5) and (9) can be rewritten as follows: During the backward tracking process, each is derived from a known . Accordingly, can be regarded as an unknown variable while is known. Thus, the numerical expression is implicit because the unknown cannot be directly calculated.
Then can be represented as follows.
Notably, the backward tracking process remains a stochastic process because the Ν(0, √dt) term is included in parameters and . By the Taylor expansion of Equation (19), is given by the following formula: (20) is defined in Equation (21) as follow: Therefore, implicit numerical governing equations are derived as Equations (16) and (20).

Case Studies with Backward Tracking Using Diffusivity Coefficients of Absi et al. [38]
The previous version of the model with diffusivity coefficients that were estimated using Rouse's formula had a fundamental problem: it provided limited information about the particle position in the z-direction when a particle had reached the water surface. When = ℎ, the diffusion coefficient and the turbulence term would be zero (Equation (9)). Hence, particles remained at the surface once they had reached it in a natural system. Therefore, the accuracy of diffusivity coefficients is improved in this section by adopting the diffusivity proposed by Absi et al. [38]. In both explicit and implicit methods, the numerical governing equations are very similar to those in Section 3.3. The non-zero value of the diffusivity of Absi et al. [38] provides a means of tracking particles in the z-direction after they have reached the surface. Furthermore, since the formula of Absi et al. [38] considers the lag between fluid particles and sediment particles, diffusivity coefficients in this version may be more appropriate than the previous version using Rouse's formula.

Comparison Among Backward and Forward Tracking Cases
Additional to backward tracking, forward tracking is also exerted in Section 3. Figures  1-7 present the results of simulations that involved scenarios, ensemble means of trajectories, probable source locations, and the distributions of probable source locations in both stream-wise and vertical directions. In those figures, a and b show backward tracking tests with the Rouse diffusivity, performed using the explicit and implicit method respectively; c and d represent backward tracking tests with the diffusivity of Absi et al. [38], conducted using the explicit and implicit method respectively; e and f show forward tracking tests carried out using the explicit method with the diffusivity of Rouse [37] and Absi et al. [ 38] respectively. In Figures 1 and 2, the red line represents the water surface at z = 0.0284 (m). Figure 1. Sample paths randomly selected from 10,000 simulations: (a) sample paths of backward tracking tests with the diffusivity of Rouse [37], using the explicit method; (b) sample paths of backward tracking tests with the diffusivity of Rouse [37], using the implicit method; (c) sample paths of backward tracking tests with the diffusivity of Absi et al. [38], using the explicit method; (d) sample paths of backward tracking tests with the diffusivity of Absi et al. [38], using the implicit method; (e) sample paths of forward tracking tests with the diffusivity of Rouse [37], using the explicit method; (f) sample paths of forward tracking tests with the diffusivity of Absi et al. [38], using the explicit method.  [37], using the explicit method; (b) mean trajectory of backward tracking tests with the diffusivity of Rouse [37], using the implicit method; (c) mean trajectory of backward tracking tests with the diffusivity of Absi et al. [38], using the explicit method; (d) mean trajectory of backward tracking tests with the diffusivity of Absi et al. [38], using the implicit method; (e) mean trajectory of forward tracking tests with the diffusivity of Rouse [37], using the explicit method; (f) mean trajectory of forward tracking tests with the diffusivity of Absi et al. [38], using the explicit method.  [37], using the explicit method; (b) positions of probable sources in backward tracking tests with the diffusivity of Rouse [37], using the implicit method; (c) positions of probable sources in backward tracking tests with the diffusivity of Absi et al. [38], using the explicit method; (d) positions of probable sources in backward tracking tests with the diffusivity of Absi et al. [38], using the implicit method; (e) final particle positions in forward tracking tests with the diffusivity of Rouse [37], using the explicit method; (f) final particle positions in forward tracking tests with the diffusivity of Absi et al. [38], using the explicit method.  [37], using the explicit method; (b) final probable source location distribution in backward tracking tests with the diffusivity of Rouse [37], using the implicit method; (c) final probable source location distribution in backward tracking tests with the diffusivity of Absi et al. [38], using the explicit method; (d) final probable source location distribution in backward tracking tests with the diffusivity of Absi et al. [38], using the implicit method; (e) final particle location distribution in forward tracking tests with the diffusivity of Rouse [37], using the explicit method; (f) final particle location distribution in forward tracking tests with the diffusivity of Absi et al. [38], using the explicit method.  [37], using the explicit method; (b) final probable source location distribution in backward tracking tests with the diffusivity of Rouse [37], using the implicit method; (c) final probable source location distribution in backward tracking tests with the diffusivity of Absi et al. [38], using the explicit method; (d) final probable source location distribution in backward tracking tests with the diffusivity of Absi et al. [38], using the implicit method; (e) final particle location distribution in forward tracking tests with the diffusivity of Rouse [37], using the explicit method; (f) final particle location distribution in forward tracking tests with the diffusivity of Absi et al. [38], using the explicit method.
(e) (f) Figure 6. Final particle location variance in the x-direction: (a) final probable source location variance in backward tracking tests with the diffusivity of Rouse [37], using the explicit method; (b) final probable source location variance in backward tracking tests with the diffusivity of Rouse [37], using the implicit method; (c) final probable source location variance in backward tracking tests with the diffusivity of Absi et al. [38], using the explicit method; (d) final probable source location variance in backward tracking tests with the diffusivity of Absi et al. [38], using the implicit method; (e) final particle location variance in forward tracking tests with the diffusivity of Rouse [37], using the explicit method; (f) final particle location variance in forward tracking tests with the diffusivity of Absi et al. [38], using the explicit method.  [37], using the explicit method; (b) final probable source location variance in backward tracking tests with the diffusivity of Rouse [37], using the implicit method; (c) final probable source location variance in backward tracking tests with the diffusivity of Absi et al. [38], using the explicit method; (d) final probable source location variance in backward tracking tests with the diffusivity of Absi et al. [38], using the implicit method; (e) final particle location variance in forward tracking tests with the diffusivity of Rouse [37], using the explicit method; (f) final particle location variance in forward tracking tests with the diffusivity of Absi et al. [38], using the explicit method.
According to Visser [39], though, the diffusivity should be evaluated at a slightly different location from for self-consistency of the process. In this section, we focus on the comparisons between diffusivities calculated by the functions of Rouse [37] and Absi et al. [38]. Therefore, we adopt the suggestion of Ross and Sharples [40], which is to output the particle positions of each time step to observe the difference between the aforementioned two numerical methods. It is shown in Figure 1a,b that particles would remain at the water surface once they reach the water surface due to the zero-vertical diffusivity of Rouse [37] at the water surface. However, this limitation can be eliminated with the diffusivity of Absi et al. [38]. From Figure 1b,c, the scenarios with the diffusivity of Absi et al. [38] showed that particles would not have to remain at the water surface, which is closer to reality. The line-like pattern at the tops of Figure 3a,b also demonstrates the above observation, as the line-like pattern consists of a great number of particles at the water surface. Furthermore, the lag between fluid particles and sediment particles is considered in the diffusivity of Absi et al. [38]. Therefore, the simulation results with the diffusivity of Absi et al. [38] would be made more meaningful than the diffusivity of Rouse [37].
In the aspect of ensemble mean of trajectories (Figure 2), no significant difference between these two diffusivities except for the case of backward tracking with diffusivity of Absi et al. [38]. Figure 2c,d shows that the final position of the mean trajectory by the explicit method is higher than that by the implicit method. It may result from the limitation of implicit methods in which the unknown cannot be calculated directly. Coupling with the lag between fluid particles and sediment particles considered in the diffusivity, the discrepancies between explicit and implicit methods of diffusivity of Absi et al. [38] consequently might be obvious.
By comparing Figures 1 and 2, one can realize the importance of utilizing a stochastic particle tracking model to assess the variability of particle movement and, subsequently, the likelihood of particle sources.
With reference to Figure 3a,b,c,d, in all cases of backward tracking with the same initial conditions and parameters, the results reveal that the source is probably at x = −1.5 and z = h. To verify whether particles are deposited on the riverbed at x = 0, two forward tracking tests are carried out to determine the trajectories of particles that are released from the aforementioned point (x = −1.5, z = h). In forward simulations, the assumptions and parameters are the same, but the initial particle positions and tracking directions are different from backward simulations. For forward simulations, the particles are released from the water surface at x = −1.5 m. Two forward tracking cases are modeled by the explicit numerical method (Equations (5) and (6)) with the diffusivity of Rouse [37] or Absi et al [38].
With reference to figures that show particle positions 10 s earlier than the time that sediment particles deposit on the riverbed at x = 0 in the backward tests (Figure 3a,b,d), most of the particles accumulated near the riverbed. Furthermore, with regard to the x-direction, the particles in the near-bed region are closer to the deposition location while those that remain at the surface are farther from that location. This fact may be attributed to the use of the logarithmic velocity distribution, which is larger near the surface, and has the potential to drift particles on the water surface farther. Since the simulation is not run for a long enough time to allow all of the particles to reach the water surface, in more than half of the 10,000 simulations, particles did not arrive at the surface. Therefore, most particles are still close to the riverbed and remain close to the deposited location. Hence, the distribution of final particle locations in the x-direction is skewed toward the origin (Figure 4a,b,d).
In the z-direction, the probabilities associated with final particle locations decrease sharply as z increases. Notably, the probability associated with the water surface is the second highest owing to the intrinsic characteristic of the governing equation in the zdirection when Rouse's diffusivity is adopted (Figure 5a,b). However, in tests in which the diffusivity of Absi et al. [38]. is used, the probability that a particle is at the water surface is high and the probability that a particle is in its neighborhood is non-zero ( Figure  5c,d). On account of the logarithmically distributed mean drift flow velocity in the horizontal direction, a higher z may be correlated with a greater distance traveled by a particle. As a result, the distribution in the x-direction has two humps because the probability that a particle is near the riverbed or the surface is higher than the probability that it is in any other region. Furthermore, particles are distributed more uniformly in the explicit backward case in which the diffusivity of Absi et al. [38] is used.
The variances of particle positions in both directions increase with time. However, the variance in the x-direction increases continuously with time while that in the z-direction tends to an asymptote. In x-direction, uncertainty is caused by the re-suspension mechanism and sediment concentration gradients in the streamwise direction. Such concentration gradients may cause the shear in the flow to mix particles out which cause the variance to grow in time. If a particle is not re-suspended, then it would stagnate at its location for some time. Therefore, the variability in the stream-wise direction continues to increase. In the vertical direction, the variance initially increases because of the turbulence term. Additionally, it tends to reach an asymptotic constant (or may even decrease if the simulation time is very long) as more particles reach the water surface and their vertical positions cease to change.
Two forward tests yield insignificantly different results. The overall trajectory is downward in the downstream direction of the channel. According to Figure 5e,f, approximately one-fifth to one-fourth of particles are near the bed. The total number of particles decreases as z increases. Since the overall trend is downward, few particles stay close to the water surface. Two peak values are observed in the distributions in the x-direction at the final time. The left peak is attributable to a large number of particles near the bed, which are transported over short distances. The right peak is attributable to particles that might have arrived from a higher region where flow is fast enough to move particles farther.

Development of Backward-Forward Stochastic Particle Tracking Model
The comparison between backward cases and forward cases demonstrates that some probable source regions may be correctly identified. A large proportion of particles are deposited at the deposition point (0, 275 × 10 ), indicating that the location (−1.5, 0.0284) might be a probable source. However, no quantitative link or relationship between each source-receptor pair exists. The possibility associated with each probable source that is estimated by backward tracking cannot be quantified even when the forward tracking process is carried out for each of them. Therefore, inspired by Lin et al. [32] , the influence function is now introduced to establish or reject each candidate "probable" source that is estimated by backward simulation in this section.

Influence Function
Uliasz et al. [33] introduced the influence function. In our study, ( , | , ) is defined as the influence function of an upstream probable source that is located at at time t, with respect to a downstream receptor that is located at at time .
In Equation (22), is the number of total particles which exist at the source at time . Thus, the spatial distribution of concentration at time should be known, so at each source could be calculated. The delta function δ( ( ) − ) represents whether a particle from the upstream source ( , ) is transported to the downstream receptor ( , ). If the particle satisfies the above conditions, then the delta function δ( ( ) − ) equals one; if it does not, the value of the delta function is zero. Equation (23) can be used to calculate the influence functions of the forward particle tracking process ( , | , ). The denominator is while the numerator is the number of particles that are transported from a certain source to the receptor during the forward particle tracking process . In the backward particle tracking process, the concept of the influence function ( , | , ) is physically the same as that in the forward particle tracking process. The total particle number at the source at time significantly affects the value of the influence function. If the process is reversible, then ( , | , ) should resemble ( , | , ). However, in numerical simulations, ( , | , ) is rarely equivalent to ( , | , ). A threshold is introduced here.
For a large number of particles, if the difference between ( , | , ) and ( , | , ) is smaller than ε, then they are regarded as equal, and the process is regarded as reversible (Equation (25)).
Notably, only the numbers of particles at ( , ) and ( , ) instead of the transport trajectories are considered in the influence function. Namely, although ( , | , ) is close to ( , | , ), it is ignored by the influence function whether all the backward trajectories have their corresponding forward trajectories.
The definition of influence function that is used during the simulation process in this study differs slightly from that used in the study of Lin et al. [32]. In the study of Lin et al. [32], the same number of particles (10,000) was released from both the source and the receptor, and the air density gradient was taken into account in modifying the influence function.
corresponds to the number of scenarios in the Monte Carlo simulation in their study. However, in this study represents not only the number of scenarios in the Monte Carlo simulations but also the total number of local particles at the sources. The numbers of local particles all exceed 50 and are statistically significant in our Monte Carlo simulations.     Figure 9b, is the total number of solid and hollow blue particles in the blue rectangle.

Simulation Process
represents the total number of particles at the source at time t. The hollow particles are traced from the receptor in the backward simulation process. 5. Use the results in step 3 with the forward stochastic particle tracking model. Release particles from the middle point of the upstream grid cell to find the number of particles that travel to the receptor , as the total number of hollow red particles in the red rectangle in Figure 9c. 6. Combine the results in step 5 ( ) with the total number of particles at probable sources ( ) to calculate forward influence functions ( ). (Figure 9d) 7. Calculate the as well as repeat step 4 to step 6 for each probable source grid identified by the pure backward simulation. 8. Map the values of the forward and backward influence functions in the domain at their corresponding source locations. 9. Compare and discuss the forward and backward influence functions to evaluate the reversibility and the possibility of probable sources in each grid cell. Determine the threshold value and find whether the difference between and is smaller than .
10. Map the most probable source regions.

A Test of Backward and Forward Influence Functions
In this section, a test of backward and forward influence functions is demonstrated with the same domain as that in Section 4.2. The parameters of the particles, flow conditions and governing equation are same as those in the previous section shown in Table 1. The number of deposited particles at the receptor ( ) is assumed to be 1,000. The receptor is at the origin as the initial condition of backward tracking. For each forward tracking process, the initial condition of particle position is at the middle point of the grid of probable sources. The simulation time for both backward and forward tracking is 10 s, which is long enough to observe the trend of the motion of sediments.
is relevant to sediment concentration. Thus, the concentration of sediment particles can be introduced below.
The numerical methods applied in this test are implicit backward tracking and explicit forward tracking, all with the diffusivity proposed by Absi et al. [38]. The particle concentrations at time (when particles are transported to and deposited at the receptor) are assumed to follow Rouse's sediment concentration profile [37]: where ( ) is the concentration at height ; ( ) is the reference concentration at height = ; ℎ is the flow depth; is the settling velocity; is the Karman constant, and * is the shear velocity.
In this test, the reference concentration ( ) = 1,000 is assumed at the bottom layer to represent the number of particles there. The sediment concentration is assumed to be uniform in the stream-wise direction. Therefore, sediment concentration is not a function of x but a function of z. The sediment concentrations within the domain at time is estimated by backward stochastic particle tracking. Then, the sediment concentrations (total numbers of sediment particles in grid cells, , within the domain) at the source are known. In Figure 10, the black dashed line denotes the numbers of particles that are calculated using Rouse's sediment concentration profile (at time ), while the blue line denotes the number of particles that are estimated by backward stochastic particle tracking (at time ). Both concentrations are independent of x but dependent on z. To confirm the accuracy of backward tracking, and are compared. Every is regarded as a true value. Then, the relative error of every is calculated and mapped. The relative error is given by the following equation.
= (27) The difference between and in each grid cell is mapped out. Figure 11 presents the influence functions that are calculated by backward tracking ( ) and forward tracking ( ). The blue regions are those in which no particle appears ten seconds before time = .
The influence function of an upstream grid cell on the receptor (0, 275 × 10 ) is shown in a lighter color. The yellow dots in Figure 11c denote the locations of 1,000 particles 10 s earlier before they deposited at the receptor. The results are calculated by backward tracking. At = −0.9, the grid cells closer to the riverbed tend to have large and small values. The highest sediment concentrations can be obtained only when locations that are determined by backward simulations are taken into account. However, these grid cells may not be the regions of the most probable sources because of their irreversibility, which is implied by the great difference between the backward and forward influence functions. To identify the probable source regions, a threshold ε is determined to decide whether the difference between a pair of backward and forward influence functions for a grid cell could be ignored. Figure 12 shows the identified probable source regions with threshold values of 1, 0.8, 0.6, 0.4, 0.2, and 0.1. It can be found that with a lower threshold value, fewer grid cells are probable source regions. Again, the results that are estimated only by backward simulations may not be precise. Nevertheless, the process of backward tracking can assist in scaling down the probable source regions.

A Case Study with Explicit and Implicit Backward Tracking Methods
The objective of this case study is to identify the probable source regions of these particles 10 s before the time (time = − 10) using the backward-forward stochastic particle tracking model. At time , sediment particles deposit on the riverbed at x = 0.

Parameters and Initial Conditions
The diffusivity applied in this case study is that of Absi et al [38]. Two numerical methods-explicit and implicit-are used during backward simulation while the numerical methods used for forward simulations are always explicit.
In this case study, the domain is the same as the one in Section 4.2. The particle parameters, channel settings, and flow conditions are mostly the same as those in Table 1. The only difference is that the number of deposited particles is 10,000. The concentration is assumed to follow Rouse's sediment concentration profile as a function of depth with 10,000 particles in the bottom layer at time . The sediment concentration is assumed to be uniform in the x-direction. The sediment concentration 10 s before (at time ) is estimated using implicit backward simulations. Figure 13 plots the sediment concentrations as functions of z at times and as blue and black curves, respectively. With the help of this concentration data, as well as the influence function for each probable source candidate can be calculated.
in each probable grid is also the number of realizations in Monte Carlo simulations.

Simulation Results and Discussions
Figures 14-16 present the simulation results obtained using explicit and implicit methods. The results obtained using the explicit method are plotted on the left-hand side while those obtained using the implicit method are plotted on the right-hand side. Figure 14 plots as white dots the 10,000 particle source locations that are estimated from 10,000 backward simulations from the receptor. As discussed in Section 3, the particles tend to cluster near the riverbed and near the water surface.  Figure 15 displays backward influence functions, forward influence functions and the relative error of influence functions. In general, regions with higher backward influence functions seem to appear along with a band from the upper-left to the lower-right. The grid cells with higher backward influence function regions determined using the explicit method significantly differ from those determined using the implicit method. Two simulations that involve the same processes yield slightly different results because of the stochastic term in the stochastic particle tracking model. This slight difference also suggests that the influence functions and concentrations may not be sufficiently accurate. This is because the concentration distribution obtained based on 10,000 backward simulations only represent a single scenario of the concentration distribution. In regions adjacent to the upstream/downstream area, the forward influence functions may be low, indicating that the influence of this area on the receptor is weak. Accordingly, the backward simulation process roughly identifies the most influential regions, but influence functions are needed to determine probable source regions.
Overall, the relative error in the implicit method is much higher than that in the explicit method, especially in the lower left regions. Therefore, an implicit method may be less accurate in most areas in this case. Notably, the spatial distributions of regions with lower relative errors, corresponding to greater accuracy, differ significantly between the explicit and implicit methods. The left-most section has the least relative error in the backward influence functions. Both the implicit and the explicit methods have a larger error near the riverbed than at the water surface. Furthermore, the regions with higher sediment concentrations computed by only backward simulations do not necessarily imply a higher likelihood of sources. While the explicit method may not identify correctly the probable source regions, the implicit method could assist in identifying such missed probable sources.
After the relative error of the influence functions is calculated, the threshold is set to different values, as depicted in Figure 16. If the threshold is 0.1, the relative error of influence functions is smaller than 10%. The most probable source regions are then identified with a threshold value of 0.1, and their corresponding influence functions are listed in Table 2.  In general, the sources that are identified by the explicit method tend to have a higher influence function, so the explicit method identifies locations that have a stronger higher influence on the receptor. However, the implicit method performs poorly in identifying locations with higher influence functions.

Discussion About Probable Sources Based on Ensemble Mean and Ensemble Standard Deviation of Spatial Sediment Concentration Distribution
Oh et al. [28] suggested that sediment concentration can be estimated by the spatial distribution of particles using the stochastic particle tracking model. As displayed in Section 3, one numerical experiment that involved 10,000 simulations can provide the ensemble mean and the ensemble variance of particle locations. However, as depicted in Figure 14, the same initial conditions can lead to different spatial concentration distributions.
Hence, the numerical experiment that involved 10,000-round Monte Carlo simulations to determine one unique spatial concentration distribution is repeated 100 times.
Then, 100 scenarios of spatial concentration distributions can be obtained. In this section, the ensemble mean and the ensemble standard deviation of spatial concentration distributions in backward tracking are computed. The position of the receptor is also located at the origin. Since the probable regions that are identified in Section 5 are of great interest, 100 spatial concentration distributions are also computed by forward tracking from these probable source regions. The ensemble statistics of the spatial concentration distributions and the influence of these sources to the receptor can thus be discussed. Table 3 presents the parameters of particles, flow conditions, simulation time and the time step. In the backward simulation, 10,000 particles are traced backward for 10 s from the receptor (0, 275 × 10 ) for each realization of spatial sediment concentration distributions. With respect to forward simulations, 10,000 particles are traced forward for 10 s from each probable source location identified in Section 5 for each realization of the spatial sediment concentration distributions. A total of 100 realizations of spatial concentration distributions are computed to estimate the ensemble means and ensemble standard deviations. As Figure 16 indicates, when the threshold value is 0.1, eight grid cells and five grid cells are suggested to be the probable source regions in explicit and implicit backward methods, respectively.  Figure 17 presents the ensemble mean and ensemble standard deviation of spatial sediment concentration distributions that are estimated using explicit and implicit methods. Most of the particles are traced back to the riverbed at about = −1. With respect to the ensemble means, the highest concentration obtained by the explicit method is farther from the receptor than that obtained by the implicit method. Another region of high concentration appears on the water surface at about = −1.5. The particles tend to be distributed more uniformly in the explicit method. The highest ensemble means of concentration that are obtained by the explicit method and implicit method are 975.93 and 4017.1 particles per grid cell, respectively. The difference is considerable, implying that at least one of these two values are far from the true value. The maximum of the ensemble standard deviations obtained by the explicit and implicit methods, are 30.36 and 26.60, respectively. The difference of the ensemble standard deviations in the grid cell with the highest concentration between the two numerical methods is quite low. However, compared to the ensemble mean value, the ensemble standard deviation in the grid cell with the highest concentration is relatively high in the explicit method.

Parameters and Initial Conditions
As depicted in Figure 17, ensemble standard deviations tend to be higher where ensemble means are higher, according to both the explicit method and the implicit method. This result is consistent with those obtained by Oh et al. [28]. This highest standard deviation may be a significant cause of why only the grid cell with the highest concentration, estimated only by the backward simulation, may not necessarily be the most probable source. However, another reason why the area of highest concentration is not the region of the most probable source may be attributed to the large difference between the estimated ensemble mean of the sediment concentration and the true value. Table 4 displays the ensemble mean and ensemble standard deviation of the spatial sediment concentration distribution from eight probable sources as determined using the explicit method in Section 5. The blue color indicates an absence of particles. A lighter color represents a higher ensemble mean or ensemble standard deviation. The red grid cell rep-resents the location of a probable source from which particles are released. For these aforementioned simulations the spatial sediment distributions cover the receptor grid cell in all simulations, demonstrating the superior performance of the explicit method.  Table 5 presents the ensemble mean and ensemble standard deviation of concentrations (represented as the number of particles) in the receptor grid cell. In all simulations, the concentration of the particles at the receptor exceeds 500 particles per grid. Namely, the probability that a particle is transported to the receptor from these sources exceeds 5%. This probability is relatively high, as the sediment is present over a large area. Moreover, the grid cells with the first or second highest concentration are receptor grid cells in most of the simulations. However, a grid cell with a higher ensemble mean concentration normally also has a higher ensemble standard deviation (variability) in forward simulations.  Table 6 presents the ensemble mean and ensemble standard deviation of the spatial sediment concentrations when particles are traced forward from the probable sources that are identified by the implicit method. The notation and representation of the results are similar to those in the explicit method. As the table shows, the receptor was covered by the sediment in all five simulations. However, in most of these simulations, the ensemble mean of the concentrations at the receptor was lower than that obtained using the explicit method. For simulations Im-2, Im-3 and Im-4, the concentration is around 200 to 300 particles per grid cell (about 2-3%); for simulation Im-1, the ensemble mean concentration is only 29.3 particles per grid cell (0.293%) ( Table 7). Only the ensemble mean of the concentration in simulation Im-5 is as high as that obtained using the explicit method. Therefore, the implicit method underperforms the explicit method in determining sources with a strong influence on the receptor. Again, a higher ensemble mean of concentration corresponds to a higher ensemble standard deviation in these simulations. None of the simulation cases is an exception to this rule. A grid cell at (−0.9, 0) is a high concentration, as estimated by either the backward explicit method or implicit method. In Section 5, neither the explicit nor the implicit method suggests this region as a probable source. Now, the ensemble mean and the ensemble standard deviation of the concentrations are computed after the particles have traveled for 10 s (Figure 18). Although the receptor grid is covered by the particles, the concentration of sediment is low in the receptor grid cell. Only 31.77 of 10,000 particles are transported to the receptor from this region. The results again demonstrate that the highest concentration, estimated by the backward simulation alone, suffices to identify the probable source. Furthermore, the probability of particle transporting to the receptor from a non-probable source is much lower than that from probable sources, no matter which numerical methods is adopted.

Conclusions and Recommendations
Source identification is useful for containment prevention or ecological conservation. It also supports our understanding of particle behaviors, which can aid in efficiently tackling sedimentation problems. However, there is limited research applying the backwardforward particle tracking model to identify sedimentation sources. In this work, a stateof-the-art Lagrangian stochastic particle tracking model is applied to identify sedimentation sources. The diffusivities proposed by Rouse [37] and Absi et al. [38] are both considered. It is found that not only the non-zero value of the diffusivity is performed at the water surface, but the lag between fluid and sediment particles can be considered by the diffusivity of Absi et al. [38].
Inspired by Lin et al. [32], a backward-forward stochastic particle tracking model is formed with the help of influence functions. The sediment concentrations in upstream regions must be known a priori to determine the influence functions. A case study of source identification using the stochastic particle tracking model is presented. In this backward-forward SD-PTM, the outcome from the backward simulation can be used as a reference in limiting the regions of probable sources. With respect to the numerical method, the explicit method performs better in identifying sources with higher influence functions than the implicit method does. In the forward simulations, the implicit method may underperform the explicit method when determining probable sources. Moreover, the probabilities are much higher comparing to simulations from non-probable sources (only about 0.3%) no matter which numerical method is applied.
The area of high concentrations also has a high standard deviation in concentrations. The results of probable source regions only identified by backward simulations may not be indicative, but they scale down the probable source regions. Both the forward and the backward influence function can be improved when calculated using the ensemble mean of particle concentrations. However, the computation time would be very large if all of the concentrations are ensemble means. In addition, the trajectories are not presented as outcomes of the backward-forward SD-PTM and inputs of influence functions. Whether there exist corresponding trajectories in pure backward and forward simulations is considered in deciding influence functions, then the backward-forward SD-PTM could be made more convincing. For example, particle trajectories could be simplified, and the similarity between patterns in the forward and backward simulations could be considered.
Finally, the computational time for backward-forward particle tracking in this preliminary study remains considerable. To solve a large spatial-temporal scale real-world engineering problem such as reservoir sedimentation, the computational demand may be met by more efficient numerical approaches and advanced computing technology.. Data Availability Statement: All data used in this study are cited in the tables in this paper.