Abstract

A sand-box experiment equipment is constructed to simulate the transportation of contaminant in unsaturated flow under dirichlet boundary condition. Based on analysis of the mechanism in unsaturated migration, taking a riverland as an example, the model of water quality which was affected by contaminant transport was built, and the calculation values coincided well with experimental data. The results showed that the lateral migration should not be neglected when seepage velocity was high. In addition, the lateral migration was faster in high moisture content areas and the solute was retained for a short time; however, the situation was reversed in low moisture content areas. Moreover, during the initial stage of solute transport, the velocity of the trailing end of the pinnate pollution zone was high under the effect of capillary action, which would lead to expansion of the range of pollution. Furthermore, the relative concentration curve of pollutants presented a single peak state, with the peak value of the curve increasing with the velocity, and the arrival time of the peak value was also gradually advanced. Therefore, this study can provide reference for the prediction and control of groundwater contaminant transport in an unsaturated zone.

HIGHLIGHTS

  • The characteristics of contaminant transport in unsaturated groundwater are explored.

  • Water content and lateral migration velocity affect the retention time of solute.

  • The relative concentration curve is an unsymmetrical single peak curve with obvious trailing phenomenon.

  • The velocity of the trailing end of the pinnate pollution zone was high during the initial stage of solute transport.

Graphical Abstract

Graphical Abstract
Graphical Abstract

INTRODUCTION

In recent years, with the rapid development of China's society and economy, the demand for water resources has been constantly increasing. The excessive exploitation of groundwater and pollution have led to the deterioration of the geological environment, making the groundwater pollution situation in individual areas more serious and some ecological damage even more difficult to reverse, so the problem of groundwater pollution has been getting more and more attention. There are some common contaminants such as heavy metal ions in industrial wastewater, pesticides and their degradation products all over the world (Cerejeira et al. 2003). In 2005 to 2008, it was reported that more than 20% of groundwater supplies were identified as being contaminated with pesticides in the Netherlands, France and Germany (Morvan et al. 2006).

Groundwater pollution is mainly manifested in various physical and chemical reactions between groundwater and solute and between solute and medium, resulting in pollutant migration in soil and aquifer. The kinematics of mass transport phenomenon in soil media are accomplished by advection, dispersion, sorption, chemical reactions and decay (Chakraborty et al. 2015). It is significant to prevent and control groundwater pollution and utilize water resources via methods to calculate and predict groundwater solute transport patterns (Hanasaki et al. 2016). The contaminants produced by human activities are mainly distributed on the earth's surface or buried in shallow soil and enter into phreatic aquifer with rainfall infiltration through the unsaturated zone (Vero et al. 2017). Because the migration and transformation of contaminants in groundwater is slow and complex, scholars have done a lot of research on the problem of unsaturated and saturated solute transport. Lee et al. (2001) discussed the relationship between simulated grid-scale and fracture scale, analyzed the main factors that influence solute transport and the basic equation of average concentration of pollutants; Scher et al. (2002) established a model of solute transport probability in heterogeneous porous medium based on the Continuous time random walk (CTRW) theory, and the model provided an effective simulation method for solute transport in geological media.

However, in the regional simulation of a stratified aquifer, the key lies in the simplification of three-dimensional solute transport equation. When the discrete analysis and numerical solution of the model are considered, it needs to figure out that how to avoid the calculation error caused by the difference of the unit scale (Bergvall et al. 2011). Zhang et al. (2016) used Monte Carlo method to simulate various possible heterogeneous spatial distributions of aquifer parameters, analyzed the random field that generated a large number of permeability coefficients, and studied the migration and transformation of pollutants in groundwater. Therefore, on the basis of the uncertainty of the parameters of the solute migration model in porous media, factors such as hydrogeological conditions, types of pollution sources and concentration of pollutants were taken into account, providing a way of thinking for the analysis of the problem (Sun et al. 2014).

