# Realistic Simulation of Dissolution Process on Rock Surface
## 4 Dissolution Simulation
### Calculation of Flow Trajectory  
In the process of dissolution, the movement track of raindrops will be affected by many factors, including the current movement direction, flow velocity, terrain slope, and fracture sequence in rocks. In this paper, these factors are comprehensively considered, and the motion trajectory of raindrops is formed by continuously calculating the position of droplets at the next moment.  
The following attributes are stored in each droplet: position \(P\), flow velocity \(V\), water volume \(W\), carbonate ion content \(M\), and carbon dioxide content \(C\). Water is reduced by evaporation, and the reduction depends on the evaporation rate. When the water volume \(W=0\), the droplet will ceases to exist. When \(W\) is not zero, droplets either flow along the rock fracture or slide along the slope in the direction of the maximum gradient.  
Assume the current position of a droplet is \(P_{k}\). When the droplet is just in a rock fracture, its movement path is completely affected by the fracture trend. According to its current movement direction \(V_{old}\), we select the point closest to the direction vector in the fracture sequence as the alternative point \(P_{k^{\prime}}\). Considering that the water droplet cannot flow to a position with an altitude higher than the current position, we can determine the position \(P_{k+1}\) of the droplet at the next moment by comparing the current height and the height of the alternative position, as follows:  
\[P_{k+1}=Max\_height(P_{k},P_{k^{\prime}}) \tag{1}\]  
If the height of the current position \(P_{k}\) is higher than the height of the alternative position \(P_{k^{\prime}}\), it means that the droplet will continue to stay at the current position for dissolution.  
When the droplet is not in the rock fracture, we calculate gradient \(G\) at \(P_{k}\), thus obtaining the next position \(P_{k+1}\). Considering that the position of the water droplet may not necessarily fall exactly on the grid point, there might be a deviation. We use the distance difference (\(u\),\(v\)) between its current position and the integer coordinate point to perform bilinear interpolation on the three grid points adjacent to \(P_{k}\), so as to obtain the gradient \(G\) of the current position, as follows:  
\[G(P_{k})=\begin{pmatrix}(P_{x+1,z}-P_{x,z})\times(1-v)+(P_{x+1,z+1}-P_{x,z+1}) \times(v)\\ (P_{x,z+1}-P_{x,z})\times(1-u)+(P_{x+1,z+1}-P_{x+1,z})\times(u)\end{pmatrix} \tag{2}\]  
where \(u\) and \(v\) are the offset of the position of the water droplet with respect to the coordinate value of the grid point, and \(P_{x,z}\) is the position of the grid point closest to \(P_{k}\).  
Considering that the water droplet has inertia, its motion is not only affected by the slope of terrain but also by its current direction. Therefore, we can obtain the new direction \(V_{new}\) via linear interpolation of the current direction \(V_{old}\) and gradient \(G\).  
\[V_{new}=V_{old}\times Inertia-G\times(1- Inertia) \tag{3}\]  
where \(Inertia\)\(\in\)(0, 1) is the user-defined coefficient. Inertia = 1 means that the raindrop will always move in the initial direction, without considering the influence of the gradient. Inertia = 0 means that the direction of the raindrop is only affected by the gradient. The next position \(P_{k+1}\) can be obtained by adding the flow direction \(V_{new}\) and the current position \(P_{k}\). The directions mentioned above are normalized. Since the flow direction is a normalized value, the distance that the raindrop moves across will not exceed one unit, so the position of the raindrop at the next moment will only be its adjacent point.