Fine Characterization of the Effects of Aquifer Heterogeneity on Solute Transport: A Numerical Sandbox Experiment

Abstract

Hydraulic conductivity (K) and the specific storage coefficient (S) are among the most important hydrogeological parameters of an aquifer. Traditionally, the hydrogeological parameters of a field aquifer system are mainly determined through a range of experiments that are both time-consuming and of poor operability. To accurately characterize aquifer heterogeneity, a synthetic sandbox is constructed using VSAFT2 (Variably Saturated Flow and Transport utilizing the Modified Method of Characteristics, in 2D) as a reference aquifer system by incorporating multilevel a priori geologic information into the sandbox configuration. The spatial distribution of the field of hydraulic conductivity (i.e., K) is inversely obtained by hydraulic tomography (HT). Then HT is compared with traditional kriging-estimated method in the fine characterization of aquifer heterogeneity, and the optimal K field is eventually selected to predict the solute transport. The influence of the number of pumping cycles on the accuracy of heterogeneity characterization is also discussed. The results show that the accuracy of the inversely obtained K field is improved with the increased number of pumping cycles. When incorporating multilevel a priori geological information, HT can characterize aquifer heterogeneity more finely than traditional kriging, and there is also a very good fitting of solute transport between the optimally estimated K field and the reference K field. Our study highlights the importance of the fine characterization of aquifer heterogeneity for the prediction of solute transport.

1. Introduction

Groundwater contamination and remediation have become very important issues. The heterogeneous nature of geologic bodies results in the heterogeneous hydrological properties of geologic media accordingly [1]. The complex spatial structure of aquifer media leads to the uncertainty of groundwater flow, thereby aggravating the uncertainty of groundwater contaminant transport [2]. Previous studies have shown that the accuracy of the spatial distribution of hydrogeological parameters directly affects the hydrogeological evaluation and prediction when solute transport is considered [3,4,5]. Therefore, it is particularly important to characterize the aquifer heterogeneity with high accuracy.

Traditional methods for determining the hydrogeological parameters of field aquifers mainly include the pumping test, injection test, slug test, and water level recovery test [6]. However, the core samples and pumping test data obtained from field sampling are a bit limited to a certain extent. In recent decades, the efficient use of field and experimental data to invert the spatial distribution of hydraulic conductivity (K) has been widely studied by hydrogeologists [7,8]. At the end of the 20th century, Yeh et al. [9] developed a hydraulic tomography (HT) method that was different from the traditional methods. By incorporating the computerized axial tomography (CAT) and the geophysical tomography concepts, HT uses the inversion model VSAFT2 [10] to obtain the spatial distribution of K field, thereby characterizing aquifer heterogeneity with high accuracy.

To date, numerous satisfactory results have been achieved in the characterization of aquifer heterogeneity using HT, mainly in the spatial estimation of K field and the specific storage coefficient (S). For example, Illman et al. [11] and Liu et al. [12] collected the hydraulic head response data for the steady state portion of a pumping test using a laboratory sandbox experiment and then estimated the parameter K using numerical simulation. Illman et al. [13] and Cardiff et al. [14] conducted pumping tests at a field test site and estimated the K value using the same method. These studies have demonstrated the advantages of HT through laboratory sandbox experiments, field tests, and numerical simulation. In addition, previous studies also have compared the equivalent homogenization method (an effective parameter/macrodispersion approach) [15] and kriging method (a heterogeneous approach using ordinary kriging based on core samples) [16] with HT in terms of the ability to characterize the heterogeneity of experimental sandboxes, and found that HT was most accurate one in predicting K and hydraulic head, but it was only comparable with kriging with slight improvement [15]. By incorporating multilevel a priori geological information into the sandbox configurations of previous studies, this study aims to explore the role of HT in the improvement of the accuracy of the K field prediction, as well as the influence of aquifer heterogeneity, with higher accuracy for the solute transport. Particularly, different from the previous studies, we focus on the smaller scales of each layer of heterogeneity. The results of this study would provide a better understanding of the fine characterization of aquifer heterogeneity for the prediction of solute transport.

2. Materials and Methods

2.1. Numerical Sandbox Construction