The above research mainly predicted quantitatively and analyzed the solute transport patterns through various mathematical models (Mumford et al. 2016; Sieczka et al. 2018), and the numerical calculation method was used to convert the contaminant transport in the unsaturated zone into vertical one-dimensional flow; In the phreatic aquifer, it was generalized to the horizontal two-dimensional flow. Nevertheless, as the contaminants buried in a shallow layer near dirichlet boundary such as fluvial boundary, the horizontal migration patterns of which entering the unsaturated zone with the rising water level of aquifer were less studied.

Groundwater in the unsaturated zone is the main carrier of solute transport, and solute is transported with the flow. Understanding and describing solute transport process in unsaturated zone is a key scientific problem in groundwater pollution control and river ecosystem restoration. In this paper, aiming at the situation of the contaminant entering the unsaturated zone near the boundary when the river level rises, considering factors such as water content, flow rate and migration time, using the indoor sand-box experiment and numerical model simulation method, respectively, the results of the experiment and numerical calculation are compared and analyzed to verify the rule of solute transport. It not only provides a scientific basis for prevention and prediction of groundwater contaminant transport in an unsaturated zone, and more clearly understand the transport mechanism, but also has practical effect and significance for later control and treatment.

MATERIAL AND METHODS

Physics experiment

From the experimental point of view, the hydrodynamic dispersion coefficient is calculated by using the soil column experiment, and the solute migration is simulated in the homogeneous and heterogeneous soil columns (Lee et al. 2015). Besides, sand-box experiment equipment is constructed to simulate solute transport in the unsaturated zone near the river boundary, which is made of transparent organic glass. The main part size is 1,200 mm × 200 mm × 1,000 mm, and the middle piece of sand-box with filling medium is 1,000 mm long, the structure of the device is shown in Figure 1. The fluctuation of river water level is simulated by sand-box head and tail, and the potassium permanganate solution at the boundary is used as contaminant to observe the migration trajectory and change of concentration.

Figure 1

The structure schematic of the experimental sand-box.

Figure 1

The structure schematic of the experimental sand-box.

The experiment simulates the contaminant transport with the flow while the water level of aquifer rises with the river level. Before the start of the experiment, adjust the regulative water tank of the head and tail to make the water level consensus, leave for one day to wait for the level of the medium to be horizontal. In the beginning, the potassium permanganate solution is injected into the boundary of the sand-box head, then open the water pump and quickly adjust the water level of the head tank to the position of potassium permanganate solution injection height, and stays the same. At last, the migration of potassium permanganate solution in the medium is recorded by camera.

Numerical model

For the results of the experiment, validation of mathematical models is necessary. It can be described by the Richards equation based on Darcy's law for the saturated-unsaturated seepage flow. The Richards equation under pressure head form (Tegnander 2001) is as follows:
formula
(1)
where h is pressure head (L); kij is the hydraulic conductivity of saturated media (LT−1); kr is the relative hydraulic conductivity which is a function of saturation (0kr1); C=dθ/dh is the moisture content, which is 0 in the saturated zone (L−1); β is 1 in the saturated zone and 0 in the unsaturated zone; Ss is specific yield (L−1); qw is source-sink term (plus is source term and minus is sink term).
It is assumed that the transverse component of the contaminant transport in the unsaturated zone near the river can't be ignored, and the horizontal component should be included in the flow equation and solute transport equation. Equation (1) is changed into the following form:
formula
(2)
Considering the effects of advection, dispersion, adsorption and degradation (Allen-King et al. 1998; Shi et al. 2012; Natarajan 2014), the equation of contaminant transport and transformation is as follows:
formula
(3)
where vi is the Darcy velocity (cm/h); Dij is the dispersion coefficient of the contaminant (cm2/h). The third term on the right-hand side is the adsorption term, ρb is the dry density of the medium (g/cm3); S is the mass component of solute in the adsorbed phase (μg/g). The fourth item on the right is degradation, where λa and λs are the first order degressive degradation coefficients of the liquid and solid phase. The last term on the right-hand side is the source-sink, cs is the contaminant concentration of source (sink) (μg/cm3).
The complete expression of the convection-dispersion equation is shown in Equation (4), where the first three terms at the right represent solute transport due to hydrodynamic dispersion, and the middle three terms represent solute transport due to water convection, the last term refers to the consideration of the effects of chemical reactions or other factors on solute transport (Cui et al. 2017).
formula
(4)
where Dxx, Dyy and Dzz are the components of hydrodynamic dispersion coefficient D in the X, Y and Z directions, respectively; vx, vy, vz are the velocity components of V in the convective motion of water flow in the X, Y and Z directions, respectively; f represents the change of solute mass caused by other reactions or factors in unit body at Δt time.
When the medium is in one-dimensional motion, the convection-dispersion equation can be simplified to Equation (5):
formula
(5)

