### 1. Introduction

^{10}ton per year [1]. The engineering problem of the repository as a source of contaminant located below the water source or above the aquifer is one of the most important challenges. Natural weak diffusive clay layer or geosynthetic layer with slight diffusive capabilities is commonly used to cover underwater contaminant sources [2], such as the stratified contaminated soil layer containing a large amount of high-concentration and highly corrosive leachate. It should be pointed out that under the coupling effect of multiple fields, the contaminant source of the repository is likely to diffuse to the adjacent soil and water with great difficulty to repair again, which seriously threatens the environmental health [3–5]. In the above cases, molecular diffusion induced by concentration differences is the key factor for the contaminant diffusion to the adjacent soil or cover layer because of its low hydraulic conductivity of the clay (e.g., < 10

^{−7}cm/s), low leachate head (e.g., < 0.3 m), and small Peclet number [6–9]. Besides, diffusive transport problem comprehensively exists in the engineering materials, such as compacted clay liner [10, 11], underwater cap [12], geosynthetic clay lining [13], and composite liners [14, 15]. In order to prevent and slow down the damage of contaminant on the surrounding environment, it is necessary to study the related problems [16–18] caused by the diffusion of contaminants.

### 2. Solution of the Contaminant Transport Problem

### 2.1. Basic Assumptions and Mathematical Model

*L*

_{1},

*h*

_{1},

*h*

_{2}and

*L*

_{2}, respectively. The origin of the coordinate is at the top surface of contaminated soil layer 1 with the positive

*z*direction being downwards. The interfaces between the contaminated soil layers and the adjacent layers are modelled as the imperfect boundaries, i.e., the concentration at the boundaries between the contaminated soil layers and the adjacent layers is continuous, but the diffusion rate is discontinuous. With this kind of imperfect boundaries, not only can various soil layers with different diffusion capabilities be handled, but also can be reduced to the classical engineering problems. For instance, when the bottom boundary has zero diffusivity and the top diffusivity is weak, it can be transformed into a boundary condition commonly used in landfill storage problems. When the top and bottom of the contaminated layers are full diffusion boundary, it can be transformed into the boundary condition used for the problem of the contaminated soil diffusing to the adjacent water source.

##### (1)

$$\frac{\partial {C}_{i}(z,t)}{\partial t}={D}_{vi}\frac{{\partial}^{2}{C}_{i}(z,t)}{\partial {z}^{2}}$$*C*

*(*

_{i}*z*,

*t*) and

*D*

*are the concentration and effective diffusion coefficient of contaminant in the*

_{vi}*i*th (

*i =*1, 2) contaminated soil layer, respectively; For the

*i*th layer, the effective diffusion coefficient

*D*

*can be further defined as the product of the apparent tortuosity factor*

_{vi}*τ*

*and the aqueous diffusion coefficient for the solutes*

_{ai}*D*

_{0}

*i.e.,*

_{i}*D*

_{vi}*= D*

_{0}

_{i}*τ*

*[31];*

_{ai}*z*is the spatial coordinate variable along the depth direction, and

*t*is the time variable.

*C*

_{0}is the initial concentration of contaminant in contaminated soil layers 1 and 2 when

*t*= 0.

*z*=

*h*

_{1}, hence the following continuity conditions hold

##### (4)

$${n}_{1}{D}_{v1}\frac{\partial {C}_{1}({h}_{1},t)}{\partial z}={n}_{2}{D}_{v2}\frac{\partial {C}_{2}({h}_{1},t)}{\partial z}$$*n*

_{1}and

*n*

_{2}are the porosities of the contaminated soil layers 1 and 2, respectively.

*z*= 0) and the bottom surface of the contaminated soil layer 2 (i.e.,

*z*=

*H*) can be written as

*D*

_{1}(

*L*

_{1}) and

*D*

_{2}(

*L*

_{2}) are, respectively, the effective diffusion coefficients (thickness) of the top capping layer and the bottom natural soil layer;

*S*

_{d}_{1}and

