Lattice-Boltzmann-Method-Based Numerical Simulation for Heavy Metal Migration Process during Deep-Sea Mining

Abstract

During deep-sea mining, heavy metal pollutants can cause contamination in the marine environment. In this paper, the multiphasic coupling model is established to describe the heavy metal migration process during deep-sea mining, which takes the effects of the convection–diffusion, adsorption–desorption, and sedimentation–resuspension of heavy metals in the aquatic environment into full consideration. Due to the advantages of the Lattice Boltzmann method, it is adopted to numerically solve the multiphasic coupling model and achieve the simulation of the heavy metal migration process during deep-sea mining. In addition, taking cadmium as an example, the concentration variations are discussed and analyzed in detail. Based on the established model and Lattice Boltzmann method, the concentration distribution of heavy metals can be accurately described to provide the reasonable bases for the evaluation of marine environmental protection.

1. Introduction

Due to the scarcity of land resources and the augmented economic growth demand for mineral resources, deep-sea mining has become an important research topic. However, the mineral resources are mainly distributed in the submarine sediments located thousands of meters deep in the ocean. The complexity of and variability in the marine environment lead to 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 marine pollution. Therefore, it is of practical significance to model and simulate the migration and transformation of heavy metals for marine environment evaluation.

In recent years, the migration and transformation of heavy metal pollutants in natural water and the laboratory have been investigated. Generally, the heavy metal pollutants are enriched in sediment particles, which are thought to be the essential carrier of heavy metal pollutants following water. Taking into account the effects of convection and diffusion on the concentration of the heavy metal pollutants particle phase, Huang et al. [5] established a mathematical model for the transport and transformation of heavy metal pollutants in fluvial rivers. In addition, the migration and transformation of heavy metal pollutants in a steady, uniform equilibrium sediment-laden flow were 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 metal pollutants is similar to that of sediment ones. Geng et al. [10] also studied the release of heavy metals with suspended sediment in Taihu Lake under hydrodynamic conditions. Based on the Lattice Boltzmann method (LBM), a numerical model of heavy metals released into overlying water under hydrodynamic conditions is established. Horvat et al. [11] developed a two-dimensional-numerical-model-suiting approach to simulate the complex flow, sediment transport, and heavy metal 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. [12] 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. [13] studied the dissolution of heavy metals in the mining transportation process and investigated its impact on seawater quality. Igoni et al. [14] estimated the levels of heavy metals Mg, Cd, and Ni using a one-dimensional transport model in dry and rainy seasons in a marine base waterfront. In Ref. [15], He analyzed the migration of heavy metal contaminants in a tailings pond and the retention behavior of a compacted bentonite engineered barrier system on the future of heavy metal contaminants. Based on the convection–dispersion equation, He et al. [16] also established a model of landfill leachate to investigate the retention behavior of a laterite–bentonite engineered barrier towards heavy metal contaminants. The concentration of heavy metals (cadmium, chromium, cobalt, copper, iron, and lead) in the water and sediment of Lake Greenwood, and in the streams and creeks that flow into the lake, was reported in [17].

However, the process of the dissolution and transport of heavy metal pollutants is very complex due to the deep-sea environment. Therefore, based on previous studies, this paper further considers the influence of multiple factors to establish the model to describe the dissolution and transport of heavy metal pollutants more accurately, and uses a new numerical method to solve the problem.

The Lattice Boltzmann method (LBM) was proposed in the mid-1980s. Based on kinetics theory, the LBM has a clear physical background and distinct features, including program simplicity, easy implementation, and natural parallelism [18,19,20]. In contrast with traditional methods, such as the finite difference method, finite element method, and spectral method, the 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, the LBM is used to study the mechanism of various complex phenomena, and it has been widely used in fluid mechanics, engineering simulation, material engineering, and environmental engineering [21,22,23,24,25]. The LBM is extensively of interest by scholars who research mathematical and physical equations, such as convection–diffusion equations, diffusion equations, and fractional Cahn–Hilliard equations [26,27,28,29,30]. Therefore, the LBM can be adopted to solve the model for describing the dissolution and transportation of heavy metal pollutants during deep-sea mining.