The above Equation (5) is the first and second terms on the right side of Equation (4) of the contaminant transport equation. In addition, the solution of the convective dispersion equation is the concentration distribution of solute. If the solution is required, some other conditions should be given, including the spatial region and time range of the study area, the distribution of the hydraulic field, the relevant parameters of the dispersion degree and the corresponding definite solution conditions.

To sum up, the groundwater flow equation is coupled to the solute transport equation by the Darcy velocity, and the mathematical model of pollutant migration in an unsaturated – saturated medium with initial conditions and boundary conditions is added. The solution process is as follows: firstly, the groundwater seepage equation should be solved, then the pressure head h is obtained, and the concentration distribution of the contaminant is obtained by using the flow velocity v.

RESULTS AND DISCUSSION

Experimental results

The model requires that the medium should have homogeneous and isotropic property. In order to satisfy this condition, the original sand sample is screened, which makes the sand sample uniform and reduces the content of mud. When filling the medium, it is uniformly scattered, layered, and filled in layers. In addition, the density of standard sand is selected from 1.47–1.61 g/cm3; the physical properties of sand samples are shown in Table 1.

Table 1

Analysis of particle sizes of sand samples

Particle sizes (mm)>1.00.5–1.00.25–0.50.075–0.25<0.075
Weight percentage 30% 66% 4% 
Particle sizes (mm)>1.00.5–1.00.25–0.50.075–0.25<0.075
Weight percentage 30% 66% 4% 
The experiment adopts the water characteristic curve of the Van Genuchten model fitting medium, which is shown as follows:
formula
(6)
where θr is the residual moisture content, θs is the saturated volume moisture content, α and n are the fitting parameters of the Van Genuchten curve, m=11/n. The hydrogeological parameters of experimental sand samples are shown in Table 2.
Table 2

Hydrogeological parameters of sand samples

Saturated hydraulic conductivity (cm/s)Specific yieldPorosityθrθsαn
0.02 0.105 0.03 0.398 0.148 0.15 3.35 
Saturated hydraulic conductivity (cm/s)Specific yieldPorosityθrθsαn
0.02 0.105 0.03 0.398 0.148 0.15 3.35 

The initial water level of the experiment is set to 300 mm, and the change of water level of the sand-box head is 200 mm. The concentration distributions of solute transport recorded in the course of the experiment are shown in Figure 2.

Figure 2

The concentration distribution of solute transport at different times. (a) 5 min. (b) 10 min. (c) 20 min. (d) 40 min. (e) 1,200 min. (f) 1,440 min.

Figure 2

The concentration distribution of solute transport at different times. (a) 5 min. (b) 10 min. (c) 20 min. (d) 40 min. (e) 1,200 min. (f) 1,440 min.