*S*

_{d}_{2}are the imperfect diffusion capacity coefficients representing the ratio of the diffusion capacity of the imperfect diffusion layers (i.e., the top capping layer and the bottom natural soil layer).

*S*

*(*

_{di}*i*= 1, 2). When the capping layer has perfect protective function (ideal condition), it can be considered that it has no ability to diffuse, i.e,

*S*

_{d}_{1}= 0. However, when the protective function of capping layer is relatively weak, it has certain ability to diffuse, i.e.,

*S*

_{d}_{1}> 0. The above description also applies to the bottom natural soil layer.

### 2.2. General Solution

#### 2.2.1. Concentration in contaminated soil

*C*

*(*

_{i}*z*,

*t*) is introduced

##### (8)

$${C}_{i}(z,t)={\phi}_{i}(z){\theta}_{i}(t),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=1,2$$*φ*

*(*

_{i}*z*)

*θ*

*(*

_{i}*t*) at the same time results in

##### (9)

$$\frac{1}{{D}_{vi}{\theta}_{i}(t)}\frac{\text{d}{\theta}_{i}(t)}{\text{d}t}=\frac{1}{\phi (z)}\frac{{\text{d}}^{2}{\phi}_{i}(z)}{\text{d}{z}^{2}},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=1,2$$*λ*, Eq. (9) can be split into two separate ordinary differential equations as

##### (13)

$${\phi}_{i}(z)={B}_{i}\hspace{0.17em}\text{cos}(\sqrt{\lambda}z)+{F}_{i}\hspace{0.17em}\text{sin}(\sqrt{\lambda}z)$$##### (14)

$$\text{d}{\phi}_{i}(z)/\text{d}z=-{B}_{i}\sqrt{\lambda}\hspace{0.17em}\text{sin}(\sqrt{\lambda}z)+{F}_{i}\sqrt{\lambda}\hspace{0.17em}\text{cos}(\sqrt{\lambda}z)$$##### (15)

$$\begin{array}{c}{C}_{1}={C}_{0}\sum _{m=1}^{\infty}{f}_{m}[\text{sin}(\frac{{\lambda}_{m}}{H}z)+\frac{{\lambda}_{m}}{{S}_{d1}}\text{cos}(\frac{{\lambda}_{m}}{H}z)]{\text{e}}^{-{\lambda}_{m}^{2}{T}_{v}},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0\le z\le {h}_{1}\\ {C}_{2}={C}_{0}\sum _{m=1}^{\infty}{f}_{m}{Q}_{m}[j{\lambda}_{m}\hspace{0.17em}\text{cos}(j{\lambda}_{m}\frac{H-z}{H})+{S}_{d2}\hspace{0.17em}\text{sin}(j{\lambda}_{m}\frac{H-z}{H})]{\text{e}}^{-{\lambda}_{m}^{2}{T}_{v}},\\ {h}_{1}\le z\le H\end{array}$$##### (16)

$${f}_{m}=\frac{2{S}_{d1}^{2}(1+\sqrt{ab}{S}_{d2}{Q}_{m})}{{\lambda}_{m}\{k({\lambda}_{m}^{2}+{S}_{d1}^{2})+{S}_{d1}+b{S}_{d1}^{2}{Q}_{m}^{2}[(1-k)({j}^{2}{\lambda}_{m}^{2}+{S}_{d2}^{2})+{S}_{d2}]\}}$$##### (17)

$${Q}_{m}=\frac{{S}_{d1}\hspace{0.17em}\text{cos}(k{\lambda}_{m})-{\lambda}_{m}\hspace{0.17em}\text{sin}(k{\lambda}_{m})}{\sqrt{ab}{S}_{d1}\{j{\lambda}_{m}\hspace{0.17em}\text{sin}[(1-k)j{\lambda}_{m}]-{S}_{d2}\hspace{0.17em}\text{cos}[(1-k)j{\lambda}_{m}]\}}$$##### (18)