In this work, the coupling model is established to describe the dissolution and transportation of heavy metal pollutants during deep-sea mining, and the LBM is adopted to solve the model numerically. This paper is arranged as follows: the transfer-transformation model of heavy metal pollutants during deep-sea mining is constructed in Section 2. In Section 3, we 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 Metal 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 laws of motion 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, as shown in Figure 1. Precisely, the physical movements of heavy metals transported in the dissolved phase represent convective and diffusive mass transfer. The suspended phase is mainly characterized by the adsorption of heavy metals from the dissolved phase by the suspended matter, as well as desorption under certain conditions. In the sediment phase, there is the adsorption of heavy metals, as well as the sedimentation and resuspension transfer of suspended matter.

2.1. Convection–Diffusion Model of Heavy Metals in Water

According to the basic principle of the 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 the x and y directions, respectively. u_x and u_y are the velocity components of water in the x and y directions, and C_w stands for the concentration of heavy metals in water.

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

$$I_x=-v_x\frac{\partial C_w}{\partial x},I_y=-v_y\frac{\partial C_w}{\partial y},$$

(2)

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

$$\frac{\partial C_w}{\partial t}+u_x\frac{\partial C_w}{\partial x}+u_y\frac{\partial C_w}{\partial y}=v_x\frac{{\partial }^2{C}_w}{\partial x^2}+v_y\frac{{\partial }^2{C}_w}{\partial y^2}.$$

(3)

where u_x and u_y represent the velocity components of water in the x and y directions. V_x and v_y represent the diffusion coefficients in the x and y directions, respectively.

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 the adsorption–desorption process is the decrease in heavy metal concentrations. A number of adsorption–desorption isotherms and kinetic models have been researched [31,32,33]. 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 a 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 water that the 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 [34]:

$$(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 the 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 metal pollutants. The coupling model for the transport and transformation of heavy metals is established by describing the transport of heavy metals with a 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_{1x}\frac{\partial C_w}{\partial x}+u_{1y}\frac{\partial C_w}{\partial y}=v_{1x}\frac{{\partial }^2{C}_w}{\partial x^2}+v_{1y}\frac{{\partial }^2{C}_w}{\partial y^2}-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_{2x}\frac{\partial C_s}{\partial x}+u_{2y}\frac{\partial C_s}{\partial y}=v_{2x}\frac{{\partial }^2{C}_s}{\partial x^2}+v_{2y}\frac{{\partial }^2{C}_s}{\partial y^2}+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, and C_d represent the concentrations of heavy metals in the dissolved phase, suspended phase, and sediment phase, respectively. u_{i =} (u_ix,u_iy)(i = 1,2) are convection coefficients, and v_{i =} (v_ix,v_iy)(i = 1,2) are diffusion coefficients. The suspended sediment concentration is denoted as s_s. p is the porosity. k₁^w and k₂^w are expressed as the adsorption rate coefficient and desorption rate coefficient of suspended solids. k₁^b and k₂^b represent the adsorption rate coefficient and desorption rate coefficient of sediment. b^w and b^b are the saturation adsorption parameters in suspended solids and sediment, respectively.

3. Lattice Boltzmann Model for Migration–Transformation of Heavy Metals

The LBM is a computational model that discretizes time, space, and the velocity space of particles. Its evolution process takes the unit time step as an independent evolution step, and two processes are included in each time step:

(1)

Collision. The particles collide at the lattice point, generating velocity to the surrounding area.

(2)

Flow. After the collision, the particles flow in accordance with the direction of the discretized velocity, which is embodied in the movement of particles at each lattice point to the surrounding neighboring lattice points.

In this section, the LBM is adapted to recover the macroscopic equations which are required in Equation (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+{\mathit{c}}_{\mathit{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 function 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. ${\epsilon }^2{F}_{ki}(X,t)$ represents the correction term added to recover the source term, and the parameter $\epsilon$ has quantities of the same order as the Knudsen in multiscale expansion. 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{\mathit{c}}_{\mathit{i}}}{f}_{1i}^{(eq)}={\mathit{u}}_1{C}_w,{\displaystyle \sum _i{\mathit{c}}_{\mathit{i}}}{f}_{2i}^{(eq)}={\mathit{u}}_2{C}_s,{\displaystyle \sum _i{\mathit{c}}_{\mathit{i}}}{f}_{3i}^{(eq)}=0,\hfill \\ {\displaystyle \sum _i{\mathit{c}}_{\mathit{i}}}{\mathit{c}}_{\mathit{i}}{f}_{1i}^{(eq)}={\Lambda }_1{C}_w+{\mathit{u}}_1^2{C}_w,{\displaystyle \sum _i{\mathit{c}}_{\mathit{i}}}{\mathit{c}}_{\mathit{i}}{f}_{2i}^{(eq)}={\Lambda }_2{C}_s+{\mathit{u}}_2^2{C}_s,{\displaystyle \sum _i{\mathit{c}}_{\mathit{i}}}{\mathit{c}}_{\mathit{i}}{f}_{3i}^{(eq)}=0,\hfill \end{array}$$

(8)

where ${\mathbf{\lambda }}_k={\mathit{v}}_k/\Delta t({\tau }_k-\frac{1}{2}), (k=1,2),{\Lambda }_k=diag({\lambda }_{kx},{\lambda }_{ky})$.

And the source term satisfies the following conditions:

$$\begin{array}{l}{\displaystyle \sum _i{F}_{1i}}=F_1(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 Equation (7), the evolution equation can be expressed as

$$f_{ki}\left(X+{\mathit{c}}_i\Delta t,t+\Delta t\right)-f_i(X,t)=\left({\mathit{c}}_i\frac{\partial f_{ki}}{\partial X}+\frac{\partial f_{ki}}{\partial t}\right)\Delta t+{\left({\mathit{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 Chapman–Enskog analysis [35], the derivatives of time and space and the distribution function can be 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 t_1}+{\mathit{c}}_i\frac{\partial }{\partial X_1}$ and substituting (11) into (10), one can obtain

$$\begin{array}{l}{f}_{ki}\left(X+{\mathit{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{\mathit{c}}_i}{f}_{1i}^{(0)}+\Delta t\left(\frac{1}{2}-\tau \right)\left(\frac{{\partial }^2}{\partial X^2}{\displaystyle \sum _i{\mathit{c}}_i}{\mathit{c}}_i{f}_{1i}^{(0)}+\frac{{\partial }^2}{\partial t\partial X}{\displaystyle \sum _i{\mathit{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{\mathit{c}}_i}{f}_{2i}^{(0)}+\Delta t\left(\frac{1}{2}-\tau \right)\left(\frac{{\partial }^2}{\partial X^2}{\displaystyle \sum _i{\mathit{c}}_i}{\mathit{c}}_i{f}_{2i}^{(0)}+\frac{{\partial }^2}{\partial t\partial X}{\displaystyle \sum _i{\mathit{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 one-dimensional macroscopic equations are accurately recovered:

$$\left\{\begin{array}{l}{f}_{11}^{(eq)}=C_w-\frac{{\lambda }_{1x}{C}_w+u_{1x}^2{C}_w}{c^2}, f_{12}^{(eq)}=-f_{13}^{(eq)}=\frac{c{u}_{1x}{C}_w+{\lambda }_{1x}{C}_w+u_{1x}^2{C}_w}{2c^2},\\ f_{21}^{(eq)}=C_s-\frac{{\lambda }_{2x}{C}_s+u_{2x}^2{C}_s}{c^2}, f_{22}^{(eq)}=-f_{23}^{(eq)}=\frac{c{u}_{2x}{C}_s+{\lambda }_{2x}{C}_2+u_{2x}^2{C}_s}{2c^2},\\ f_{31}^{(eq)}=C_d/3, f_{32}^{(eq)}=f_{33}^{(eq)}=0,\\ F_{1i}=-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),\\ F_{2i}=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},\\ F_{3i}=k_1^b{C}_w\left(b^b-C_d\right)-k_2^b{C}_d, i=1,2,3.\end{array}\right.$$

(14)

Additionally, the finite difference method (FDM) is applied to compare with the LBM. $\Delta x$ and $\Delta t$ are the space sizes and time sizes. $x_i=i\Delta x$ and $t_j=j\Delta t$, and the difference scheme of [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.$$

(15)

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 metal 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. Spatiotemporal Distribution of Heavy Metal Concentrations

The following assumptions are introduced to simplify the problem.

(1)

The biological enrichment and decay processes of heavy metals are not considered.

(2)

The diffusion of heavy metals in the sedimentary phases is small and negligible.

(3)

The exchange between flowing waters and the consequent exchange of heavy metals are negligible.

(4)

The specificity of the water body itself is not taken into account.

In this research, cadmium (Cd) is selected as an example, and some physical parameters are shown in

Table 1. The physical parameters.

Physical Quantities	Parameters	Physical Quantities	Parameters
ss	1.8378 kg/m3	k1w	7.6 × 10−3 L/(mg·s)
k1b	7.6 × 10−4 L/(mg·s)	k2w	8.4 × 10−4 L/s
k2b	8.4 × 10−5 L/s	bb	5.34 g/m2
K	1.1 × 10−6	bw	0.534 mg/g
Kw	9.0 × 10−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 =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 developed coupling model and LBM, the overall trend distribution of the simulated spatiotemporal variation in cadmium concentration in the dissolved, suspended, and sedimentary phases is shown in Figure 2a, Figure 3a and Figure 4a. The concentration of cadmium decreases along the x-axis when the time is fixed, as shown in Figure 2b, which shows the cross-sections of Figure 2a. The pollution range increases with time. The concentration of cadmium rapidly increases to a certain value, and then gradually tends to a balanced value.

The spatiotemporal concentration distribution of cadmium in the suspension phase is shown in Figure 3a. As shown in Figure 3b, the concentration of cadmium first increases and then decreases, when the time is fixed. At the same location, the change rule of concentration in the suspension is similar to that in the dissolved phase. However, the position where the concentration reaches the highest point gradually moves forward due to the spread of polluted flow, and the concentration is also constantly augmented.

The spatiotemporal concentration distribution of cadmium in the sediment phase is shown in Figure 4. Obviously, the change rule of concentration is simple. The concentration of cadmium decreases along the x-axis when the time is fixed, but the concentration of cadmium increases 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 the LBM are compared with those simulated by the FDM. From Figure 5, the concentration variations of cadmium by the 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 the LBM and FDM, and the errors are very small. However, the time consumption of the LBM is less than that of the FDM. This is because the time step $(\Delta t=0.2)$ of the 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 6. Obviously, the concentration of cadmium decreases with the increase in distance in the dissolved phase. However, the decrease trend becomes slower near the peak, which is caused by the rapidly moving and slower settling of the 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 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 highest, and the concentration decreases with the increase in distance from the upstream.

The concentration variations of cadmium vs. time are shown in Figure 7, and the observation location is 20 m away from the upstream. The concentration of cadmium is 0 mg/L when the polluted water has not spread to this position. The concentration of cadmium increases as the polluted flow spreads over time. The concentration of cadmium in the dissolved phase gradually increases and then reaches equilibrium, which is similar to that in the suspension phase, but the concentration of cadmium in the suspension phase slowly increases. This is because the polluted water flows continuously, and the heavy metals absorbed in the suspended phase migrate by following the polluted water. In addition, the concentration of cadmium in the sediment phase increases, 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 investigated by the proposed model. Clearly, the pollution range in the dissolved phase expands 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 51 m at 50 s, as shown in Figure 8a. The concentration of cadmium in the suspension phase does not exceed 0.01 mg/L within 150 s, as shown in Figure 8b. 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. Conclusions

In this paper, the convection–diffusion, adsorption–desorption, and 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 the FDM, the LBM takes less time, and the errors between the LBM and FDM are very small. The spatiotemporal 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 metal pollutants during deep-sea mining. Moreover, the absorption of heavy metals by marine flora and fauna will be systematically investigated and modeled in our future research.