As shown in Figure 2, when the water level of the sand-box rises to the position of potassium permanganate solution, the solute is transported with unsaturated seepage. In the absence of infiltration supplies and due to the convection of seepage flow, the contaminant is in lateral migration in the unsaturated zone without entering the saturated zone. In the initial stage of contaminant transport, the velocity of the trailing end of the pinnate pollution zone is high under the effect of capillary action, and makes the pollution dispersed. Dissolved phase plumes typically form at contaminated sites, which has been studied in some literatures (Robertson et al. 2016). Contaminant transport in the unsaturated zone can be divided into two different areas: the velocity of contaminant transport is relatively fast in the high moisture content area, and as the water level rises, the contaminant is transferred to some areas with lower moisture content by the convective effect of seepage; in the low moisture content area, the contaminant transport rate is slower, and lasts for a long time. The above characteristics are basically consistent with the conclusion of some papers (Sheng et al. 2011; Zhang et al. 2018). In this experiment, the contaminant in the unsaturated region with high moisture content reaches the end of the river at 40 min, and then rapidly flows into the river with the seepage (Figure 2(d)); Nevertheless, the contaminant in the area with low moisture content can be detected after 1,440 min.

Experimental results demonstrated that the lateral migration component of the solute was larger in the unsaturated zone with high moisture content. When the phreatic aquifer flow speed was faster, the unsaturated zone solute transport could not be completely generalized to one-dimensional vertical motion.

Analysis of model application

Taking the interfluve as an example, it is assumed that the landmass and soil homogeneity are isotropic, which extends infinitely in the lateral space with a horizontal water bottom plate. The river cuts through the whole phreatic aquifer, and then the water level rises instantaneously on the left side of the river, which has increased Δh and remains unchanged later, Δh0,t=h0,t-h0,0; the initial free surface of phreatic water is horizontal, and the hydraulic gradient is zero, h=h(x, 0), and the main infiltration in the surface or evaporation is zero too. It ignores the adsorption and degradation of contaminants, and the calculated coordinates and contaminant locations are shown in Figure 3.

Figure 3

The sketch of unsaturated contaminant transport in an interfluve.

Figure 3

The sketch of unsaturated contaminant transport in an interfluve.

In the experiment of solute transport in groundwater, the flow state in the sand-box is analyzed according to the measured hydraulic gradient and velocity data, so as to reveal the concentration change of solute transport in the sand layer. When the water flow is laminar, the hydraulic gradient is linearly correlated with the corresponding velocity, which conforms to Darcy's law in the process of one-dimensional porous medium migration. With the increase of the flow rate, Darcy's law is no longer applicable to the experimental conditions when the flow is turbulent. At this time, the liquid flow follows the nonlinear seepage equation, and its fitting curve is shown in Figure 4.

Figure 4

Nonlinear relationship between flow velocity and hydraulic gradient.

Figure 4

Nonlinear relationship between flow velocity and hydraulic gradient.

In this case, the solubility of the contaminant is large, and the absorption is small. The equation plus initial conditions and boundary conditions are in the above situation, so that the mathematical flow model of the problem can be described as:
formula
(7)
formula
(8)
formula
(9)
formula
(10)
The solute mathematical model can be described as:
formula
(11)
formula
(12)
formula
(13)
formula
(14)
formula
(15)
where C* is the relative concentration, C*=C/C0, C0 is the initial concentration of the contaminant. With regard to the model of Equations (1) and (3), some scholars have proposed different solutions and compared the advantages and disadvantages of different algorithms (Kheirabadi et al. 2017), which are also suitable for the model of Equations (2) and (3).

In this paper, we adopted the methods of combining forward weighted difference and implicit finite difference to solve the model of Equations (5)–(8) and Equations (9)–(13). In the above equations, L=1,000 mm, H=1,000 mm, h0,0=300 mm, Δh=200 mm, the soil parameters are the actual values of the sand-box experimental medium. The change rules of contaminant transport are simulated in unsaturated and saturated zones with water flow. When the flow velocity is varied, the curve of solute relative concentration over time is shown in Figure 5. The comparison between the simulation results and the calculated values of the one-dimensional migration model is shown in Figure 6.

Figure 5

Solute relative concentration distribution curves at different flow rates. (a) v1=0.789 cm/s (b) v2=1.421 cm/s. (c) v3 = 2.112 cm/s (d) v4 = 2.685 cm/s.