$$a=\frac{{n}_{2}{D}_{v2}}{{n}_{1}{D}_{v1}};b=\frac{{n}_{2}}{{n}_{1}};k=\frac{{h}_{1}}{H};j=\sqrt{{D}_{v1}/{D}_{v2}}=\sqrt{b/a};{T}_{v}=\frac{{D}_{v1}t}{{H}^{2}}$$*λ*

*(*

_{m}*m*= 1,2,3,…) is the positive root of the following characteristic equation

##### (19)

$$\begin{array}{c}\sqrt{ab}({\lambda}_{m}+{S}_{d1}\hspace{0.17em}\text{tan}(k{\lambda}_{m})\{j{\lambda}_{m}\hspace{0.17em}\text{tan}[j(1-k){\lambda}_{m}]-{S}_{d2}\}\\ -[{S}_{d1}-{\lambda}_{m}\hspace{0.17em}\text{tan}(k{\lambda}_{m})]\{j{\lambda}_{m}+{S}_{d2}\hspace{0.17em}\text{tan}[j(1-k){\lambda}_{m}]\})=0\end{array}$$*λ*

*in Eq. (19) can be solved using the dichotomization method. Besides, the series in Eq. (15) converge quickly, and the accuracy of the series can be ensured when setting*

_{m}*m*= 100.

*i*th (

*i*= 1, 2) layer caused by diffusion in contaminated soil can be expressed in terms of the first Fick’s law as

*J*

*(*

_{i}*z*,

*t*) is the mass flux in the

*i*th contaminated soil layer.

##### (21)

$$\begin{array}{c}{J}_{1}(z,t)=-{n}_{1}{D}_{v1}{C}_{0}\sum _{m=1}^{\infty}{f}_{m}[\frac{{\lambda}_{m}}{H}\text{cos}(\frac{{\lambda}_{m}}{H}z)-\frac{{\lambda}_{m}^{2}}{{S}_{d1}H}\text{sin}(\frac{{\lambda}_{m}}{H}z)]{\text{e}}^{-{\lambda}_{m}^{2}{T}_{v}},\\ 0\le z\le {h}_{1}\\ {J}_{2}(z,t)=-{n}_{2}{D}_{v2}{C}_{0}\sum _{m=1}^{\infty}{f}_{m}{Q}_{m}[\frac{{j}^{2}{\lambda}_{m}^{2}}{H}\text{sin}(j{\lambda}_{m}\frac{H-z}{H})-\frac{{S}_{d2}j{\lambda}_{m}}{H}\text{cos}(j{\lambda}_{m}\frac{H-z}{H})]{\text{e}}^{-{\lambda}_{m}^{2}{T}_{v}},\\ {h}_{1}\le z\le H\end{array}$$#### 2.2.2. Average degree of diffusion

*D*

*is defined as the ratio of the average concentration diffused to the adjacent layer to its initial concentration as [7].*

_{a}##### (22)

$$\begin{array}{l}{D}_{a}=\frac{{n}_{1}{\int}_{0}^{{h}_{1}}({C}_{0}-{C}_{1})\text{d}z+{n}_{2}{\int}_{{h}_{1}}^{H}({C}_{0}-{C}_{2})\text{d}z}{{C}_{0}({n}_{1}{h}_{1}+{n}_{2}{h}_{2})}\\ =\frac{k{D}_{d1}+b(1-k){D}_{d2}}{k+b(1-k)}\\ =1-\sum _{m=1}^{\infty}{f}_{m}\frac{1+\sqrt{ab}{S}_{d2}{Q}_{m}}{{\lambda}_{m}[k+b(1-k)]}{\text{e}}^{-{\lambda}_{m}^{2}{T}_{v}}\end{array}$$##### (23)

