Lattice Boltzmann Method-based Numerical Simulation for Heavy Metal Migration Process during Deep-sea Mining

1. Introduction

Due to the scarcity of land resources and the augmented economic growth demand for mineral resources, the deep-sea mining has become a hot research topic. However, the mineral resources are mainly distributed in the submarine sediments located thousands of meters deep in the ocean. The complexity and variability of the marine environment lead to the higher technical requirements for the detection and mining of mineral resources as well as environmental assessment [1,2,3,4]. During deep-sea mining, the subsidence, migration and transformation of submarine sediments can cause the migration of heavy metals, and the heavy metals are dissolved in water, which can lead to the marine pollution. Therefore, it is of momentous practical significance to model and simulate the migration and transformation of heavy metals for the marine environment evaluation.

In the last years, the migration and transformation of heavy metals pollutants in natural water and laboratory have been investigated. Generally, the heavy metals pollutants are enriched in sediment particles, which are thought to be the most essential carrier of heavy metals pollutants following water. Taking into account the effects of convection and diffusion on the concentration of heavy metals pollutants particle phase, Huang et al. [5] established a mathematical model for the transport and transformation of heavy metals pollutants in fluvial rivers. In addition, the migration and transformation of heavy metals pollutants in a steady, uniform equilibrium sediment-laden flow was also simulated in [6]. Considering the adsorption process of heavy metals on sediment, He et al. [7] constructed an equilibrium model and a non-equilibrium model, respectively. By means of the hydrodynamic model and sediment transport model, Periáñez [8] built a model considering the interaction of water, sand and metal to study the diffusion laws of heavy metals. Moreover, Wang et al. [9] developed a coupled model based on environmental fluid dynamics code and applied it to Taihu Lake, and the results show that the transformation trend of the heavy metals pollutants is similar with that of sediment ones. Horvat et al. [10] devolved a two-dimensional numerical model suiting approach to simulate the complex flow, sediment transport and heavy metals transport conditions in natural watercourses. By assessing the performance of the distributed hydrological model for simulating the transport of various heavy metals in rivers, Bouragba et al. [11] utilized a hydrological model to numerically calculate the migration of multiple heavy metals, i.e. Pb, Hg, Cr, and Zn, in the Harrach Rivera. Zeng et al. [12] studied the dissolution of heavy metals in the mining transportation process, and investigated its impact on seawater quality. Igoni et al. [13] estimated the levels of heavy metals Mg, Cd, and Ni using a one-dimensional transport model in dry and rainy seasons in the Marine-Base waterfront. However, the dissolution and transportation of heavy metals pollutants are very complex due to the deep-sea environment. Therefore, it is necessary to establish an accurate model to describe the dissolution and transportation of heavy metals pollutants.

Lattice Boltzmann method (LBM) was proposed in the mid-1980s. Based on the kinetic theory, LBM has the clear physical background and distinct features, including program simplicity, easy implementation, and natural parallelism [15]. In contrast with traditional methods, such as finite difference method, finite element method and spectral method, LBM is a mesoscopic numerical simulation method. It regards the motion of all particles as a whole, and their motion characteristics are expressed by distribution functions. So, LBM is used to study the mechanism of various complex phenomena, and it has been widely used in the engineering simulation, fluid mechanics, material engineering and environmental engineering [16,17,18,19]. LBM particularly is extensively concerned by scholars who research the mathematical and physical equations as convection-diffusion equations, diffusion equations, fractional Cahn-Hilliard equations [20,21,22,23]. Therefore, LBM can be adopted to solve the model for describing the dissolution and transportation of heavy metals pollutants during deep-sea mining.

In this work, the coupling model is established to describe the dissolution and transportation of heavy metals pollutants during deep-sea mining, and the LBM is adopted to solve the model simulation numerically. The paper is arranged as follows: the transfer-transformation model of heavy metals pollutants during deep-sea mining is constructed in Section 2. In Section 3, we will present the Lattice Boltzmann model to recover macroscopic equations by Chapman-Enskog analysis. Then Section 4 focuses on the numerical simulation. Finally, the conclusions are given in Section 5.

2. Transfer-Transformation Model of heavy metals Pollutants in Deep-Sea Mining

