Abstract
In this work a fast computational particles tracer model is developed based on Particle-In-Cell method to estimate the sediment transport in the access zone of a river port area. To apply the particles tracer method, first it is necessary to calculate the hydrodynamic fields of the study zone to determine the velocity fields in the three directions. The particle transport is governed mainly by the velocity fields and the turbulent dispersion. The mechanisms of dispersion and resuspension of particles are based in stochastic models, which describes the movement through a probability function. The developed code was validated using two well known cases with a discrete transformation obtaining a max relative error around 4.8% in both cases. The simulations were carried out with 350,000 particles allowing us to determine under certain circumstances different hydrodynamic scenarios where the zones are susceptible to present erosion and siltation at the entrance of the port.
1. Introduction
The analysis of sediment transport is extremely important for predicting siltation and contaminant transport in rivers or sea ports. Sediment transport is a crucial process to understand the environmental state of aquatic systems [1], because contaminants in aquatic systems are usually cohesive or hydrophobic. Therefore, it is necessary to estimate accurately the dynamics of the sediment transport to obtain a reliable estimation of the contaminants in the water.
Numerical simulations have shown to be efficient, economic, and practical tools for predicting sediment transport in complex hydrodynamic systems. Nevertheless, to obtain a reliable numerical simulation results, it is necessary to take into account as many parameters as possible, to predict correctly the erosion and deposition. For this work, the number of particles and their size contained in every cell are considered to determine if the concentration, erosion, and deposition are governed by probabilistic functions.
Erosion is the flux of particles from a sediment bed into the overlying water, and deposition is the flux of particles back to the sediment bed. The movement of the water composed by the advection, the turbulent dispersion, and settling are the principal phenomena that affect the sediment transport. For this reason, it is necessary to calculate the velocities field of the water to obtain accurate results.
The main objective of this work is to validate a three-dimensional computational model for the sediment transport, developed specifically to obtain results in a short time period. The model is supported in the Particle-In-Cell method and was implemented in MATLAB; it is very fast, consumes low computer resources, and yields good results. The selected study case is a small port terminal in the Magdalena river in Colombia.
As all the numerical models for the sediment transport, it is necessary to use a hydrodynamic model to estimate the velocities field. The hydrodynamic model used has been proven to be robust and precise, and many research has been carried out with this model and they can be consulted in [2–4].
The numerical computational model developed in this work for the transport of the particles is governed by the Lagrangian approach, where the particles are located following a concentration exponential law or randomly located in the space. The advantage of using Lagrangian models to estimate sediment transport and some temporal changes in the morphology of the bottom lies in the computational speed of using a previously calculated hydrodynamic field for the movement of particles. These allow the approximation of the temporal concentration of sediments contemplating the density of the material and using the PIC method to quantify the sediment transport associated with the displacement of the particles near the bottom.
The Lagrangian approach is widely used in the study of the trajectories of movement of solid particles in fluid environments. This is due to the fact that it is possible to track the movement for each specific particle in more detail, in comparison with determining average concentration for grid cells. However, in practice, the most difficult part to use a Lagrangian based method is the strong dependence on the performance of computational resources, such as the amount of memory required for the particles.
One of the pioneering jobs in the application of the Lagrangian technique regarding particle tracking is introduced by Hunter [5], where the method is described and its superiority with finite-difference and finite-element methods is compared for the prediction of dispersion of contaminant in fluid. The method approximates with greater exactitude the nonlinearity of the advective term in the transport equation. For this reason, it is commonly applied in coastal oceanography [6]. In fact, the Lagrangian method is widely used in the dispersion of cluster of particles, having its origins in atmospheric modelling or thermal engineering, and today is used as a well proven method, which can be reviewed in very recent applications [7, 8].
For particle tracking modelling specifically for the sediment transport research, it is necessary to take into account erosion, transport, settling, deposition, and resuspension. In general, the submodule of tracking for sediments is decoupled of the hydrodynamic model, an overview how to perform particle tracking modelling and setup could be consulted in [9]. To carry out the sediment transport study, it is necessary to first calculate the hydrodynamics and subsequently couple the module of sediments, such as, for example, the module developed for sediments for the Environmental Fluid Dynamics Code provided by EPA [9], and it is common for the modelling of sediment to be done this way, such as DEDPLUME-RW or SEDTRAIL-RW [10].
Basically, to the modelling of particle tracking hydrodynamics, the field used corresponds to a stable field and with convergence in time; then the velocity and the turbulent parameters can be considered constant or periodic. The temporal variation could be bigger than the used for the hydrodynamic module. Thus, the particles transport simulation could be carried out for bigger simulation times than the obtained in the hydrodynamic. Therefore, the simulation of the particle transport could use the hydrodynamic field repeatedly until the desired simulation time period is reached. The velocities that move the particles are located in the sides of a three-dimensional cell that build a mesh. It is possible to approximate the velocities in any point inside the cell using a linear interpolation. In fact, this modelling process is applied for any pollutant discharge governed by the transport equation [11] as thermal plume discharges [4, 12].
The Lagrangian mathematical approximations, based on methods of random movement, are well established tools for the calculus of pollutant discharges into aquatic environments. The discharges are treated as a finite number of particle; these particles move under the influence of the previously established flow field. The amount of contaminant represented by each particle in the model can be decremented through the time to represent processes of disappearance. With this method, it is relatively easy to estimate the concentration that is required at a given volume. It is necessary to just carry out a simple calculation to relate the number of particles contained in a fixed space with the number of particles included in a mesh of known concentration [13].
Basically in a Lagrangian approach consisting in determining the next accurate location for each particle in the computational domain, a particle located at the position is moved to the position at time , where , , and , are expressed as follows:where the drift terms are governed for river flows by the advection-diffusion and the random displacements by some probabilistic functions. For this reason, in this research the advection of a diffusion is tested under controlled cases to determine the precision with different numbers of particles and a probabilistic function is proposed.
Additionally, among the benefits of the Lagrangian models over Eulerian Advection-Diffusion models is the computational efficiency [14]. However, it requires the use of mobile coordinate system that significantly increases the computational resources [15]. In addition, it is necessary to have a good solution in order to have a large number of particles that generally demand large memory resources and introduce an additional calculation time for the determination of the trajectories at every time step. A Lagrangian based particle tracking method, versus a Eulerian method, has the advantage of realistic subgrid scale motion and best approximates the movement of particles [16]. However, there are still currently some studies which are carried out by solving the transport equation using a Eulerian approach [17].
Sensibility and accuracy of the Lagrangian Particle Tracking for sediment transport have been well proven. To exemplify, in this research [18] the research first by resolving the hydrodynamic using the nonconservative form of the shallow water equations is carried out. Nonetheless, for the particle movement, a fourth-order Runge-Kutta is used which is computationally heavy. For this reason, OpenMP is used to parallelize computation over memory shared systems [19]. For this reason, some studies are performed avoiding the resolution for ordinary differential equations to lighten the use of computing resources [20]; this is similar to what is proposed in this paper.
For the particle movement, it is considered the specific weight of each particle, as well as the drop velocity [21]. This analysis is very important to calculate the number of particles in suspension. The Lagrangian model has the capacity to be applied after the hydrodynamic fields were calculated or at the same time that the hydrodynamics is simulated. Nevertheless, it is preferable to apply the particle model once the hydrodynamics is estimated for computational resources saving reasons.
Separating the hydrodynamic calculation of the motion simulation of particles allows developing independent transport simulations. Thus, we can simulate particle movement with different locations and sources, various lengths simulations, different transport parameters, and different physical properties of particles (specific weight and diameter), all this supported by the hydrodynamic flow simulation previously calculated. In comparison with other models this is convenient because the hydrodynamic module and the transport module are separated and the erosion and deposition are calculated with particles instead of estimation formulae as part of the hydrodynamic module [1].
Finally, according to some current review papers regarding the Lagrangian modelling of Saltating Sediment Transport [22, 23], a model for the transport of sediment has to include mainly the motion of saltating grains, diffusion of particles, and calculation of bedload transport rate and to improve the motion of particles representing more natural shapes. And, for rivers, it has to consider the nonuniform character of sand, distribution of particle saltation lengths, and excursion lengths may be more important in determining the morphodynamic behavior of the channel bed than the average particle motion. Therefore, future research should focus on relating particle travel distances of individual particles to their shape and physical parameters.
This paper is organized as follows. Section 1 focuses on the introduction. In Section 2 are presented the materials and methods where the hydrodynamic model is introduced, and the description of the particle motion model proposed for the sediment transport. Section 3 shows validation tests of the code proving the advection and the diffusion with different number of particles and are subsequently described; then simulations are carried out in the study area under two dominant scenarios, dry and rainy seasons, and finally the discussion is presented.
2. Materials and Methods
2.1. Hydrodynamic Model
The governing equations for the velocity field are the shallow waters equations, where the model is simplified with the assumption that the vertical scale is smaller than the horizontal scales with the hypothesis of the hydrostatic pressure and considering the Reynolds approximation [2]:where is the coefficient of effective viscosity obtained adding the turbulent velocity coefficient and the molecular velocity coefficient . The coefficient could be approximated as follows [24]:where the vertical longitude is defined as for and for ; is the von Karman constant with typical value of 0.41, is the distance from the river bottom, is the thickness of the boundary layer, and is a constant with a value of 0.09.
2.2. Sediment Transport Model
The numerical model for the particle transport developed in this research is based on a Lagrangian approach (Particle-In-Cell). For this model, the particles are located or placed following an exponential law of concentrations. In addition, it is possible to place the particles randomly in the domain [25].
For the particles movement, it is considered a three-dimensional stochastic model (see Figure 1). For an accuracy simulation, it is necessary to consider the specific weight of each particle, as well as the fall velocity. The movement of the particles in the next time step is calculated from the previous step as follows:where , , and represent the position of the particle in the time () and are the average velocities located in . is the coefficient of turbulent viscosity, is the differential time, and is the sediment fall velocity. The movement equations of the particles 6 are taken from [26] but are modified adding the random movement (Brownian motion).
The particles tracking is governed by the motion equation (4); therefore the particles are subjected to a displacement magnitude of . In either direction of the domain, the sign is positive or negative depending on the direction of the location. In each differential time, the velocity field acts over every particle; therefore the velocity field is the main force that moves the particles. The term represents a positive scalar field, with the turbulent intensity, the shear stress along its deep is expressed as follows [27]:where is the water density, is the kinematic viscosity coefficient, corresponds to the mean flow velocity, and is the double correlation at the bottom. The critical shear stress term acting over the particles is calculated as follows [27]:where is the solid density, the water density, the gravity acceleration, and the diameter of the particles that is greater than the 50% of specific size.
The following equation determines the probability function for the deposition of particles [27]:and the probability function for resuspension of the particles is determined as follows [27]:
To determine the fall velocity of the sediment particles, they are considered to have a nonspherical shape [28]. The form of the particle effects the velocity, particularly relatively large particles (>μm); expressions that determine the magnitude of the fall velocity are given below.where is the diameter of the sieve, is the specific gravity and is the kinematic viscosity coefficient, and the acceleration of gravity.
The particle fall velocity is altered by the presence of other particles [29]; the effect of the other particles causes a greater fall velocity, at higher concentrations. To correct this effect, the following equation is proposed:where is fall velocity in suspension, is the fall velocity of the particle in clean liquid, is the concentration by volume of sediment, and is dimensionless coefficient (with a value from 4 to 5 for particles with diameters between 50 and 500 ).
The discrete-continuous transformation for particle concentration is used to obtain the field of sediment concentration from the Lagrangian modelling of particles [30, 31]. The concentration in a cell is obtained by dividing the total mass of the particles contained in the cell by the volume contained inside.
3. Results
3.1. Particle Model Validation
The model was validated using two cases with known analytic solutions; the first is an instantaneous release of a conservative substance or contaminant in a stagnant environment, in which the mass is moved only by diffusion. The following equation gives the three-dimensional analytical solution [26]:The resulting distribution of concentration was obtained after 50 seconds of real simulation, with initial masses of 1500, 15000, and 50,000 particles with the same diameters and weight. The distribution results at 8 seg, 15 seg, and 50 seg are shown in Figure 2, and the comparison between profiles is depicted in the Figure 3.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
The second validation example corresponds to a continuous release of infinite duration, where the analytical solution corresponds to a conservative one-dimensional substance and is expressed as follows [26]:where is the complementary error function [32]. For this case, 1500, 15000, and 50,000 particles with the same diameters and weight particles were released every 0.05 seconds for 50 seconds. The concentration profiles obtained are depicted in the Figure 4.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
The results obtained from Figures 3 and 4 show that the method used is quite precise and accurate and that keeping the number of particles constant reduces the error. It must be taken into account that increasing indiscriminately the number of particles only increases the use of computational resources.
3.2. Simulation in the Study Zone
The study area is located in the Magdalena River, Colombia (see Figure 5); it is a port terminal that has a dock of approximately 441,417 m2, designed for river operations and small vessels, into two marginal wharves, each with an operational depth of 12 feet. This dock is connected to the Magdalena River, which receives significant sediment supply and presents an exchange of water flow, causing seasonal low flow conditions.
The characterization of total sediment transport in a river current is composed of two types of charge and wash bedding material. The first type of load consists of materials produced by the erosion of the basin, which is characterized by a very fine grain size and fall velocity below the magnitude of the turbulent velocity fluctuations of the current fall, so they are always suspended. The second type of load consists of materials originating from the bed, which in itself is in part moved by entrainment at the river bottom and moves in suspension, depending on availability and flow rate.
The study zone is discretized in a regular mesh (see Figure 6) with a constant spacing in both directions of 10 m ( = = 10 m), giving 320 cells in the direction and 180 cells in the direction.
The simulations were carried out under two scenarios designed in function of the siltation conditions recorded in the entrance of the dock, completely open and partially closed. For this reason, the months of March and October were selected to perform the simulations during 30 days with an increment of time of 0.2 seconds.
The first step is to determine the hydrodynamics of the zone. For this purpose, a well proven numerical hydrodynamic model was used [2, 4, 33, 34] to generate the velocities field for the months of March and October under several initial conditions at the entry of the dock. Flow rates for the years from 2005 to 2013 for the specific months of March and October are mentioned in Table 1.
The mesh used to calculate the hydrodynamic is spatially constant with therefore with 180 discrete elements in the direction and 320 discrete elements in the direction, giving a total of 57600 elements. The results of simulations are depicted in the Figures 7 and 8 for the months of March and October, respectively. The results of the simulations show different behaviors and current patterns that affect the flow exchange between the dock and the river, which consequently affect the dynamics of the sediments.
(a) Open
(b) Partially open
(a) Open
(b) Partially open
The current velocity inside the dock is considerably lower than the flow in the river, being in the order of 3 m/s, being a typical velocity value within these bodies of water. For the dry season (March), the results of the hydrodynamic show that, for both cases (partially open and open dock), the velocities within the basin are relatively low. However, for the scenario when the dock is completely open, there is greater flow exchange between the river and the basin. When the dock is partially open, it is observed that the current of the river forms a homogeneous stream pattern within the basin. While the dock is completely open, the current of the river does not enter directly, but vortices are generated.
For the wet season, the river flow approximately doubles the flow compared with March. Likewise, the velocities are significantly higher, so the flow of the river dominates the behavior and the dynamics within the basin. Additionally, it can be observed in the simulations that the variation of the speed during this month is quite low. These two situations of the wet season (higher velocities) influence significantly the behavior of sediments and suspended material.
Once the simulations of the hydrodynamic conditions for March and October corresponding to the years from 2005 to 2013 are generated, the velocities field is obtained and we can proceed to apply the particle model to estimate the sediment transport.
The dock is clearly influenced by the hydrodynamics of the Magdalena River, due to the important exchanges of water flow. The characteristic diameters of sediment at the bottom were taken from three different sections of the study area. With this information, the sediment model can be set. The diameters are listed in Table 2. The material is 87% sand and the remaining is 13% bulk and volley material.
The initial number of sediment particles at the bottom of the channel and in the inner harbor was established using the PIC method approach. The simulations of the sediment transport were carried out for the two months under study during 720 hours, determining the areas of accumulation of sediment that means the silting.
3.3. Simulation for the Dock Entrance Partially Open
For March, in Figure 9 is shown 24 hrs, 240 hrs, and 720 hrs of simulation with an initial concentration of 395,000 particles distributed in 13 different diameters based on a logarithmic function.
(a) Sediment distribution at 24 hrs
(b) Sediment distribution at 240 hrs
(c) Sediment distribution at 720 hrs
From Figure 9, we can observe constant incipient sediment transport on the river bed, where the blue areas are particles that are deposited randomly depending on the fall velocity and on the velocity gradients in the water column in each cell of the calculation mesh.
In incises (b) and (c) of Figure 9 are presented the simulations of the sediment distribution results after 10 and 31 days, respectively; we estimated the increase in sediment zones. The final simulation results show an average of 75% to 80% of the channel bottom with potential sediment areas or retain a considerable concentration of particles at the bottom.
For the month of October, the river sediment is on average double compared with March. The results show an immediate increase of the magnitude of sediment, and consequently the precipitation of particles tends to decrease towards the bottom. In Figure 10 a reduction in the number of particles in the sediment zones is appreciated.
(a) Sediment distribution at 24 hrs
(b) Sediment distribution at 240 hrs
(c) Sediment distribution at 720 hrs
Similarly, in the incises (b) and (c) of Figure 10, the results for 10 days and 31 days of simulations are shown, respectively. We found that the behavior is similar to the first 24 hours and a decrease of sedimentation is slightly remarkable as compared to March.
3.4. Simulation for the Dock Entrance Completely Open
For the second scenario, when the entrance to the dock is open, for March and October, the results are modified by two issues. In the first instance, bathymetry has slight modifications. This is due to the date when the data were acquired, but the same deposits are observed, and hydrodynamics suffers slight modifications, which significantly alter the behavior of sediment due to cutting with the river bottom.
In Figure 11, the obtained result for March is shown for the first 24 hours of simulation. The same initial particle concentration was used and we can observe that the behavior of the particles in areas where the accumulation occurs produced slight alterations, mainly at increasing the range of percentage concentration around 5%.
(a) Sediment distribution at 24 hrs
(b) Sediment distribution at 240 hrs
(c) Sediment distribution at 720 hrs
In the incises (b) and (c) of Figure 11 are depicted the results for 10 and 31 days of simulation, respectively. The simulation showed an increase of silt deposition; the particles concentration is considerably higher than in the firsts 24 hours, mainly in the dock.
The sediment distribution depicted in the Figure 11 showed that, compared to the closed basin, there is an increase in siltation areas of 15%; this is mainly due to the shape of the river bottom and to the amending of the speed in cells near the river.
The results for October with the dock entrance completely open showed no significant variation. Movement and distribution of bottom particles are practically the same compared with the previous scenario (see Figure 12). The same initial conditions for the semiopen dock were used.
(a) Sediment distribution at 24 hrs
(b) Sediment distribution at 240 hrs
(c) Sediment distribution at 720 hrs
3.5. Quantification of the Results
Quantification of the sediment transport is necessary to validate the simulations results. We studied the months of March and October because they are characteristic of the sediment phenomena in the zone. Therefore, several parameters were recorded, such as, average size of the cross-section, velocity fields, and granulometry.
In Table 3 the total particle count for each cell of the domain is presented, obtaining an approximation of the solid sediment for each month under study.
The method used to corroborate the quantification of sediment transport was developed by Van Rijn [27] and amended by Maza and García [35]. Table 4 shows the comparisons of results obtained with the simulations for the open dock conditions versus the obtained with the method of Van Rijn.
Finally, the quantification of the number of particles in each cell at the entrance of the dock is depicted in Figures 13 and 14, in which sedimentation profiles for the two months are calculated.
4. Discussion
The model has been validated with a series of simple analytical cases showing a good agreement and has been used to simulated a real case of sediment transport in a section of Magdalena river in Colombia. The introduction of random motion to the equations of the particle movement produces good results by simplifying the equations of particle motion, thus producing a model that can withstand a large number of particles without consuming too many computational resources.
The sediment transport model developed for the simulation of the sediment transport in the study zone is robust and precise; considering the river dimensions, flow rates, sediment deposits, and the flow contribution of the left branch of the river, the analysis and simulation are a challenge.
The sediment displacement is calculated for March and October, because these are characteristic months when the pier is open or closed in a similar magnitude, 0.003 Kg/s for March and 0.155 Kg/s for October.
The distributions of areas of siltation or concentration of sediment particles in the river bottom are variables, both between stages and between months, giving an average of potential sediment areas in the riverbed of 75% for March and of 60% for October.
Regarding the behavior of the dock, sediment transport considerably alters the shape of the river bottom, 44% of sedimentation occurs in March (involving an open and closed entrance), caused mainly by decreasing the volume of incoming and outgoing water between the dock and the river itself. However, for the month of October, the balance of incoming and outgoing volume improves drastically. On average 30% of the total surface has potential siltation, which is distributed at the entrance and in the central area.
Therefore, it may be concluded that sediment transport allows natural closure of the entrance mainly in the dry season, recovering an important area when it rains and, therefore, the increase in sediment deposit and speed of the river.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Authors’ Contributions
Israel E. Herrera-Díaz, Franklin M. Torres-Bejarano, Jatziri Y. Moreno-Martínez, and C. Rodriguez-Cuevas contributed equally to this work.
Acknowledgments
This study was carried out under the financial support of the University of Guanajuato, Mexico. C. Couder-Castañeda acknowledges EDI Grant by SIP/IPN. This research was partially supported by Projects SIP/IPN nos. 20160105 and SIP 20160576.