$$\begin{array}{l}{D}_{d1}=1-\frac{\frac{1}{{h}_{1}}{\int}_{0}^{{h}_{1}}{C}_{1}\text{d}z}{{C}_{0}}=1-\sum _{m=1}^{\infty}{f}_{m}\frac{{S}_{d1}[1-\text{cos}(k{\lambda}_{m})]+{\lambda}_{m}\hspace{0.17em}\text{sin}(k{\lambda}_{m})}{{S}_{d1}k{\lambda}_{m}}{\text{e}}^{-{\lambda}_{m}^{2}{T}_{v}}\\ {D}_{d2}=1-\frac{\frac{1}{{h}_{2}}{\int}_{{h}_{1}}^{H}{C}_{2}\text{d}z}{{C}_{0}}=1-\sum _{m=1}^{\infty}{f}_{m}\frac{{S}_{d2}\{1-\text{cos}[j(1-k){\lambda}_{m}]\}+{\lambda}_{m}\hspace{0.17em}\text{sin}[j\text{(1}-k\text{)}{\lambda}_{m}]}{(1-k)j{\lambda}_{m}}{Q}_{m}{\text{e}}^{-{\lambda}_{m}^{2}{T}_{v}}\end{array}$$*D*

*(*

_{di}*i*= 1, 2) is the average degree of diffusion of contaminant in the

*i*th contaminated soil layer;

*D*

*is the average degree of diffusion, which describes the overall degree of diffusion of contaminants in the contaminated soil layers.*

_{a}*S*

*= 0 and ∞ (*

_{di}*i*= 1, 2) are the ideal boundaries, which correspond to the full diffusion and non-diffusion boundaries, respectively. When

*S*

*are equal to the finite value, the boundaries are really the imperfect boundaries. When*

_{di}*S*

_{d}_{1}is equal to 0 and the value of

*S*

_{d}_{2}represents the diffusion capacity of the boundary between the liner and the contaminated soil layer, the calculation model can be transformed into the corresponding engineering problem in the landfill [11, 13]. When

*S*

_{d}_{2}is equal to 0 and the value of

*S*

_{d}_{1}represents the diffusion capacity of the boundary between the capping layer and the contaminated soil layer, the calculation model can be transformed into the engineering problem corresponding to the underwater capping layer [12]. In addition, the diffusion capacity of the boundary layer increases with the increase of

*S*

*.*

_{di}### 3. Results and Discussion

### 3.1. Verification of the Present Solution

*H*= 1.5 m,

*D*

*= 9.4 × 10*

_{v}^{−10}m

^{2}/s,

*n*= 0.45 and retardation factor

*R*

*= 43.3. It is noted that Li and Cleall [31] considered the effect of retardation factor, hence we need to modify*

_{d}*D*

*as*

_{v}*D*

*/*

_{v}*R*

*to keep the same condition. Hence, the input parameters in the present calculation are:*

_{d}*h*

_{1}=

*h*

_{2}= 0.75 m,

*D*

_{v}_{1}=

*D*

_{v}_{2}= 2.1709 × 10

^{−11}m

^{2}/s,

*n*

_{1}=

*n*

_{2}= 0.45. For general applications, the normalized contaminant concentration is defined as

*C*

^{*}=

*C*/

*C*

_{0}, in which

*C*denotes the concentration of contaminant in contaminated layers 1 and 2. Fig. 2 shows the comparison of contaminant concentration vs. normalized depth between the reduced present solution and Li and Cleall’s solution. It can be seen from Fig. 2 that the results from the present solution are very close to those from Li and Cleall [31] for different times.

_{4}and NaCl contaminant sources with the same initial contaminant concentration. The effective diffusion coefficients of Zn

^{2+}(i.e.,

*D*

_{v}_{1}) and Cl

^{−}(i.e.,

*D*

_{v}_{2}) are employed in the following calculation [35]. Four cases with different combination of

*S*

_{d}_{1}and

*S*

_{d}_{2}are: Case 1 with

*S*

_{d}_{1}= ∞,

*S*

_{d}_{2}= ∞; Case 2 with

*S*

_{d}_{1}= ∞,

*S*

_{d}_{2}= 1; Case 3 with

*S*

_{d}_{1}= 15,

*S*

_{d}_{2}= 0; Case 4 with

