Assessment of the Impact of Sand Mining on Bottom Morphology in the Mekong River in An Giang Province, Vietnam, Using a Hydro-Morphological Model with GPU Computing

: Sand mining, among the many activities that have signiﬁcant e ﬀ ects on the bed changes of rivers, has increased in many parts of the world in recent decades. Numerical modeling plays a vital role in simulation in the long term; however, computational time remains a challenge. In this paper, we propose a sand mining component integrated into the bedload continuity equation and combine it with high-performance computing using graphics processing units to boost the speed of the simulation. The developed numerical model is applied to the Mekong river segment, ﬂowing through Tan Chau Town, An Giang Province, Vietnam. The 20 years from 1999 to 2019 is examined in this study, both with and without sand mining activities. The results show that the numerical model can simulate the bed change for the period from 1999 to 2019. By adding the sand mining component (2002–2006), the bed change in the river is modeled closely after the actual development. The Tan An sand mine in the area (2002–2006) caused the channel to deviate slightly from that of An Giang and created a slight erosion channel in 2006 ( − 23 m). From 2006 to 2014, although Tan An mine stopped operating, the riverbed recovered quite slowly with a small accretion rate (0.25 m / year). However, the Tan An sand mine eroded again from 2014–2019 due to a lack of sand. In 2014, in the Vinh Hoa communes, An Giang Province, the Vinh Hoa sand mine began to operate. The results of simulating with sand mining incidents proved that sand mining caused the erosion channel to move towards the sand mines, and the erosion speed was faster when there was no sand mining. Combined with high-performance computing, harnessing the power of