In reference to Zhao [17], the sandbox is divided into 18 layers, and the working areas of the corresponding layers are drawn using VSAFT2 (Figure 1), which is two-dimensional in the vertical cross-section. There is a total of 48 wells, and their positions are determined according to Liu [12]. VSAFT2 is used to mesh the constructed sandbox area into finite elements, and the entire model consists of 73 rows and 163 columns of 1-cm square unit elements. The thickness of the sandbox is set as 10.2 cm, and the more detailed information about the sandbox is introduced in Zhao [17], including sand types and initial K of 18 layers. When in a forward simulation of VSAFT2, the first step is that the homogeneous initial K is assigned for the entire working area (here, we focus on the smaller scales, which means the initial K in a specific layer is homogeneous, but varies for different layers, as shown in Figure 1). Then random seeds, variance, and correlation scale are all set to the same under various conditions. Among them, the uncertainty of variance is slightly higher and, in addition, the uncertainty set by different experimenters can vary noticeably. Next, a reference K field with 18-layered randomness is produced in the middle process. Finally, the synthetic reference H (head) field in 8 pumping tests is obtained by the VSAFT2, and it should be noted that, for every pumping test, we could obtain a group of head data at 46 observation ports (excluding the injection and pumping wells), and finally there are 8 groups of observations of head data that are used in the kriging and hydraulic tomography. In contrast, in an inverse simulation, the spatial distribution of K field is derived from a given H field, and here we discuss the best one (note that initial S is 3.94 x 10⁻⁵ cm⁻¹ for each layer, but we do not focus on S field in this study). The variance in this study is set to 0.01 cm²/s² because the hydraulic conductivity of the layers in the sandbox itself is approximately 0.1 cm/s [17]. Considering the a priori geological information, the correlation scale at a small level is used as the correlation scale for the entire area. That is, the correlation lengths in the horizontal and vertical direction are set at 20 cm and 8 cm, respectively. As the inner heterogeneity of hydraulic conductivity of each layer is studied, there exists prior information for 18 layers compared with the unique K in each layer. In other words, we could know finer characterization for the sandbox at the first instance so that we call it “multilevel prior information”.

2.2. Pumping Tests

In the forward model, a zero-flow boundary condition is set at the periphery. Because the model is at a vertical section to exclude other factors, such as the effects of model errors and gravity [18,19], all initial heads are assumed as the total head (i.e., 100 cm). A total of 8 pumping test conditions are designed, and each group contains 46 observation heads, which are regarded as the reference heads. The locations of pumping/injection wells in each group of pumping tests are selected according to Illman [16]. Well spacing optimization is outside of the scope of this study [20,21,22]. The rate is 0.48 cm³/s for both pumping and injection for all 8 pumping tests, and it is set to be negative for pumping and positive for injection.

It is worth noting that head response data are obtained from the pumping test. However, because signals of head are easily transmitted through aquifers, it is not sufficient to only examine the degree of hydraulic head fitting; solute transport must also be predicted. During the simulation of tracer solute transport in the present study, a solute at a concentration of 1513 mg/mL was first injected for 10 s, followed by continuous injection of water until 6000 s, and finally the solute transport was quantitatively analyzed.

2.3. Methods

In this study, not only does the entire sandbox area have a large-scale heterogeneity, but each layer also has a smaller-scale heterogeneity of its own. As a priori geological information is assigned to the 18 layers, we first discuss the equivalent homogenization method. The corresponding homogeneity model is different from the previous theoretical homogeneity model in that each layer is treated as homogeneous (i.e., the K of each layer is a constant) [16] in the test, and the initial value of K is determined according to the hydraulic conductivity of different sand layers in Appendix A of Zhao [17] to examine the effect of the number of pumping cycles on the accuracy [23].

The hydrogeological parameters of aquifers (e.g., K) are random variables changing with the spatial position, which can be analyzed by the geological statistical methods. The core of geostatistics is “Kriging”, which is also called Kriging spatial interposition. Kriging is a stochastic linear estimation method based on the observed values of variables at different positions or some known parameters, as well as the a priori information on variable spatial statistical structures [24,25]. It obtains the optimal unbiased estimation of geological variables and quantitatively analyzes the uncertainty of the corresponding estimates. The value at an unknown point can be estimated by the weighted sum of known observation values. However, the Kriging is not suitable for the limited core samples, as well as the smaller scales.