*S*

_{d}_{1}= 5,

*S*

_{d}_{2}= 0. Unless otherwise specified, the basic parameters used in the following calculations (i.e., Fig. 3 to Fig. 8) are:

*H*= 1m,

*h*

_{1}= 0.5 m,

*T*

*= 0.129 or 0.279,*

_{v}*C*

_{0}= 0.1 mol/L,

*D*

_{v}_{1}= 3.3 × 10

^{−10}m

^{2}/s,

*D*

_{v}_{2}= 6.8 × 10

^{−10}m

^{2}/s,

*n*

_{1}= 0.25,

*n*

_{2}= 0.5. It should be pointed out again that

*S*

*= 0, ∞ and finite values (*

_{di}*i*= 1, 2) represent the non-diffusion, full diffusion and imperfect diffusion boundaries.

*D*

*in the*

_{ri}*i*th layer can be defined as [36].

*D*

*is the mechanical dispersion coefficient in the*

_{mi}*i*th layer (

*i*= 1, 2), and it can be defined as

*D*

*=*

_{mi}*h*

_{i}*v*

*/10 with*

_{i}*h*

*and*

_{i}*v*

*being the thickness and seepage velocity of the*

_{i}*i*th layer, respectively.

*K*

*is the harmonic average of the permeability coefficient in each layer [36];*

_{a}*H*

*and*

_{w}*n*

*are, respectively, the water head height and porosity of the*

_{i}*i*th layer (

*i*= 1, 2).

*v*

_{i}*h*