During deep-sea mining, heavy metals exist in seawater, suspended sediment and seabed sediment in a certain proportion. And the migration of heavy metals depends on the motion laws of the dissolved phase, suspended phase and sediment phase. The dissolution and transportation of heavy metals are affected by the convection-diffusion process, adsorption-desorption process and sediment settlement and resuspension.

2.1. Convection-diffusion model of heavy metals in water

According to the basic principle of water environment, the convective migration flux can be express as

$$f_x=u_x{C}_w, f_y=u_y{C}_w,$$

(1)

where f_x and f_y are the convective migration flux of heavy metals in x and y directions, respectively. u_x and u_y are the velocity components of water in x and y directions, and C_w stands for the concentration of heavy metals in water.

Based on the Fick 's first law, the mass flux of diffusion is proportional to the gradient of concentration,

$$I_x=-E_x\frac{\partial C_w}{\partial x},I_y=-E_y\frac{\partial C_w}{\partial y},$$

(2)

where I_x and I_y represent the diffusive flux of heavy metals in x and y directions, respectively. E_x and E_y are the diffusion coefficients. Based on the conservation law of mass, the convection-diffusion equation is

$$\frac{\partial C}{\partial t}+v\Delta C=E\Delta C.$$

(3)

2.2. Model of adsorption-desorption process

Most heavy metals are adsorbed by suspended sediment, and the desorption process is an inverse process of adsorption. The final trend of adsorption-desorption process is the decrease of heavy metals concentrations. The quantity of adsorption-desorption isotherm and kinetic models have been researched [24,25,26]. The typical Langmuir kinetic model is adopted in this paper,

$$f\left(N,C_w,t\right)=k_1{C}_w\left(b-N\right)-k_{2}N,$$

(4)

where N is the adsorption content of heavy metals of unit weight of sediment. k₁ represents the coefficient of adsorption rate, and k₂ represents the coefficient of desorption rate. b is the maximum absorption capacity of sediment.

2.3. Model of sediment settlement and resuspension

It is a natural purification process of natural water that the suspended particles suspended in water are deposited to form sediment due to gravity. However, due to the change of water environment, sediment may be suspended in water again, which causes secondary pollution. Therefore, it is of great significance to study the settlement and resuspension of suspended solids in environmental protection research. According to the principle of sediment engineering, the sedimentation process of heavy metals in suspension solids and the resuspension of sediment are described as follows [27],

$$(K-1)w_0\frac{\partial C_s}{\partial z}=\frac{K{u}^m{C}_d-K_w{u}^{-n}{C}_s}{H_i}$$

(5)

where H_i is the average depth, and u denotes the velocity of water. K_w represents the comprehensive influence constant. w₀ means velocity of sedimentation, and K stands for the resuspension coefficient. m=4, and n=2.

The suspended sediment is the carrier for the transport and transformation of heavy metals pollutants. The coupling model for transport and transformation of heavy metals is established by describing the transport of heavy metals with dissolved phase, suspended phase and sediment phase. For the dissolved and suspended phases affected by the migration of water, the adsorption and desorption of heavy metals are determined by suspended sediment and seabed sediment. For the sediment phase, the diffusion effect of heavy metals can be ignored, and the adsorption and desorption are considered. Based on [6], the improved coupling model for the migration process of heavy metals during deep-sea mining can be expressed as,