HT estimation is a method of using the successive linear estimate (SLE) to estimate the spatial distribution of hydrogeological parameters [26]. It is an extension to the spatial random linear estimation (kriging or cokriging) [9]. In essence, HT collects nonredundant information through a limited number of wells, and some previously uncollected information could be added during each pumping test [27,28,29]. Assimilation of such data series with HT can provide a more accurate description of the aquifer. When estimating the hydraulic conductivity K based on observed head response data yielded from pumping tests, the SLE considers the nonlinear relationship between hydrogeological parameters (such as K) and system responses (such as H) [30,31]. Its core algorithm is as follows [32]:

$$\begin{array}{l}{\widehat{f}}^{(r+1)}(x_m)={\widehat{f}}^{(r)}(x_m)+{\displaystyle \sum _{j=1}^{n_h}{\omega }_{mj}^{(r)}}\left({\varphi }_j^{\ast }(x_j)-{\varphi }_j^{(r)}(x_j)\right)\\ ={\widehat{f}}^{(r)}(x_m)+{\displaystyle \sum _{j=1}^{n_h}{\omega }_{mj}^{(r)}}\left[\left(h_j^{\ast }(x_j)+H(x_j)\right)-\left(h_j^{(r)}(x_j)+H(x_j)\right)\right]\\ ={\widehat{f}}^{(r)}(x_m)+{\displaystyle \sum _{j=1}^{n_h}{\omega }_{mj}^{(r)}}\left(h_j^{\ast }(x_j)-h_j^{(r)}(x_j)\right)\end{array}$$

(1)

where ${\widehat{f}}^{\left(r\right)}$ is the estimated value at $x_m$ from the r^th iteration;${\varphi }_j^{\ast }\left(x_j\right)$ is the observed head at $x_j$; ${\varphi }_j^r\left(x_j\right)$ is the head value calculated using the current ${\widehat{f}}^{\left(r\right)}\left(x\right)$; $H\left(x_j\right)$ is the average head calculated using the unconditional mean K value with no added information; $h_j^{\ast }\left(x_j\right)$ and $h_j^r\left(x_j\right)$ are the disturbance terms for observed and calculated heads, respectively; and ${\omega }_{mj}^{\left(r\right)}$ is the weight of the difference between the observed and estimated values at point $x_j$ for the estimated f at point $x_m$ in the r^th iteration.

In the solute transport analysis method, the displacement of solute particles is determined using the concept of spatial moments. The zeroth moment $M_0$ represents the total mass, and the ratio of the first to the zeroth moment, ${\mu }_x$, represents the mean position the solute reaches at this time. Solute transport can be quantitatively evaluated according to the conservation of mass concentration and the deviation of displacement at the central point. In this study, the total mass concentrations (mg/mL) of the reference field and the estimated field are estimated by using the zeroth moment, i.e., by multiplying the concentration at each point by time and then calculating the integral. Note that the effective porosity was set to 0.36 for all cases. The computed value of the longitudinal dispersivity is 52.08 cm. The transverse dispersivity values are zero.

3. Results and Discussion

3.1. Equivalent Homogeneous Conceptual Model

The equivalent homogenization method only considers the fitting of the head and can be used to evaluate the effect of the number of pumping cycles on the accuracy. Figure 2a–h show the fitting between the observed and estimated head values of 46 wells as the number of pumping cycles increases from one to eight. These figures indicate that the scatter points are more concentrated on the 45° line as the number of pumping cycles increases and that the correlation coefficient R² gradually increases. However, R² has reached 0.99 upon the completion of the 6th pumping test, after which there is no notable increase in the accuracy. Therefore, when we determine the best accuracy, we do not need more pumping tests, which can save more energy and is less time-consuming.

3.2. Heterogeneity Characterization by Kriging