Introduction
Sand is an essential resource in the construction and maintenance of river delta surfaces. In the process of increasing urbanization, the demand for extensive construction material is becoming a concern for infrastructure. To meet the demands for construction material, several sand mines on the river have been operating and have produced many impacts on the environment and river morphology [1][2][3][4][5]. Numerous studies have demonstrated that sand mining in alluvial rivers causes disturbances within geomorphic processes, bottom revolutions, and their functions [2,4,6]. Removing a large volume of sediment alters the hydraulic properties of the flow and the morphodynamic evolution of the riverbed [6][7][8][9][10].
The Mekong River is the 10th largest river in the world, with a basin of about 795,000 km 2 in size and about 4800 km in length. The main portion of the Mekong River flows through six countries including China, Myanmar, Thailand, Laos, Cambodia, and Vietnam. [11]. The lower Mekong morphology is created by bedrock-controlled and alluvial reaches [12]. There are many causes of continuous and severe erosion in the area, such as excessive sand mining and sediment imbalances leading to erosion [13,14]. The annual extraction, according to the estimation of Bravard et al. (2013), is 34 Mm 3 (approximately 54 MT/year for an aggregate density of 1.6 t/m 3 ), featuring mostly (90%) sand and, to a smaller degree, gravel (8%) and pebbles (1%). This utilization is focused on mining in Cambodia (21 MT/year), Vietnam (8 MT/year), Thailand (5 MT/year), and Laos (1 MT/year) [15]. It has been estimated that sand mining has removed 200 Mm 3 of bedload since 1998. Hackney's research indicated that sand in the Mekong riverbed is minimal and that sand transport to the delta is under threat from increased bed sediment extraction [16]. The total sand flux entering the Mekong delta (6.17 MT/year ± 2.01 MT/year) is far less than the current sand extraction rates (50 MT/year). As a result, it takes only 10 to 18 years to scour the riverbed sufficiently at these current rates to change the riverbed and sediment transport [17].
In a geological survey study along the Tien River, the characteristics of the Tien Riverbed in Hong Ngu, An Long, Sa Dec, Mo Cay, Ben Tre, and Ba Tri were shown to be clay (d < 0.0015 mm), while Tan Chau and My Thuan were characterized by clay mud (0.0015 < d < 0.003 mm) [18]. Compared with the Hjulstrom-Sundborg diagram, a flow with a velocity higher than 0.005 m/s is capable of carrying moving particles from 0.002 to 0.075 mm, and if the speed is higher than 0.3 m/s-0.4 m/s, the sediment particles with the above dimensions are likely to be separated from the bottom by the flow [19]. Currently, 11 sand mines are operating in An Giang area according to Decision No. 1697/QÐ-UBND on 18 July 2018, An Giang Province. The voids created by the mining (the pits themselves and the incised channels) trap sediment transported into the reach from upstream, reducing sediment loads downstream of the pit, thereby inducing incision from the flow. This matter leads to bed changes in the river [20].
Many studies have used numerical models to simulate sediment transport and bed changes, including MIKE ( [21,22], DELFT [23], TELEMAC [24,25], CCHE 2,3D [26,27]. In general, these studies focus on explaining the current situation. However, for the effects under the scenarios of flow and sediment changes, there is no research demonstrating the impact of sand mining on the bottom changes in this area. Moreover, in these studies, the hydrodynamic model is the primary tool. A model was developed by our team and successfully applied to several other regions, including the Mekong Delta. Bay et al. (2018) simulated the litho-dynamic processes and bottom morphology in the coastal area, such as the flow, sediment transport, and bed changes for Can Gio and Cua Lap, Vietnam [28,29]. This model is based on the Reynolds equation system coupled with the sediment transport and bedload continuity equations [6,10,[30][31][32][33][34]. The authors applied the above model to calculate the bottom changes for the waters of Can Gio and Cua Lap. The results from the model prove that the calculations and developments are suitable in reality. In 2019, the authors improved the hydraulic model with a roughness coefficient that changes with the bed level. Moreover, the alluvial boundary condition at the liquid boundary in the sediment transport module was developed using a characteristic line, allowing sediment to move out of the computational area. Thus, the movements of suspended sediments and the bottom transport in the near-boundary area are more accurate. The study area in this research was a part of the Tien River that flows through the Tan Chau area of the Mekong Delta [35]. However, the results were still inaccurate compared to reality due to sand exploitation in the area.
For bed change due to sand mining, Binoy Aliyas Mattamana (2013) studied the transportation of sand in the river using numerical modeling and applied sand inflow to Periyar River, India. The present study deals with Meyer-Peter's equations to estimate bedload transport. However, this study only focused on calculating the sand flow to the sand mine area to assess the recovery time of the sand mine. The amount of sand lost due to exploitation and riverbed evolution during sand mining has not been studied [36]. Other sand mining studies also focused on assessing the impact of the sand mining concept on river hydrodynamics [37] and environmental water quality [38,39].
In another respect, although simulations of the morphological changes and other hydrodynamic phenomena can be modeled via partial differential equations, this method is implicitly impractical due to the large temporal and spatial grid. Hence, it is not possible to solve within a certain amount of time [40][41][42][43][44][45]. In the 15th Pacific Conference on Computer Graphics and Applications, the United States, Wu and Eftekkarian (2011) emphasized that reducing the calculation speed and improving in-situ visualization will be necessary for integration hydrodynamic models [45][46][47].
Our primary motivation for this improvement is that in the current era of high-performance computing, numerous multi-core co-processors can help accelerate simulation processes at different scales (Intel MIC, Graphics Processing Units) in which the data are laid out as a rectilinear grid. There is little or no dependency between the results of the current and other time steps. At present, increasingly more artificially intelligent applications whose core components are neural networks are using graphics processing units (GPUs) to accelerate their training and inference phases. Specifically, we developed our numerical model using NVIDIA Compute Unified Device Architecture (CUDA) C/C++ to reduce the timing execution of the simulation. This library helps our code communicate with the GPUs' parallel processing units (physical threads) via the abstraction of CUDA logical threads.
Based on the above analysis, in this paper, we proposed a sand mining component for the numerical model established by Bay et al. [35] to achieve the most accurate bottom evolutions. Furthermore, the numerical model is combined with high-performance computing using graphics processing units. The objectives of this study are as follows: (1) First, the computational times on a single-core CPU machine are too slow to meet the simulation's on-demand requirements. We thus parallelize our solver using massively multi-core GPUs to accelerate the computing speed, which is evaluated with a strong-scaling factor. (2) Second, we assess the changes in bottom morphology in the Mekong River in An Giang province under sand mining.

Numerical Model
This model is based on a system of four governing equations averaged according to depth as follows. In particular, the ς-axis is facing up, and the standard "0" is placed on the static surface, as shown in Figure 2.

Numerical Model
This model is based on a system of four governing equations averaged according to depth as follows. In particular, the ς-axis is facing up, and the standard "0" is placed on the static surface, as shown in Figure 2.
Water 2020, 12, x FOR PEER REVIEW 4 of 25 In this area, there are two sand mines, including Tan An (2002)(2003)(2004)(2005)(2006) and Vinh Hoa (2014-2019). The Tan An sand mine was invested in by Vinh An Private Enterprise (An Giang). This mine operated from 2002 to 2006; its annual output is 200,000 m 3 , and its mining area is 91 ha [49]. From 2014 until the present, the Vinh Hoa mine has operated (According to License No. 04/GP-UBND 23 April 2013) with an annual output of 200,000 m 3 and a mining area of 27.9 ha.

Numerical Model
This model is based on a system of four governing equations averaged according to depth as follows. In particular, the ς-axis is facing up, and the standard "0" is placed on the static surface, as shown in Figure 2.
The model adopted is a 2D surface model where Ox and Oy represent the study area's length and width.

Suspended Sediment Transport Equation
where S in this equation stands for the erosion and deposition rates (m/s). The source S is defined after Van Rijn [51] with three different states: When τ b > τ e , S = E (Erosion rate): When τ b < τ d , S = D (Deposition rate): When τ e ≥ τ b ≥ τ d : Here, S sm is the sand mine source and the standing sand mining rate (m/s). In this paper, the sand component is integrated into the bed load continuity equation. After sand mining, sand is removed from the calculation area; hence, S sm is not involved in the suspended sediment transport equation. Sand mining will stop when the bottom elevation reaches the government's permitted mining level, where q b = (q bx , q by ).
These equations are solved by the Alternating Direction Implicit (ADI) method. The ADI method's central concept is to split different finite equations into two: one with the x-derivative and the other with the y-derivative, both taken implicitly [52].

Graphics Processing Unit (GPU)
In the current era of high-performance computing, numerous multi-core co-processors can help accelerate the simulation processes at different scales (Intel MIC, Graphics Processing Units) where the data are laid out as a rectilinear grid. There is little or no dependency between the results of the current and other time steps. Presently, increasingly more artificially intelligent applications whose core components are neural networks make use of GPUs to accelerate their training and inference phases. Programming such accelerators requires a medium to interact with the physically parallel cores that handle the data. For example, Intel MIC and ATI GPUs can be accelerated by Open Computing Language (OpenCL); NVIDIA GPUs ship with NVIDIA CUDA C/C++. In particular, these libraries/programming languages help our code communicate with the parallel processing units of GPUs (physical threads) via the abstraction of logical threads. Running our numerical model on GPUs can reduce the timing execution of simulation. We chose to implement our model on NVIDIA GPUs since these GPUs achieve higher Floating Point Operation Per Second (FLOP) when dealing with numerically precise data. NVIDIA GPU schedulers help map many logical threads to many (but limited) physical threads with zero-cost. In other words, switching between the connections of logical and physical threads will not waste time and energy. In NVIDIA GPUs, the logical threads are laid out in warps, blocks, and grids (in a coarse-to-fine manner). Each warp contains 32 threads, and each block can hold a maximum of 1024 threads. All of the threads in the same block are scheduled to be executed concurrently, and all 32 threads in the same warp are designed to run a set of instructions simultaneously. This is why GPUs can support running a program in parallel with many computational operations (i.e., an embarrassingly parallel workload).
A computational algorithm such as the above numerical model is a procedure in which each step is a correctly defined state and executed by computers. For a given problem, multiple algorithms are commonly used to solve it. Some may require a lower computational cost than others. Some may allow execution in parallel at a larger scale than others; some are more stable, and some require less memory consumption. Unfortunately, no solution satisfies all of the above four criteria. Given a problem and a list of constraints to run that problem, we usually have to choose a direction to compensate for the above four aspects under a specific physical hardware system. Some algorithms have certain advantages/disadvantages in maintaining certain numerical precisions. In contrast, others have to sacrifice accuracy to scale up the solution to handle big data such as our hydrological data.
The theoretical solution for this 2D hydrological problem is common to use Finite Difference Methods (FDM) with the ADI scheme. At the first half of the time step, our computational model starts to scan the grid along the horizon to compute the velocity u and the depth ς implicitly, and then explicitly solves for the velocity ν. Similarly, at the second half of the time step, the model scans the grid along the vertical axis to compute the velocity ν and the depth ς implicitly. Explicitly, the model applies the FDM to numericalize the equation and rearrange variables; it will then form a system of linear equations (more concretely, a tridiagonal linear system). Some numerical solvers for this, such as the Thomas algorithm [53], include a forward scan and a backward substitution. However, there are more dependencies in this algorithm that hinder its full parallelization. Therefore, leveraging the Thomas algorithm for solving tridiagonal linear systems for our problem is simple but cannot exploit the massive parallelism of a GPU or a GPU's high memory bandwidth. Various algorithms have been developed to parallelize banded matrix solvers [53]. After examining different parallel matrix solvers, we decided to adopt the parallelized Thomas algorithm with data marshaling, developed in [54] with open source code provided in [55]. To avoid memory transfer from the device to the host, we exploit a CUDA feature called dynamic parallelism, where the device kernels can be launched from other kernels. By parallelizing the Thomas algorithm and performing data marshaling, we achieved a higher level of parallelism and higher memory bandwidth, increasing the performance significantly.
The Alternating Direction Implicit (ADI) scheme is known to be naturally suitable for parallelization, since, in each sweep, each row/column can be scanned independently. Therefore, each row or column can be assigned to one thread to be solved individually. However, due to the GPU's hierarchical memory, the computation needs to be distributed into small chunks to process the tridiagonal linear system. Hence, there are many more optimization steps in programming design strategies to access different levels of GPU memories, such as the local registers, shared memory, cache, and global memory. By exploiting these, we reduced the computational overhead from the data movement in between the GPU memories and improved the performance of the GPU cores.

Study Area Mesh
The topography was collected from the Department of Investment and Construction Project of Tan Chau town area on 6 October 1999 at The Southern Institute Of Water Resources Research in 2006 and 2014 and from the "Research to identify causes, mechanisms and propose feasible technical and economical solutions to reduce erosion, sedimentation for the Mekong river system (2017-2020)" project, code No. KHCN-TNB.ÐT/14-19/C10, in 2019.
The features (water level ς(t), discharge Q(t), and total suspended sediment C(t)) at Tan Chau station were collected from 1999 to 2019 ( Figure 1). The study area has two boundaries whose features were identified based on their correlation function (f(ς), f(Q), f(C) with Tan Chau station. Measurement data were taken from the boundaries of the "Research to identify causes, mechanisms and propose feasible technical and economical solutions to reduce erosion, sedimentation for the Mekong river system (2017-2020)" project, code No. KHCN-TNB.ÐT/14-19/C10, from 10:00on 6 June 2018 to 10:00 on 13 June 2018.
The study area mesh is a rectangular grid with 981 × 845 elements, including the land and riverbed, with ∆x = ∆y = 10 m for the whole region. Sand mine sources were added to the model by matrices S sm (i, j) (as Figure 3), in which S sm (i, j) is the sand mining rate at i, j cells (the locations of these sand mines are described in Figure 1). The sand mining rate at the i, j cell corresponding to scenarios (SC1, SC2A, SC2B, SC3, SC4A, SC4B), the annual output, the square of Tan An, and the Vinh Hoa sand mines, is outlined in Table 1.
The speed and location of the sand mining of Tan Chau and Vinh Hoa sand mines are shown in Figure 3: GPU's hierarchical memory, the computation needs to be distributed into small chunks to process the tridiagonal linear system. Hence, there are many more optimization steps in programming design strategies to access different levels of GPU memories, such as the local registers, shared memory, cache, and global memory. By exploiting these, we reduced the computational overhead from the data movement in between the GPU memories and improved the performance of the GPU cores.

Study Area Mesh
The topography was collected from the Department of Investment and Construction Project of Tan Chau town area on 6 October 1999 at The Southern Institute Of Water Resources Research in 2006 and 2014 and from the "Research to identify causes, mechanisms and propose feasible technical and economical solutions to reduce erosion, sedimentation for the Mekong river system (2017-2020)" project, code No. KHCN-TNB.ĐT/14-19/C10, in 2019.
The features (water level ς(t), discharge Q(t), and total suspended sediment C(t)) at Tan Chau station were collected from 1999 to 2019 ( Figure 1). The study area has two boundaries whose features were identified based on their correlation function (f(ς), f(Q), f(C) with Tan Chau station. Measurement data were taken from the boundaries of the "Research to identify causes, mechanisms and propose feasible technical and economical solutions to reduce erosion, sedimentation for the Mekong river system (2017-2020)" project, code No. KHCN-TNB.ĐT/14-19/C10, from 10:00on 6 June 2018 to 10:00 on 13 June 2018.
The study area mesh is a rectangular grid with 981 × 845 elements, including the land and riverbed, with Δx = Δy = 10 m for the whole region. Sand mine sources were added to the model by matrices Ssm(i, j) (as Figure 3), in which Ssm(i, j) is the sand mining rate at i, j cells (the locations of these sand mines are described in Figure 1). The sand mining rate at the i, j cell corresponding to scenarios (SC1, SC2A, SC2B, SC3, SC4A, SC4B), the annual output, the square of Tan An, and the Vinh Hoa sand mines, is outlined in Table 1.
The speed and location of the sand mining of Tan Chau and Vinh Hoa sand mines are shown in Figure 3:

Initial Conditions
In the model, if starting from t 0 = 0, the hydraulic module is tied to a static state (Figure 2), and the sediment transport module is set to an initial constant basal concentration. Where the problem is calculated from a time t = t 1 , the initial condition will be the velocity fields u, ν (x, y), and concentration C (x, y) at time t 1 across the computational domain.

Boundary Conditions
The Open Boundary

Hydrodynamics conditions:
The upstream boundary is Q(t), from Q(t), the velocity recalculated as follows [56]: Furthermore, because up and down the fluctuations of the water levels vary from time to time (and due to the moving boundaries problem (flooding and drying fronts)), the study area is classified into calculation grid cells. The depth in each element/cell was monitored, and the elements were classified as dry, partially dry, or wet.

Sediment transport conditions:
When water flows into the computational domain, C = C o . In cases where the flow flows out of the domain, C is calculated through the advective transport process (where the diffusion process is ignored), solved by the characteristic line method [35]. By this method, sediment is loaded in and out of the calculation area, making the simulation results more accurate.

Solid Wall Boundary
Hydrodynamics conditions: u n = 0 Sediment transport conditions: ∂C ∂n p = 0 For the calibration and validation model, the location of the calibration is Tan Chau station. In this model, the calibration time is 1999, and 2018 (7-day measurement data from 10:00 6 June 2018 to 10:00 13 June 2018) is the validation time. The input data, including flow Q(t), water level ς(t), and sediment concentration C(t) for 1999 and 2018, are shown in Figure 4. The sediment boundary downstream is not used for the model because it flows out of the domain.

Hydraulic Model
To calibrate and validate the models, as well as for comparison purposes, the Nash-Sutcliffe efficiency coefficient (NSE) and the coefficient of determination (R 2 ) were used to measure the model performance [57][58][59].
The verification parameter in this model is the roughness coefficient, which was changed in inverse proportion to the water depth. The calibration results show that there was little difference between the simulation results and measurements of the discharge and water level at Tan Chau station between the observation and simulation results from 11 July 2002 to 30 July 2002 ( Figure 5). The graphical results during calibration indicated adequate calibration and validation across the range of discharge and water levels. However, the calibration results of the water level showed a better match than those of the discharge. The NSE values for the calibration of discharge and the water level reached 0.61 and 0.77, respectively, while the R 2 values were 0.98 and 0.99. According to Moriasi's research [59], these values indicate that the hydraulic model performance achieved outstanding values. Thus, it is believed that the hydrodynamic model is well-calibrated and that the predicted results are close to actual water movements.
Water 2020, 12, x FOR PEER REVIEW 11 of 25

Hydraulic Model
To calibrate and validate the models, as well as for comparison purposes, the Nash-Sutcliffe efficiency coefficient (NSE) and the coefficient of determination (R 2 ) were used to measure the model performance [57][58][59].
The verification parameter in this model is the roughness coefficient, which was changed in inverse proportion to the water depth. The calibration results show that there was little difference between the simulation results and measurements of the discharge and water level at Tan Chau station between the observation and simulation results from 11 July 2002 to 30 July 2002 ( Figure 5). The graphical results during calibration indicated adequate calibration and validation across the range of discharge and water levels. However, the calibration results of the water level showed a better match than those of the discharge. The NSE values for the calibration of discharge and the water level reached 0.61 and 0.77, respectively, while the R 2 values were 0.98 and 0.99. According to Moriasi's research [59], these values indicate that the hydraulic model performance achieved outstanding values. Thus, it is believed that the hydrodynamic model is well-calibrated and that the predicted results are close to actual water movements. The roughness coefficient was recorded after calibration with a range of 0.055 from 0.005 to 0.06 corresponding to the bed level, varying from 41 to 0.1 m. During the experimental process, n2 (0.003) corresponding to h 2 (25 m) was found to make the roughness coefficient in the study area more suitable. Hence, the roughness value is changed from 0.005 to 0.03 when the bed level ranges from 0.1 to 25 m and from 0.003 to 0.06 when the bed level ranges from 25 to 41 m.
Below is a rough map for the first 500 h starting on 15 May 1999 in Figure 6. Consequently, the roughness coefficient in the domain was computed, as shown based on the bed level at that time. The roughness coefficient was recorded after calibration with a range of 0.055 from 0.005 to 0.06 corresponding to the bed level, varying from 41 to 0.1 m. During the experimental process, n 2 (0.003) corresponding to h 2 (25 m) was found to make the roughness coefficient in the study area more suitable. Hence, the roughness value is changed from 0.005 to 0.03 when the bed level ranges from 0.1 to 25 m and from 0.003 to 0.06 when the bed level ranges from 25 to 41 m.
Below is a rough map for the first 500 h starting on 15 May 1999 in Figure 6. Consequently, the roughness coefficient in the domain was computed, as shown based on the bed level at that time.
After calibration, the model was validated against 7-day measurement data from 10:00 6 June 2018 to 10:00 13 June 2018. Similar to the calibration, although the calibration result of the water level showed a better match than those of discharge (Figure 7), the evaluation criteria values were also very good for both features. While NSE values were 0.92 and 0.98 for discharge and water level, the R 2 values were 0.99 for both these features. Water 2020, 12, x FOR PEER REVIEW 12 of 25 After calibration, the model was validated against 7-day measurement data from 10:00 6 June 2018 to 10:00 13 June 2018. Similar to the calibration, although the calibration result of the water level showed a better match than those of discharge (Figure 7), the evaluation criteria values were also very good for both features. While NSE values were 0.92 and 0.98 for discharge and water level, the R 2 values were 0.99 for both these features.   After calibration, the model was validated against 7-day measurement data from 10:00 6 June 2018 to 10:00 13 June 2018. Similar to the calibration, although the calibration result of the water level showed a better match than those of discharge (Figure 7), the evaluation criteria values were also very good for both features. While NSE values were 0.92 and 0.98 for discharge and water level, the R 2 values were 0.99 for both these features.   After calibration, the model was validated against 7-day measurement data from 10:00 6 June 2018 to 10:00 13 June 2018. Similar to the calibration, although the calibration result of the water level showed a better match than those of discharge (Figure 7), the evaluation criteria values were also very good for both features. While NSE values were 0.92 and 0.98 for discharge and water level, the R 2 values were 0.99 for both these features.

Sediment Transport Model
Similar to the hydraulic model, the sediment transport simulation was calibrated and validated for the two periods from 11 July 2002 to 30 October 2002 and for the 7-day measurement data from 10:00 6 June 2018 to 10:00 13 June 2018. According to evaluation criteria values, the sediment transport model achieved excellent results in terms of the sediment concentration for calibration and validation, with 0.67 and 0.52 for the NSE values and 0.96 and 0.81 R 2 values. The graphical results during calibration and validation at Tan Chau station are shown in Figure 8. After calibration and validation, numerous studies have applied bed parameters such as dispersion, critical shear stress for deposition, and acute shear stress for the erosion of bed layers to simulate the processes of sediment transportation, erosion, and deposition [60]. In this study, these parameters were used for calibrating the sediment simulation. The most relevant parameters in the model are summarized in Table 2.

Improved Computing Speed When Combining GPUs
We conducted the experiment by directly comparing the accuracy, running times, and the scalability of the proposed solution with and without GPUs. Our computational workstation was as follows: Intel i7 7000 CPU, 16 GB RAM, equipped with an NVIDIA Titan V GPU 12 GB VRAM. Consequently, running in the same environment will produce a direct and fair comparison and a concise result for acceleration (Figures 9 and 10). After calibration and validation, numerous studies have applied bed parameters such as dispersion, critical shear stress for deposition, and acute shear stress for the erosion of bed layers to simulate the processes of sediment transportation, erosion, and deposition [60]. In this study, these parameters were used for calibrating the sediment simulation. The most relevant parameters in the model are summarized in Table 2.

Improved Computing Speed When Combining GPUs
We conducted the experiment by directly comparing the accuracy, running times, and the scalability of the proposed solution with and without GPUs. Our computational workstation was as follows: Intel i7 7000 CPU, 16 GB RAM, equipped with an NVIDIA Titan V GPU 12 GB VRAM. Consequently, running in the same environment will produce a direct and fair comparison and a concise result for acceleration (Figures 9 and 10). As can be seen from Figure 11, the results obtained from the CPU (drawn by Surfer) and GPU versions (drawn by an in-situ Python script) align with each other and thus guarantee the correctness of the model (the source of error from numerical precision is considered insignificant). Moreover, the running times of the simulations (Table 3) show that harnessing the power of GPUs can accelerate performance from 10 to 25 times.  As can be seen from Figure 11, the results obtained from the CPU (drawn by Surfer) and GPU versions (drawn by an in-situ Python script) align with each other and thus guarantee the correctness of the model (the source of error from numerical precision is considered insignificant). Moreover, the running times of the simulations (Table 3) show that harnessing the power of GPUs can accelerate performance from 10 to 25 times. As can be seen from Figure 11, the results obtained from the CPU (drawn by Surfer) and GPU versions (drawn by an in-situ Python script) align with each other and thus guarantee the correctness of the model (the source of error from numerical precision is considered insignificant). Moreover, the running times of the simulations (Table 3) show that harnessing the power of GPUs can accelerate performance from 10 to 25 times.  In this subsection, we explain and discuss the methodology that parallelizes the new programming language of CUDA C/C++ using GPUs to accelerate the simulation of 2D hydrological models. We chose and combined multiple CPUs approaches and usedc current state-of-the-art tridiagonal linear system solvers on the GPU to achieve the same results within a shorter period of time in both coarse and fine grid decomposition. The results clearly show that leveraging GPU acceleration can accommodate complex grids in a large-scale manner. Table 3. A direct comparison of running times on the CPU (a) and GPU (b) for 20 m-and 10 m-grid simulations.

The Hydraulic Simulation of the Tien River Segment in Mekong Delta
The simulation period for this study was continuous from 1999 to 2019. A velocity field of the peaking flood at 21:00 on 29 September 2002 was extracted, as shown in Figure 12.  In this subsection, we explain and discuss the methodology that parallelizes the new programming language of CUDA C/C++ using GPUs to accelerate the simulation of 2D hydrological models. We chose and combined multiple CPUs approaches and usedc current state-of-the-art tridiagonal linear system solvers on the GPU to achieve the same results within a shorter period of time in both coarse and fine grid decomposition. The results clearly show that leveraging GPU acceleration can accommodate complex grids in a large-scale manner.

The Hydraulic Simulation of the Tien River Segment in Mekong Delta
The simulation period for this study was continuous from 1999 to 2019. A velocity field of the peaking flood at 21:00 on 29 September 2002 was extracted, as shown in Figure 12.  The velocity field results of the peaking flood at 21:00 on 29 September 2002 show that the following: At site 1 (Figure 12), due to a combination of initially deep terrain (−12 m) and the suddenly narrowing section compared to the upstream section, the velocity of the flow was quite large (about 3 m/s). The conflux flow after passing through Chinh Sach islet tended to push towards An Giang (about 2.3 m/s, at site 2). The velocity at the position of the policy tail (site 3) was relatively small (0.13 m/s), which is smaller than the deposition speed of the sediment particles at the Tien river segment (about 0.34 m/s) [19], so this position tended to accrete.
Similarly, at the tail of the mudflat on the Dong Thap side (at site 4), the bathymetry was quite shallow (about −1 m). When the flow velocity was small (0.12 m/s), the sediment particle was deposited gradually, leading to the mudflat being widened.
At the curved section (site 5), there was a combination of a narrowed curved flow and a deep erosion pit, so the flow velocity here was very large (about 3.56 m/s) leading to severe bed erosion.

Results of Bottom Changes and Analysis
According to the bottom changes in the study area calculated from 1999 to 2019, four scenarios happened in reality, as follows:  Figure 13b). Meanwhile, the two banks were gradually deposited. The islet tended to erode upstream and settle downstream gradually. The erosion channel on the Dong Thap side moved gradually towards An Giang.
At the cross-section 1.1 (Figure 13d), the initial bottom depth of 1999 at the deepest position was −18 m (near the riverbank of Dong Thap). However, by 2002, it moved towards An Giang (the average rate was 80 m/year), reaching a value of −16 m. Meanwhile, close to the banks of Dong Thap and An Giang, the simulation results show that there was an accretion phenomenon, with an average sedimentation speed of about 0.4 m/year (see Figure 13c). In this area, in 2002, the People's Committee of An Giang province allowed the Tan An sand mine to begin operation [49].  It can be seen that sand mining was the cause of the displacement of the erosion channel. The channel flowing through the sand mine position and the lack of sand made this area further eroded. Sand mining made the erosion region widen compared to the scenarios without sand exploitation, which was the predominant erosion trend for sand extraction (see Figure 14e,f). The results show that the riverbed morphology from the sand mine position to the downstream changed between the two scenarios, with no change observed towards the upstream.
Water 2020, 12, x FOR PEER REVIEW 18 of 25 It can be seen that sand mining was the cause of the displacement of the erosion channel. The channel flowing through the sand mine position and the lack of sand made this area further eroded. Sand mining made the erosion region widen compared to the scenarios without sand exploitation, which was the predominant erosion trend for sand extraction (see Figures 14e,f). The results show that the riverbed morphology from the sand mine position to the downstream changed between the two scenarios, with no change observed towards the upstream. During this period, there was no sand mining in the study area (see Table 1). According to the amount of alluvium flowing from the upstream area of the study area, the erosion channel in the mining area (Tan An) gradually accreted, and the average speed reached 0.25 m/year (see Figures  15d,e). Compared to when the sand mines were operating, the sedimentation speed in this area was relatively slow. It should also be noted that from 2002 to 2006, the erosion speed here reached nearly 3 m/year.
The results show that the accretion process prevailed over a period of 8 years from 2006 to 2014. Stream erosion tended to move back toward Dong Thap (Figure 15e), and the islet was deposited towards the downstream (see Figure 15d). During this period, there was no sand mining in the study area (see Table 1). According to the amount of alluvium flowing from the upstream area of the study area, the erosion channel in the mining area (Tan An) gradually accreted, and the average speed reached 0.25 m/year (see Figure 15d,e). Compared to when the sand mines were operating, the sedimentation speed in this area was relatively slow. It should also be noted that from 2002 to 2006, the erosion speed here reached nearly 3 m/year. SC4A: During this period, alluvial sediment was reduced by two-thirds compared to the previous period [14], and, accordingly, the sediment boundary was also reduced. The calculation results show that the whole study area tended to erode, even at the site of Tan An mine (which has stopped operating), which gradually eroded (see Figures 16f,h).
The results were compared between the observations and simulations at cross-section 2.2, which recorded a small error in 2019 (Figure 16g). The sand mining in this area made stream erosion flow towards Vinh Hoa sand mine and connected to the stream erosion downstream (Figures 16c,h).
At the cross-section 1.1, the riverbed eroded from a depth of −20 to about −22 m (the speed erosion of 0.4 m/year). Moreover, in the Vinh Hoa communes, An Giang Province, the Vinh Hoa sand mine was in operation (see Table 1). The calculation results show that there was an eroded channel here (the bed elevation changed from −12 (2014) to −18 m (2019)), with an erosion speed of 1.2 m/year. The bed change here caused the flow to change its path towards Dong Thap. SC4A: During this period, alluvial sediment was reduced by two-thirds compared to the previous period [14], and, accordingly, the sediment boundary was also reduced. The calculation results show that the whole study area tended to erode, even at the site of Tan An mine (which has stopped operating), which gradually eroded (see Figure 16f,h).
The results were compared between the observations and simulations at cross-section 2.2, which recorded a small error in 2019 (Figure 16g). The sand mining in this area made stream erosion flow towards Vinh Hoa sand mine and connected to the stream erosion downstream (Figure 16c,h).
At the cross-section 1.1, the riverbed eroded from a depth of −20 to about −22 m (the speed erosion of 0.4 m/year). Moreover, in the Vinh Hoa communes, An Giang Province, the Vinh Hoa sand mine was in operation (see Table 1). The calculation results show that there was an eroded channel here (the bed elevation changed from −12 (2014) to −18 m (2019)), with an erosion speed of 1.2 m/year. The bed change here caused the flow to change its path towards Dong Thap. In the case in An Giang, the erosion canal moved closer to shore (the sand mining site was near the shore). The erosion channel near the riverbank causes the slope of the riverbank in An Giang to be steeper. As a result, the stability of the bank decreases; this is one of the causes of a river bank collapse. This matter directly affects the socio-economic development planning of the local area and construction along the river in An Giang Province. although Tan An mine was not in operation, the riverbed in this place recovered quite slowly with a small accretion rate (0.25 m/year). The operation of the Vinh Hoa sand mine has created a deep erosion channel here, and made the main flow move completely over Vinh Hoa sand mine during the period of 2014-2019. Moreover, a decrease in the amount of alluvium caused Tan An sand mine to erode again (the average rate of erosion was 0.4 m/year). Sand mining made the erosion channel move towards sand mines, with a faster erosion speed in the absence of sand mining. In the case in An Giang, the erosion canal moved closer to shore (the sand mining site was near the shore). The erosion channel near the riverbank causes the slope of the riverbank in An Giang to be steeper. As a result, the stability of the bank decreases; this is one of the causes of a river bank collapse. This matter directly affects the socio-economic development planning of the local area and construction along the river in An Giang Province.

Conclusions
The results of simulating the sedimentation progress of the study area using the model were appropriate, compared to the measurement bathymetry of specific years (2006, 2014, and 2019) when Tan An and Vinh Hoa sand mines were added. Due to the impact of sand mining, the main flow changed, producing a significant effect on the hydrodynamic regime, as well as the morphology of the Water 2020, 12, 2912 20 of 24 riverbed in the study area. The analysis showed that sand mining made the erosion channel move towards the sand mines, after which the erosion speed increased when there was no sand mining.
The simulation period was quite long (1999-2019), making it sufficient to prove the good applicability of the model to other river areas where sand mines are operating. The study results will help policymakers in the Tien river area (An Giang, Dong Thap) in planning bank protection works, agricultural planning, and sand mining management.
In addition, harnessing the power of accelerators such as GPUs can help run the numerical simulations up to 23x times faster. This helps quickly verify the practical data collected coarsely and reduces the risk of critical situations, such as resident evacuations, when riverbank failure happens. In future work, we plan to extend our implementation to a distributed GPU system to study larger areas of rivers (or even oceans) using the numerical model.  Chezy coefficient, Ch = 18 log 12 H K s K s Bed roughness, K s = 3D 90 D 90 Diameter of particle that is equal to or less than 90% of the mass of the particles present S Standing for erosion and deposition rates (m/s). D Mean diameter of the particle (m)