Figure 5

Solute relative concentration distribution curves at different flow rates. (a) v1=0.789 cm/s (b) v2=1.421 cm/s. (c) v3 = 2.112 cm/s (d) v4 = 2.685 cm/s.

Figure 6

The comparison of concentration distribution between calculated values and experimental values (t = 40 min).

Figure 6

The comparison of concentration distribution between calculated values and experimental values (t = 40 min).

As can be seen from Figure 5, under different hydraulic conditions, the relative concentration curve of pollutants presents a single-peak curve. The peak value of the curve increases with the increase of velocity, and the arrival time of the peak value gradually advances with the increase of velocity. In addition, the relative concentration curves of pollutants under the four working conditions have obvious trailing phenomenon. When the flow velocity is low, the main part of the curve is wider, and with the increase of the flow velocity, the main part of the curve is gradually narrowed. The above research is consistent with the conclusions of relevant literature (Ye et al. 2005; Ji et al. 2017).

Obviously, the simulated results coincide well with the experimental data. Figure 6 shows that the one-dimensional migration model transforms the unsaturated zones' contaminant transport into one-dimensional vertical motion, which ignores the transverse component of seepage in unsaturated zone and calculates mainly the distribution of a contaminant in aquifers, and it is inconsistent with reality. In addition, due to the velocity of the aquifer being higher than that of the unsaturated zone, the contaminant spreads faster, and the concentration of contaminants in the same section is lower than the actual value.

In summary, the flow velocity is faster in the phreatic aquifer, and the lateral convection is stronger with high moisture content, the lateral migration of contaminant unsaturated zone cannot be ignored in the upper infiltration condition.

CONCLUSIONS

  • (1)

    In this paper, the mathematical model of unsaturated groundwater and contaminant solute transport is established, which is verified by laboratory sand-box experiment, and the simulation results are basically consistent with the experimental data. In the initial stage of solute transport, the velocity of the trailing end of the pinnate pollution zone was high under the effect of seepage capillary action in the unsaturated zone, which made the pollution expand.

  • (2)

    The lateral migration of contaminants in the unsaturated zone cannot be ignored when the flow velocity is fast in the phreatic aquifer. In the region with higher water content, the lateral migration velocity is faster and the contaminant retention time is shorter, while in the region with lower water content, the lateral migration velocity is slower and the contaminant retention time is longer.

  • (3)

    In the process of solute transport, the relative concentration curve of pollutants presents a single peak curve, but there are some differences among different hydraulic conditions. The peak value of the curve increases with the velocity, and the arrival time of the peak value is slowly advanced with the increase of the velocity. The main part of the curve is also gradually narrowed with the velocity.

ACKNOWLEDGEMENTS

This research work was supported by the key project of the Natural Science Foundation of Anhui (Grant No. KJ2017A409), partly supported by the National Natural Science Foundation of China (Grant No. 51309071), and supported by the Central Government Guided Local Special Funds for Science and Technology projects (Grant No. 2016080802D118). We would also like to thank the editor and reviewers for their valuable and useful comments and suggestions.

CONFLICTS OF INTEREST

The authors declare no conflict of interest.

DATA AVAILABILITY STATEMENT

All relevant data are included in the paper or its Supplementary Information.

REFERENCES