Linear estimation was performed by kriging using the reference K values (hard data) corresponding to the 48 ports. We assume that the core sample data can be directly extracted by VSAFT2. Two scenarios were conducted, one without and the other with multilevel a priori geological information. The results are shown in Figure 3. Figure 3a shows the distribution of the reference K field, and Figure 3b,c represent the distribution of the K field estimated by the two scenarios. A comparison of the two scenarios reveals that the inverted K field provides clearer hierarchical information and better characterization of heterogeneity, indicating that incorporation of multilevel, small-scale a priori geological information can improve the accuracy of the estimated K field. It should be noted that because kriging is an estimation method based on interpolation and extrapolation, the resulting estimated K field has a distribution with a relatively smooth hierarchical outline [33,34].

In addition, it can be observed from the scatter plots that the correlation coefficient R² is 0.25 and the slope is 0.32 in Figure 3e and R² is 0.21 and the slope is 0.29 in Figure 3d. The scatter points in Figure 3e are more concentrated and the blue line is closer to the 45° line, suggesting that kriging method which considers smaller scales has a better overall fitting. In summary, the incorporation of multilevel a priori geological information enables more accurate characterization of heterogeneity, indirectly proving the importance of hierarchical information.

3.3. Heterogeneity Characterization by HT

The estimated value of K for the entire work area is inverted with the VSAFT2 model by inputting the 8 groups of observed head data under forward modeling into the HT model. Figure 4a,b show the spatial distributions of the reference K field and the HT-estimated K field, respectively. It can be seen that the interlayer information is more significant and the hierarchical outline is more clearly characterized; the spatial distribution trends of K values agree very well, with one-to-one correspondence in the high-K and low-K regions. Moreover, the slope of the scatter plot in Figure 4c reaches 0.74, which is very close to the 45° line, and the correlation coefficient R² reaches 0.58 (compared with the maximum correlation coefficient of 0.25 by kriging method), indicating that the HT-estimated K field is quantitatively closer to the reference field. Therefore, HT remarkably improves the simulation of the K field and hydraulic head compared with the kriging estimation.

3.4. Tracer Test

3.4.1. Comparison of Concentration Contours in the Reference Field and the Estimated Field