*/*

_{i}*D*

*≤ 1 [37]. The additional selected parameters used in numerical example B are:*

_{vi}*H*

*= 0.3 m,*

_{w}*K*

*= 1.27 × 10*

_{a}^{−10}m/s,

*v*

_{1}= 6.6 × 10

^{−10}m/s,

*v*

_{2}= 3.3 × 10

^{−10}m/s.

*C*

^{*}in numerical example B is smaller than that in numerical example A. Moreover, as the imperfect diffusion capacity coefficients increases, the convection and mechanical dispersion contributes less to transport of contaminant.

### 3.2. Influence of Imperfect Diffusion Boundaries

*S*

_{d}_{1}and

*S*

_{d}_{2}). It can be observed from Fig. 4 that

*S*

*(*

_{di}*i*= 1, 2) have great influence on the concentration distribution along the depth direction. When

*S*

*= 0 (i.e., ideal non-diffusion boundary), the normalized contaminant concentration*

_{di}*C*

^{*}is equal to 1, which means that no contaminant diffuses to the adjacent layer. However, this is very difficult to realize in the engineering practice. When

*S*

*= ∞ (i.e., ideal full diffusion boundary), the concentrations at*

_{di}*z*= 0 and

*H*are equal to zero, indicating that the contaminants at the top surface of the contaminated soil layer 1 and at the bottom surface of the contaminated soil layer 2 can freely transport to the top capping and bottom natural soil layers. Hence, in this case (i.e.,

*S*

*= ∞),*

_{di}*C*

^{*}is the smallest when

*t*= 6.0155 years. As

*S*

_{d}_{1}or

*S*

_{d}_{2}decreases,

*C*

^{*}shows marked increase within the entire depth range. This indicates that the diffusion capacity of the capping and bottom natural soil layers decreases markedly with decreasing

*S*

_{d}_{1}or

*S*

_{d}_{2}, and more contaminants are kept in the contaminated soil layers. It can be also seen from Fig. 4 that there exists obvious change of curve slope at the interface of two contaminated soil layers. The greater the difference between the values of

*S*

_{d}_{1}and

*S*

_{d}_{2}, the more obvious the curve slope changes at the interface.

*T*

*is relatively small, the contaminant concentration at fixed level (i.e.,*

_{v}*z*= 0.5 m) decreases greatly with increasing

*S*

*, especially when*

_{di}*S*

*= ∞. It can be also observed from Fig. 5 that the decreasing ratio of the contaminant concentration in the contaminated soil decreases with increasing time factor. This phenomenon is in accordance with the basic law, i.e., the diffusion velocity of contaminant is affected by the concentration gradient.*

_{di}*x*-axis is in a base-10 logarithmic scale. It can be observed from Fig. 6 that the average degree of diffusion

*D*

*shows marked increase as time increases. As the time further increases, the changing rate of the average degree of diffusion decreases due to the decrease of the concentration gradient (see the upper four curves in Fig. 6). It can be also seen from Fig. 6 that for a given time, the larger the value of the combination of*

_{a}*S*

*, the greater the value of the average degree of diffusion*

_{di}*D*

*. When*

_{a}*S*

_{d}_{1}=

*S*

_{d}_{2}= ∞, the top and bottom surfaces of the contaminated soil layers are ideally full diffusion boundaries, and completing the diffusion process needs only 23.32 years. However, when

*S*

_{d}_{2}= 0, due to the influence of the ideally non-diffusion boundary of contaminated soil layer 2, the contaminant can only diffuse through the top capping layer. Accordingly, it needs at least 69.95 years to complete the entire diffusion process. It can be concluded from Fig. 4~Fig. 6 that the imperfect diffusion capacity coefficients have great effect on the diffusion of the contaminant. Therefore, in engineering practice, the intrinsic diffusion capacity of the capping layer should be carefully designed for the better environmental conservation.

*z*= 0 and

*z*=

*H*level. It can be observed from Fig. 7 and Fig. 8 that the imperfect diffusion boundaries have great influence on the variation of the mass flux with increasing time. In the initial stage of the diffusion (i.e.,

*T*

*is relatively small), the mass flux at the top and bottom level of the contaminated soil layers increases with the increase of*

_{v}*S*

_{d}_{1}or

*S*

_{d}_{2}. However, when the time factor

*T*

*further increases, the corresponding mass flux decreases with the increase of*

_{v}*S*

_{d}_{1}or

*S*

_{d}_{2}. The main reason for this phenomenon is that in the initial stage, the contaminant concentration in the soil decreases quickly for a larger

*S*

_{d}_{1}or

*S*

_{d}_{2}value, which induces the decreasing concentration difference in contaminated soil. As a result, when the time is relatively big, the mass flux corresponding to larger

*S*

_{d}_{1}or

*S*

_{d}_{2}value is smaller than that corresponding to smaller

*S*

_{d}_{1}or

*S*

_{d}_{2}value.

### 4. Conclusions

*S*

_{d}_{1}and

*S*

_{d}_{2}) on the spatial/time distribution of contaminant concentration, average degree of diffusion and the mass flux is analyzed. The theory developed in this paper reveals the diffusion mechanism of the contaminant from the double-layered contaminated soil layers to the adjacent layers considering the influence of the imperfect diffusion boundaries, which can provide useful guidance for the engineering practice and could be further served as benchmarks for future numerical solutions. The main conclusions are as follows:

*S*

_{d}_{1}and

*S*

_{d}_{2}.

*S*

_{d}_{1}or

*S*

_{d}_{2}decreases, contaminant concentration shows marked increase within the entire depth range for a given time. There also exists obvious change of curve slope at the interface of two contaminated soil layers due to the layering of the soil. When time factor

*T*

*is relatively small, the contaminant concentration at a fixed depth of the contaminated soil layer decreases greatly with increasing*

_{v}*S*

_{d}_{1}or

*S*

_{d}_{2}for a given time.

*S*

_{d}_{1}and

*S*

_{d}_{2}have great influence on the average degree of diffusion and mass flux. The average degree of diffusion increases with increasing

*S*

_{d}_{1}or

*S*

_{d}_{2}for a given time. The mass flux decreases greatly with increasing

*S*

_{d}_{1}or

*S*

_{d}_{2}in the initial stage, while it increases with increasing

*S*

_{d}_{1}or

*S*

_{d}_{2}when time factor

*T*

*is relatively large.*

_{v}