Bergvall
M.
Grip
H.
Sjostrom
J.
Laudon
H.
2011
Modeling subsurface transport in extensive glaciofluvial and littoral sediments to remediate a municipal drinking water aquifer
.
Hydrology and Earth System Sciences
15
(
7
),
2229
2244
.
Cerejeira
M. J.
Viana
P.
Batista
S.
Pereira
T.
Silva
E.
Valerio
M. J.
Silva
A.
Ferreira
M.
Silva-Fernandes
A. M.
2003
Pesticides in Portuguese surface and ground waters
.
Water Research
37
,
1055
1063
.
Chakraborty
R.
Ghosh
A.
Ghosh
S.
Mukherjee
S.
2015
Evaluation of contaminant transport parameters for hexavalent chromium migration through saturated soil media
.
Environmental Earth Sciences
74
,
5687
5697
.
Hanasaki
N.
Yoshikawa
S.
Kakinuma
K.
Kanae
S.
2016
A seawater desalination scheme for global hydrological models
.
Hydrology and Earth System Sciences
20
,
4143
4157
.
Ji
S. S.
Liu
J. G.
Wu
Y.
Peng
D. W.
Zhang
B. L.
Wang
M. X.
2017
Comparative study on characteristics of tracer curve laboratory experiment of karst pipeline and karst fracture
.
Site Investigation Science and Technology
4
,
11
14
.
Lee
S. H.
Loudh
M. F.
Jensen
C. L.
2001
Hierarchical modeling of flow in naturally fractured formations with multiple length scales
.
Water Resources Research
37
(
3
),
443
455
.
Lee
C. P.
Wu
M. C.
Tsai
S. C.
Liu
C. Y.
Tsai
T. L.
Pan
C. H.
Yang
M. S.
2015
Numerical analysis of transport and retardation for cesium in crushed granite using multi-stage advection–dispersion column experiments
.
Journal of Radioanalytical and Nuclear Chemistry
304
,
377
386
.
Mumford
K. G.
Mustafa
N.
Gerhard
J. I.
2016
Probabilistic risk assessment of contaminant transport in groundwater and vapour intrusion following remediation of a contaminant source
.
Stochastic Environmental Research and Risk Assessment
30
,
1017
1031
.
Robertson
W. D.
Van Stempvoort
D. R.
Roy
J. W.
Brown
S. J.
Spoelstra
J.
Schiff
S. L.
Rudolph
D. R.
Danielescu
S.
Graham
G.
2016
Use of an artificial sweetener to identify sources of groundwater nitrate contamination
.
Groundwater
54
,
579
587
.
Scher
H.
Margolin
G.
Berkowitz
B.
2002
Towards a unified framework for anomalous transport in heterogeneous media
.
Chemical Physics
284
(
1/2
),
349
359
.
Sheng
F.
Wang
K.
Zhang
R. D.
Liu
H. H.
2011
Modeling preferential water flow and solute transport in unsaturated soil using the active region model
.
Environmental Earth Sciences
62
,
1491
1501
.
Shi
X. Q.
Wu
J. C.
Wu
J. F.
Jiang
B. L.
Xu
H. X.
2012
Effects of the heterogeneity of multiple correlated random parameters on solute transport
.
Advances in Water Science
23
(
1
),
509
515
.
Sun
F. S.
Cheng
P.
Zhang
B.
2014
Physical process based risk assessment of groundwater pollution in the mining area
.
Environmental Science
35
(
4
),
1285
1289
.
Vero
S. E.
Healy
M. G.
Henry
T.
Creamer
R. E.
Ibrahim
T. G.
Richards
K. G.
Mellander
P. E.
McDonald
N. T.
Fenton
O.
2017
A framework for determining unsaturated zone water quality time lags at catchment scale
.
Agriculture, Ecosystems & Environment
236
,
234
242
.
Ye
H.
Qian
J. Z.
Huang
X. C.
Li
Y. L.
Dong
H. X.
2005
Application of projection pursuit model to evaluation of groundwater quality
.
Hydrogeology and Engineering Geology
5
,
9
11
.
Zhang
B.
Li
G. X.
Cheng
P.
Sun
F. S.
Wang
H. Y.
2016
Groundwater environment risk assessment based on stochastic theory
.
Advances in Water Science
27
(
1
),
100
106
.
Zhang
Y. H.
Xu
X. F.
Liu
Y.
Si
G. H.
Wang
Q. H.
2018
Influence of water content on permeability and dispersion characteristics of unsaturated clay and radionuclide migration
.
Journal of Engineering Geology
26
(
S
),
657
664
.