Figure 5 shows the positions of solute transport simulated at output time points of 10 s, 70 s, 270 s, 1000 s, 2000 s, and 6000 s under the reference K field in a group of pumping/injection tests (by pumping at well #7 and injecting solute–water mixture at well #42). Figure 6 shows the positions of simulated solute transport under the optimal estimated K field. In the figure, the concentrations have been normalized so that a darker color represents a higher concentration. A comparison of Figure 5 and Figure 6 reveals that the positions corresponding to each time point are approximately the same and the transport direction and the diffusion range of the tracer coincide well, indicating that the HT-derived optimal estimated K field can sufficiently predict the transport direction and concentration of solute.

3.4.2. Mass Conservation of Solute Concentration

Table 1. Total mass concentration at each time point.

Time Point (s)	Concentration of Reference Field (mg/mL) Multiplied by Time (s)	Concentration of Estimated Field (mg/mL) Multiplied by Time (s)
10	12,617.05	12,348.89
70	13,389.39	13,184.79
100	13,530.88	13,437.63
270	13,965.98	14,346.23
500	14,235.86	14,513.34
700	14,176.88	14,257.88
1000	14,047.74	14,160.73
2000	14,030.94	14,016.92
3000	13,845.20	13,828.36
5000	12,698.78	12,609.27
6000	11,431.97	11,283.76

presents the total mass of the reference and estimated fields at different time points. The correlation coefficient between the total mass concentrations of the reference field and the estimated field is greater than 0.96, indicating that the total mass concentration was properly estimated. It should be noted that before the solute was pumped out of well #7, the theoretically calculated total mass concentration of the solute should have been equal to the total mass concentration of the injected solute (1513 mg/mL) multiplied by the injection time (10 s). However, the estimated concentrations shown in Table 1 are not equal to the theoretical concentrations, and such differences might be attributed to computer simulation errors.

3.4.3. Route of Solute Mean Position

The center positions of the tracer at different time points are estimated using the ratio of the first moment to zeroth moment, ${\mu }_x$ The direction of solution transport is recognized by determining the displacement of the mean position, and the position coordinates of the tracer at different time points are extracted and plotted on a plane, as shown in Figure 7. The curve and direction of solute displacement under the estimated field basically coincide with those under the reference field, respectively, indicating that the simulated K field can properly characterize the heterogeneity of the sandbox used in our experiment.

Table 2. Displacement and velocity of tracer center at different time points.

Time (s)	Displacement of Reference Field (cm)	Velocity of Reference Field (cm/s)	Displacement of Estimated Field (cm)	Velocity of Reference Field (cm/s)
10	0.23	0.0227	0.48	0.0478
70	2.35	0.0310	3.38	0.0385
100	3.14	0.0263	4.25	0.0289
270	6.36	0.0175	7.75	0.0191
500	9.27	0.0118	10.58	0.0113
700	11.28	0.0101	12.53	0.0097
1000	14.03	0.0089	15.37	0.0094
2000	23.48	0.0094	25.34	0.0099
3000	34.48	0.0110	36.86	0.0115
5000	51.35	0.0070	53.18	0.0064
6000	56.03	0.0047	57.00	0.0038

shows the approximately estimated displacements of the tracer center after different periods of time and the transport velocity in each period of time. It is indicated that the transport velocity of the solute gradually decreases with the passage of time, which may be due to the decrease in the concentration gradient, the heterogeneity of the aquifer, and the radial flow in the vicinity of the wells. According to the displacement data, it can be roughly estimated that the tracer was initially pumped out at a time point between 1000 s and 2000 s. The displacement data are accumulated starting at 10 s, and it is found that the total displacement at 1000 s and 2000 s is approximately 84.95 cm and 110.29 cm. In addition, as shown in Figure 5 and Figure 6, the coordinate of the outer boundary is (x,y) = (110 cm, 25 cm) at 10s, (40 cm, 50 cm) at 1000s, and (12 cm, 70 cm) at 2000s. However, the distance between the well #7 and the outer boundary of solute plume (at 10s) is almost exactly 100 cm. That is, the time when the solute moves to the mouth of well #7 should be between 1000 s and 2000 s.

Figure 8a,b show the concentration range of the estimated field at 1000 s and 2000 s, respectively. As seen, the solute (tracer) has noticeably not arrived at well #7 at 1000 s, but the lightest blue area has already passed over well #7 at 2000 s, consistent with the conclusion of the displacement data analysis. Considering the previous analysis, it is effectively proved that HT is able to not only fit head signals very well but also accurately predict solute transport, and it also offers a good alternative for the treatment of pollutants.

4. Conclusions

Based on a synthetic sandbox incorporated with multilevel a priori geological information, we compared the roles of kriging and HT in the fine characterization of aquifer heterogeneity. The solute transport under the optimal K field is examined, and the effects of the number of pumping cycles on the accuracy of aquifer heterogeneity are investigated in this study. The following conclusions can be drawn:

• As the number of pumping cycles increases, the accuracy of the inversely obtained K field is improved, although this improvement is no longer significant after the 6th pumping test. For specific regional hydrogeological studies, we can use a sufficient number of pumping cycles during a field experiment and determine the number of pumping cycles beyond which the effect of this number is no longer significant, thus providing a reference for determining the number of pumping cycles needed for different hydrogeological points of the same region.
• Incorporation of multilevel a priori geological information improves the accuracy in the estimated K field. When small-level heterogeneity is considered, the accuracy in the characterization of the K field heterogeneity of the aquifer inversely obtained by HT is generally higher than that by kriging.
• The HT-derived optimal estimated K field is able to not only fit head signals very well but also more accurately predict solute transport, suggesting that fine characterization of aquifer heterogeneity is of great importance to the prediction of solute transport.

The sandbox used in this study is directly synthesized on a computer using VSAFT2, lacking experimental data for practical validation and containing systematic errors in the numerical model. And it should be noted that the study is 2D, which is not realistic, especially in the typical, local scenario of HT. In addition, in this study, only the analysis of spatial moments is performed in the prediction of solute transport. Future work can combine the breakthrough curves of other observation wells to investigate the accuracy in the prediction of solute transport using time moments. In addition, we will attempt to combine the time moment approach with the sandbox test data to verify the importance of fine characterization of aquifer heterogeneity in actual applications.