$$\left\{\begin{array}{l}\frac{\partial C_w}{\partial t}+u\nabla C_w=v\Delta C_w-s_s\left(k_1^w{C}_w\left(b^w-C_s\right)-k_2^w{C}_s\right)-\frac{1-p}{H_i}\left(k_1^b{C}_w\left(b^b-C_d\right)-k_2^b{C}_d\right),\hfill \\ \frac{\partial C_s}{\partial t}+u\nabla C_s=v\Delta C_s+k_1^w{C}_w\left(b^w-C_s\right)-k_2^w{C}_s+\frac{K{u}^m{C}_d-K_w{u}^{-n}{C}_s}{H_i},\hfill \\ \frac{\partial C_d}{\partial t}=k_1^b{C}_w\left(b^b-C_d\right)-k_2^b{C}_d-\frac{K{u}^m{C}_d-K_w{u}^{-n}{C}_s}{H_i},\hfill \end{array}\right.$$

(6)

where C_w, C_s, C_d represent the concentrations of heavy metals in dissolved phase, suspended phase and sediment phase, respectively. u=(u_x, u_y, u_z) is convection coefficient, and v=(v_x, v_y, v_z) is diffusion coefficient. The suspended sediment concentration is denoted as s_s. p is the porosity. k₁^w and k₂^w are expressed as adsorption rate coefficient and desorption rate coefficient of suspended solids. k₁^b and k₂^b represent adsorption rate coefficient and desorption rate coefficient of sediment. b^w and b^b are the saturation adsorption in suspended solids and sediment, respectively.

3. Lattice Boltzmann model for migration-transformation of heavy metals

In this section, the LBM is adapted to recover the macroscopic equations which are required in Eq. (6).

3.1. Lattice Boltzmann model

Applying the LBGK collision operator, the evolution equation with a source term can be written as

$$f_{ki}\left(X+c_i\Delta t,t+\Delta t\right)-f_{ki}(X,t)=-\frac{1}{{\tau }_k}\left[f_{ki}(X,t)-f_{ki}^{eq}(X,t)\right]+{\epsilon }^2{F}_{ki}(X,t),$$

(7)

where $f_{ki}(X,t)$ is the distribution functions in the $i^{th}$ direction of the $k^{th}$ equation. $f_{ki}^{eq}(X,t)$ is the equilibrium distribution function. ${\tau }_k$ denotes the dimensionless relaxation time. And $f_{ki}^{eq}(X,t)$ satisfies the following conditions,

$$\begin{array}{l}{\displaystyle \sum _i{f}_{1i}}={\displaystyle \sum _i{f}_{1i}^{(eq)}}=C_w,{\displaystyle \sum _i{f}_{2i}}={\displaystyle \sum _i{f}_{2i}^{(eq)}}=C_s,{\displaystyle \sum _i{f}_{3i}}={\displaystyle \sum _i{f}_{3i}^{(eq)}}=C_d,\hfill \\ {\displaystyle \sum _i{c}_i}{f}_{1i}=u_1{C}_w,{\displaystyle \sum _i{c}_i}{f}_{2i}=u_2{C}_s,{\displaystyle \sum _i{c}_i}{f}_{3i}=0,\hfill \\ {\displaystyle \sum _i{c}_i}{c}_i{f}_{1i}={\lambda }_1{C}_w+u_1^2{C}_w,{\displaystyle \sum _i{c}_i}{c}_i{f}_{2i}={\lambda }_2{C}_s+u_2^2{C}_s,{\displaystyle \sum _i{c}_i}{c}_i{f}_{3i}=0,\hfill \end{array}$$

(8)

where ${\lambda }_k=v_k/\Delta t({\tau }_k-\frac{1}{2}),k=1,2,3$.

And the source term satisfies the following conditions,

$$\begin{array}{l}{\displaystyle \sum _i{F}_{1i}}=F_{1i}(C_w,C_s,C_d),{\displaystyle \sum _i{F}_{2i}}=F_2(C_w,C_s,C_d),{\displaystyle \sum _i{F}_{3i}}=F_3(C_w,C_s,C_d).\hfill \end{array}$$

(9)

3.2. The chapman-Enskog analysis

Applying the Taylor expansion to Eq. (7), the evolution equation can be expressed as,

$$f_{ki}\left(X+c_i\Delta t,t+\Delta t\right)-f_i(X,t)=\left(c_i\frac{\partial f_{ki}}{\partial X}+\frac{\partial f_{ki}}{\partial t}\right)\Delta t+{\left(c_i\frac{\partial f_{ki}}{\partial X}+\frac{\partial f_{ki}}{\partial t}\right)}^2\frac{{(\Delta t)}^2}{2!}+O\left(\Delta t^3\right).$$

(10)

Through the Chapman-Enskog analysis, the derivatives of time and space, and distribution function are expanded as,

$$\frac{\partial }{\partial t}=\epsilon \frac{\partial }{\partial t_1}+{\epsilon }^2\frac{\partial }{\partial t_2},\frac{\partial }{\partial X}=\epsilon \frac{\partial }{\partial X_1},f_{ki}=f_{ki}^{(0)}+\epsilon f_{ki}^{(1)}+{\epsilon }^2{f}_{ki}^{(2)}+O\left({\epsilon }^3\right),$$

(11)

denoting $D_{1i}=\frac{\partial }{\partial t1}+c_i\frac{\partial }{\partial X_1}$ and substituting (11) into (10), one can obtain

$$\begin{array}{l}{f}_{ki}\left(X+c_i\Delta t,t+\Delta t\right)-f_{ki}(X,t)\\ =\Delta t\left(\epsilon D_{ki1}+{\epsilon }^2\frac{\partial }{\partial t_2}\right)\left(f_{ki}^{(0)}+\epsilon f_{ki}^{(1)}+{\epsilon }^2{f}_{ki}^{(2)}+O\left({\epsilon }^3\right)\right)\\ +\frac{\Delta t^2}{2!}{\left(\epsilon D_{ki1}+{\epsilon }^2\frac{\partial }{\partial t_2}\right)}^2\left(f_{ki}^{(0)}+\epsilon f_{ki}^{(1)}+{\epsilon }^2{f}_{ki}^{(2)}+O\left({\epsilon }^3\right)\right)+O\left({\epsilon }^3\right)\\ =-\frac{1}{{\tau }_k}\left[\epsilon f_{ki}^{(1)}+{\epsilon }^2{f}_{ki}^{(2)}+O\left({\epsilon }^3\right)\right]+{\epsilon }^2{F}_{ki}(X,t).\end{array}$$

(12)

On the basis of (12), we can derive the following equation in order of ${\epsilon }^0,{\epsilon }^1,{\epsilon }^2$,

$$\left\{\begin{array}{l}\frac{\partial }{\partial t}{\displaystyle \sum _i{f}_{1i}^{(0)}}+\frac{\partial }{\partial X}{\displaystyle \sum _i{c}_i}{f}_{1i}^{(0)}+\Delta t\left(\frac{1}{2}-\tau \right)\left(\frac{{\partial }^2}{\partial X^2}{\displaystyle \sum _i{c}_i}{c}_i{f}_{1i}^{(0)}+\frac{{\partial }^2}{\partial t\partial X}{\displaystyle \sum _i{c}_i}{f}_{1i}^{(0)}\right)=\frac{{\epsilon }^2}{\Delta t}{\displaystyle \sum _i{F}_{1i}}(X,t),\hfill \\ \frac{\partial }{\partial t}{\displaystyle \sum _i{f}_{2i}^{(0)}}+\frac{\partial }{\partial X}{\displaystyle \sum _i{c}_i}{f}_{2i}^{(0)}+\Delta t\left(\frac{1}{2}-\tau \right)\left(\frac{{\partial }^2}{\partial X^2}{\displaystyle \sum _i{c}_i}{c}_i{f}_{2i}^{(0)}+\frac{{\partial }^2}{\partial t\partial X}{\displaystyle \sum _i{c}_i}{f}_{2i}^{(0)}\right)=\frac{{\epsilon }^2}{\Delta t}{\displaystyle \sum _i{F}_{2i}}(X,t),\hfill \\ \frac{\partial }{\partial t}{\displaystyle \sum _i{f}_{3i}^{(0)}}=\frac{{\epsilon }^2}{\Delta t}{\displaystyle \sum _i{F}_{3i}}(X,t).\hfill \end{array}\right.$$

(13)

Combined with the conditions satisfied by the equilibrium distribution function, the macroscopic equations are accurately recovered. Additionally, the finite difference method (FDM) is applied to be compared with LBM. $\Delta x$ and $\Delta t$are the space sizes and time sizes. Denote $x_i=i\Delta x$, $t_j=j\Delta t$ and the difference scheme of Ref. [6] can be derived as

$$\left\{\begin{array}{l}{C}_{w,i}^{j+1}=a_1{C}_{w,i-1}^j+b_1{C}_{w,i}^j+c_1{C}_{w,i+1}^j+\Delta t{s}_s\left(k_1^w{C}_{w,i}^j-k_2^w\right)C_{s,i}^j+\Delta t\frac{1-p}{H_i}\left(k_1^b{C}_{w,i}^j-k_2^b\right)C_{d,i}^j,\hfill \\ C_{s,i}^{j+1}=a_2{C}_{s,i-1}^j+b_2{C}_{s,i}^j+c_2{C}_{s,i+1}^j+\Delta t{k}_1^w{C}_{w,i}^j\left(b^w-C_{s,i}^j\right)+\Delta t\frac{K{u}^m{C}_{d,i}^j}{H_i},\hfill \\ C_{d,i}^{j+1}=b_3{C}_{d,i}^j+\Delta t{k}_1^b{C}_{w,i}^j\left(b^b-C_{d,i}^j\right)+\frac{\Delta t{K}_w{u}^{-n}{C}_{s,i}^j}{H_i},\hfill \end{array}\right.$$

(14)

where

$$\begin{array}{l}{a}_1=\frac{u_1\Delta X\Delta t+v_1\Delta t}{\Delta X^2},b_1=1-u_1\frac{\Delta t}{\Delta X}-2v_1\frac{\Delta t}{\Delta X^2}-\Delta t{s}_s{k}_1^w{b}^w-\frac{\Delta t(1-p)k_1^b{b}^b}{H_i},\\ a_2=\frac{u_2\Delta X\Delta t+v_2\Delta t}{\Delta X^2},b_1=1-u_2\frac{\Delta t}{\Delta X}-2v_2\frac{\Delta t}{\Delta X^2}-\Delta t{k}_2^w{b}^w-\frac{\Delta t{K}_w{u}^{-n}}{H_i},\\ b_3=1-\Delta t{k}_2^b-\frac{\Delta tK{u}^m}{H_i},c_1=\frac{v_1\Delta t}{\Delta X^2},c_2=\frac{v_2\Delta t}{\Delta X^2},a_3=c_3=0.\end{array}$$

4. Numerical simulations for the migration-transformation of heavy metals in deep-sea mining

In this section, the numerical simulations for the migration-transformation of heavy metals pollutants in the dissolved phase, suspended phase and sediment phase are presented, and the proposed method is compared with the FDM. In addition, the concentration variation of heavy metals as time and distances is analyzed.

4.1. Spatial-temporal distribution of heavy metals concentrations

In this research, the Cadmium (Cd) is selected as an example, and some physical parameters are shown in Table [6].

Table 1. The physical parameters.

Physical quantities	parameters	Physical quantities	parameters
ss	1.8378 kg/m3	k1w	0.0076 L/(mg·s)
k1b	7.6e-4 L/(mg·s)	k2w	0.00084 L/s
k2b	8.4e-5 L/s	bb	5.34 g/m2
K	1.1e-6	bw	0.534 mg/g
Kw	9.0e-5	p	0.5

Table 1. The physical parameters.

Physical quantities	parameters	Physical quantities	parameters
ss	1.8378 kg/m3	k1w	0.0076 L/(mg·s)
k1b	7.6e-4 L/(mg·s)	k2w	0.00084 L/s
k2b	8.4e-5 L/s	bb	5.34 g/m2
K	1.1e-6	bw	0.534 mg/g
Kw	9.0e-5	p	0.5

In the simulations, $N=100, \Delta x=L/N, \Delta t=2$ and $c=\Delta x/\Delta t.$ The relaxation time $\tau =\mathrm{1.38.}$The boundary concentration $c_0=1mg/L,$which indicates the concentration of continuous discharge of polluted water flow. The diffusion and convection coefficients are $[v_1 v_2 v_3]=[0.29 0.29 0]$ and $[u_1 u_2 u_3]=[1.04 1.04 0]$, respectively.

Based on the established coupling model and LBM, the simulated concentration distribution of Cadmium in the dissolved phase, suspended phase and sediment phase are shown in Figure 1, Figure 2 and Figure 3. As shown in Figure 1 (b), the concentration of Cadmium is decreased along x axis when the time is fixed. The pollution range is increased with time. The concentration of Cadmium rapidly increases to a certain value, and then gradually tend to balance value.

The spatial-temporal concentration distribution of Cadmium in the suspension phase is shown in Figure 2 (a). As shown in Figure 2(b), the concentration of Cadmium is firstly increased and then decreased when the time is fixed. At the same location, the change rule of concentration in the suspension is similar with that in the dissolved phase. However, the position where the concentration reaches the highest point is gradually moving forward due to the spread of polluted flow, and the concentration is also constantly augmented.

The spatial-temporal concentration distribution of Cadmium in sediment phase is shown in Figure 3. Obviously, the change rule of concentration is simple. The concentration of Cadmium is decreased along x axis when the time is fixed, but the concentration of Cadmium is increased at the same location with more time for the continuous spread of polluted flow.

4.2. Compare with finite difference schemes

The concentration variations of Cadmium in different phases simulated by LBM are compared with those simulated by FDM. From Figure 4, the concentration variations of Cadmium by LBM and FDM are nearly consistent. In addition, three location spots of the dissolved phase (66 m), suspended phase (72 m) and sediment phase (24 m) are selected to compare the concentrations of Cadmium between LBM and FDM, and the errors are very small. However, the time consumption of LBM is less than that of FDM. This is because that the time step $(\Delta t=0.2)$ of FDM is smaller to satisfy the stability conditions.

4.3. Analysis for the concentration variation of heavy metals

The concentration variations of Cadmium vs. distance are shown in Figure 5. Obviously, the concentration of Cadmium is decreased with increases of distance in the dissolved phase. However, the decrease trend becomes slower near the peak, which is caused by rapidly moving and slower settling of dissolved phase. The concentration of Cadmium initially increases and then gradually decreases in the suspended phase, and the concentration of Cadmium near the peak is always at the highest point. This is because that the suspended sediment moves along with the flow, and the heavy metals have enough time to be absorbed by the suspended sediment. In the sediment phase, the sediment particles nearly do not move with the flow, so the concentration of sediment deposited at the entrance is the highest, and the concentration decreases with the increase of distance from the upstream.

The concentration variations of Cadmium vs. time is shown in Figure 6, and the observation location is 20 m away from upstream. The concentration of Cadmium is 0 mg/L, when the polluted water has not spread to this position. The concentration of Cadmium is increased as the polluted flow spreads over time. The concentration of Cadmium in the dissolved phase is gradually increased and then reaches equilibrium, which is the similar with that in the suspension phase, but the concentration of Cadmium in the suspension phase is slowly increased. This is because that the polluted water flows continuously, and the heavy metals absorbed in the suspended phase migrates follow polluted water. In addition, the concentration of Cadmium in the sediment phase is increased, because of the low migration rate of sediment and the little influence of convective diffusion.

The effect of heavy metals on the water quality is simply investigated by the proposed model. Clearly, the pollution range in the dissolved phase is expanded with time. When the concentration of Cadmium at the entrance is 0.02 mg/L, the concentration of Cadmium in the dissolved phase is less than 0.01 mg/L at the location beyond 51m at 50s, as shown in Figure 7 (a). The concentration of Cadmium in the suspension phase is not exceeding 0.01mg/L within 150 s, as shown in Figure 7(b). Although the concentration of Cadmium is more than 0.01 mg/L in the sediment phase, there is almost no diffusion of Cadmium in sediment. So the concentration of Cadmium at the location beyond 93 m in the sea can satisfy the water quality standard. Therefore, the pollution of heavy metals should be controlled according to the mining time and the initial concentration of heavy metals, which can guarantee the water quality to meet the standards within a certain range during deep-sea mining.

5. Conclusion

In this paper, the convection-diffusion, adsorption-desorption, and the sedimentation-resuspension of heavy metals are analyzed during deep-sea mining, and the coupling model is established to simulate the migration of heavy metals in the dissolved phase, suspended phase and sediment phase. The LBM is used to numerically solve the coupling model. Compared with FDM, the LBM takes less time, and the errors between LBM and FDM are very small. The spatial/temporal distribution and movement regularity of Cadmium are investigated by the developed model. The concentration variations of Cadmium with time and distance are discussed and analyzed in detail, which can provide the reasonable bases for the evaluation of Cadmium pollution. In addition, the developed model can be used to describe the migration of other heavy metals pollutants during deep-sea mining.