Estimation Parameters of Soil Solute Transport Processes by Using the Electric Resistivity Method

Abstract

Preferential solute transport is a common phenomenon in soil, and it is of great significance to accurately describe the mechanism of pollutant transport and water and soil environmental governance. However, the description of preferential solutes still relies on applying solute breakthrough curves for model parameters fitting. At present, most of the solute breakthrough curves are obtained indoors, and with some limitations. Therefore, this study established a method for securing solute breakthrough curves based on the electrical resistivity method. The research results show that the change in soil concentration during the tracer infiltration process can be captured by establishing the fitting relationship between soil resistivity and solute concentration. Then the solute breakthrough curve can be found. Through a time moment analysis, the difference between the breakthrough curve parameters obtained by the traditional method and the resistivity method is slight; the average error is less than 10%. On this basis, the sensitive response of the parameters of the “mobile–immobile” model to concentration was elucidated through different concentration tracer experiments, among which β and D are more sensitive, and w is less sensitive. The suitable tracer concentration range should be 50–120 mg/L. Therefore, the established method could obtain the breakthrough curves and describe the transport of preferential solutes at the field scale.

1. Instruction

Preferential flow is a common phenomenon of water and solute transport in soil, and the pollutants it carries can reach deep soil and aquifers [1,2,3]. Due to the short interaction time between preferential water and solutes and the soil matrix, soil pollutants cannot be wholly absorbed and degraded by the soil matrix, thus increasing the risk of polluting groundwater and affecting biological growth and human health [4,5]. Therefore, the transport mechanism of solutes in soil with preferential channels is key to comprehensively revealing the soil solute pollution transportation process and then providing schemes to protect soil and water security [6,7].

Many scholars [8,9,10,11] have detailed research on the mass flow, diffusion, and chemical coupling in soil solute transport. The soil convection–diffusion equation (Convection Dispersion Equation) proposed by Nielson et al. [12,13,14,15,16] is widely used in the numerical simulation of solute transport. The equation could give a good description of convection and dispersion in the process of solute transport under homogeneous conditions. However, due to its limited interpretation of the early breakthrough and trailing phenomena in breakthrough curves (BTCs), the mobile–immobile water model (MIM) can improve the effectiveness of the CDE equation in the physical and chemical non-equilibrium description of the solute transport process [17,18,19].

However, both the CDE and the MIM need to use the breakthrough curves for parameter estimation. The breakthrough curves are generally measured by an indoor saturated soil column test, which will disturb the soil, thereby affecting the accuracy and limiting the acquisition of field-scale parameters [20,21,22,23]. Therefore, the current methods still have gaps in obtaining regional-scale soil solute transport parameters. At present, it is urgent to break through the limitations of soil solute parameter acquisition, to solve the problem that the existing preferential flow model is only suitable for indoor soil column research and lay the foundation for the regional-scale preferential flow solute transport simulation. However, high-density electrical resistivity tomography (ERT) is a non-destructive detection method that has recently become widely used in acquiring soil moisture and solute transport parameters due to its flexible detection scale and suitability for three-dimensional dynamic monitoring [24,25,26]. Many scholars [27,28,29] have used the method to study water and solute transport parameters in soil and achieved fruitful results, laying a theoretical foundation for acquiring soil-water parameters. However, there is still insufficient data to support the use of this method to obtain field BTCs.

Consequently, this study intends to use ERT to get unsaturated solute breakthrough curves with a specially designed lysimeter. The purpose is to carry out a tentative analysis of the solute breakthrough curves in the field. This will provide a new method for the acquisition of soil solute transport parameters at the regional scale, thus making it possible to accurately describe soil solute transport. Since it is a field method for obtaining soil solute transport parameters, and is not affected by soil physical properties, it can be applied to various soil types worldwide. At the same time, in order to evaluate the accuracy of the obtained parameters, the characteristics of the MIM’s parameters are intended to be studied to provide data support for the description of the field-scale solute transport in the preferential flow.

2. Experimental Theory and Methods

2.1. Device Design and Experimental Scheme

The undisturbed soil of this experiment was collected in the Changhua River Basin, which is located at the source of the Qiantang River in Zhejiang Province. The primary soil type was paddy soil, the average proportion of silt, clay, and sand were 45.35%, 38.61%, and 16.04 respectively, and the preferential flow was developed [30]. The basin belongs to the subtropical monsoon climate zone, with four distinct seasons and uniform distribution of precipitation seasons. The annual precipitation is between 800 m and 1600 m. The yearly average temperature in this area is 15.3 degrees, the average yearly evaporation is 1150 mm, and the average annual rainfall is 1703.01 mm. According to the investigation, the thickness of the vadose zone is about 100 cm. Therefore, this study selected the area with developed channels (0–100 cm) as the primary research area [31].

The experiment was completed in a self-designed special lysimeter. The design of the lysimeter is shown in Figure 1. It was a cube with a length, width, and height of 1 m. It was reinforced with steel rings, and stainless steel was installed on the four sides. Sixty-four electrodes were installed on each surface, with 256 electrodes for the entire cube. All electrodes were connected to an E60DN resistivity meter (Geopen, China) using copper wires to measure soil electric resistivity. The top of the cube soil column related to a constant water head device as the upper boundary for infiltration under different water head conditions. To minimize the disturbance to the undisturbed soil, the lysimeter was directly inserted into the undisturbed soil in the field. The surrounding area was excavated. Then, the bottom plate of the lysimeter was inserted to maintain the integrity of the soil core.

To obtain the breakthrough curve of the soil and its spatial variation characteristics, the soil core was conceptually divided into four parts in this study. Each piece was installed with soil moisture sensors (TDR100, CAMPBELL) at 20 cm, 40 cm, 60 cm, and 80 cm from the soil surface to the bottom for collecting the initial soil moisture. In addition, four small vacuum pumps were installed at the bottom of each part, and the lower boundary condition with a certain suction was set during the experiment. Four vacuum pumps controlled four areas respectively; the vacuum pump was connected to the collection bottles, then the collected solution volume, concentration, and other parameters were determined in the laboratory.

Soil samples were collected from the corresponding regions and brought back to the laboratory to obtain the soil’s physical parameters. The soil texture information is shown in

Table 1. The measured soil physical parameters for 4 sites.

Sites	Clay (%)	Silt (%)	Sand (%)
#1	40.21	43.19	16.6
#2	28.51	50.72	20.77
#3	38.47	47.64	13.89
#4	47.25	39.85	12.9

and was measured by laser particle sizer (S3500, Manufactured by Microtrac MRB). The soil texture differences in the four regions are apparent. The subregional parameter acquisition method can effectively capture the differences caused by soil heterogeneity, thereby laying a foundation for further description of soil preferential solute transport. The overall scheme of this experiment is shown in Figure 2, which is generally divided into two parts. The first part is establishing a BTCs acquisition scheme based on the electric resistivity method. During the solution infiltration test, the concentration and resistivity of the collected soil samples were measured and a relationship established to obtain BTCs. The other part is to collect solution samples through a specially designed vacuum pump device and use traditional methods to get BTCs. The primary purpose is to verify the feasibility of the established BTCs acquisition scheme based on the electrical resistivity method by comparing it with the traditional scheme. Consequently, the uncertainty analysis of the proposed project is carried out. The parameters of MIM are determined in terms of obtained BTCs, which provide theoretical support for the field application of the electric resistivity method to get the solute transport parameters of the vadose soil.

2.2. High-Density Electric Resistivity Tomography Method

As shown in Figure 1, the horizontal and vertical electrode distances were 14 cm; the electrodes were divided into eight columns on each face, 32 columns were engaged with the resistivity measurement using the dipole-dipole array; 62,496 measurements were conducted in total. To reduce the uncertainty of the measurement, the backward error analysis method was applied to estimate the error. By calculating the difference between the normal and back measurements, errors more significant than 15% were excluded from further calculation. Reasonable data were further calculated using the inversion method proposed by Lu et al. [32,33] to obtain the soil resistivity value. However, the process will not be discussed here due to space limitations. In addition, the resistivity of the soil samples was determined by the four-electrode method in the laboratory, and the chloride ion concentration was determined by the national standard method.

2.3. MIM Parameters Determination Using BTCs

The traditional convection–dispersion model cannot accurately describe the preferential flow. According to previous research, the “mobile–immobile model” (MIM) has certain advantages for describing soil solute transportation for preferential flow, which divides the soil into a moving water area (A) and a non-moving water area (B). There is no water exchange but solute diffusion between the two areas, and it only depends on the concentration difference between the two areas. Therefore, MIM was used to fit the obtained breakthrough curves. The dimensionless form of MIM can be described as [34]:

$$\beta R\frac{\partial C_A}{\partial T}=\frac{1}{P}P\frac{{\partial }^2{C}_A}{\partial Z^2}-\frac{\partial C_A}{\partial Z}-\omega \left(C_A-C_B\right)$$

(1)

$$\left(1-\beta \right)R\frac{\partial C_B}{\partial T}=\omega \left(C_A-C_{B }\right)$$

(2)

$$\beta =\frac{\theta +f{\rho }_b{k}_d}{\theta +{\rho }_b{k}_d}$$

(3)

$$\omega =\frac{\alpha \left(1-\beta \right)RL}{v}$$

(4)

$$R=1+({\rho }_b{k}_d/\theta )$$

(5)

$$P=v_{m}L/D_m$$

(6)

The initial condition is the given concentration boundary, and the upper boundary is:

$$C_A\left(0,T\right)-\frac{1}{P}\frac{\partial C_A\left(0,T\right)}{\partial Z}=C_{0}T$$

(7)

The lower boundary is the fixed water head boundary in the above equations: C, Z, T are dimensionless concentration, depth, and time respectively; P is the peclet number; f is the proportion of the mobile water adsorption area; β is the proportion of the mobile water area; α is the mass exchange coefficient between the two areas; L is the length of the soil column; C_A, C_B are the dimensionless concentration values of the mobile water area and the immobile water area, respectively; k_d is the adsorption coefficient; ${\rho }_b$ is the soil bulk density; D_m is the moving area dispersion coefficient; v_m is the pore flow rate in the moving area; R is the resistance factor; v is the pore water flow rate; θ is the volumetric water content. There are D, P, R, β, and ω in the model that need to be fitted by BTCs. In this study, Levenberg–Marquart’s algorithm [35], widely used in nonlinear fitting, was used to fit the measured BTCs.

3. Results

3.1. Establishment of R-cl⁻ Relationship

Determining the quantitative relationship between resistivity and tracer chloride ions is crucial for obtaining BTCs by the electrical resistivity method. In this study, samples were taken from four areas in the undisturbed soil core. The soil resistivity under different ion concentrations is measured in the laboratory, and the relationship was determined by the nonlinear fitting method. However, the initial water content of the soil has a significant impact on the relationship, so, referring to the actual initial water content of the soil, the resistivity and ion concentration under different initial water content was calibrated.

The initial soil content measured shows a small change. Figure 3 shows the calibrated relationship between resistivity and chloride ion concentration fitted by the power and exponential functions, respectively. According to the fitting results of the judgment parameter R² and the Nash efficiency coefficient in

Table 2. Resistivity–chloride fitting results and error assessment parameters.

Method	Result	R2	Nash Coefficient
Power function	$R=159.7C^{-0.3683}$	0.9353	0.95
Exponential function	$R=53.18e^{-0.005369C}$	0.7679	0.82

, the power function is better than the exponential function. From the overall trend, the resistivity decreases with the increase of chloride ion concentration. This is because the electrical conductivity increases with the growth of chloride ions, which makes the resistivity tend to fall. In addition, due to the nonlinear relationship between resistivity and chloride ion concentration, the error distribution range based on the measured and simulated data was calculated. Figure 4 shows that when the chloride ion concentration is between 50–160 mg/L, the error between the measured value and the simulated value is small, less than 10%; when the ion concentration is less than 50 mg/L, the error is between 10–20%. However, when the ion concentration is greater than 160 mg/L, the relative error is significant, ranging from 20% to 60%. This is due to the error caused by the reduced sensitivity of the resistivity to high concentrations of chloride ions.

3.2. The Determination of the BTCs

In this study, chloride ions with an initial concentration of 70 mg/L were used for infiltration experiments. According to the determined resistivity–chloride ion concentration relationship, the resistivity of the four regions of the soil column during the infiltration process was monitored, and the inversion calculation was carried out. The obtained resistivity value of each area at different times was converted to chloride ion concentration to estimate the breakthrough curve. Figure 5 shows the obtained BTCs for four regions determined by the ERT method. As shown, the MIM has a good fitting result for the breakthrough experimental data in which the correlation coefficients are greater than 0.9 while all R² are above 0.8 in each site. Hence, the presented method could capture the front, peak, and tail features of the breakthrough process. Overall, the four regions show similar characteristics: single peak, asymmetry, and tailing. However, there are still some differences, such as the arrival time of peak concentration. The difference of peak concentration arrival time in sites 1, 3 and 4 is small, and site 2 reaches the peak concentration first. Thus, it can be inferred that there are preferential channels in the four regions, and certain heterogeneity happens in all sites. In addition, the BTCs concentration value obtained by the resistivity method is generally low, which may be due to the systematic nature of the resistivity measurement errors.

3.3. Parameter Characteristic Analysis of Breakthrough Curve

This study investigated the effect of different tracer concentrations on the fitting parameters of breakthrough curves. Figure 6 shows the dispersion parameters β, D, and w fitted by ERT and traditional methods at three different concentrations (concentration 1: 15 mg/L; concentration 2: 85 mg/L; concentration 3: 140 mg/L) respectively. The results show that parameter β is more sensitive to concentration, as the error of the low and medium concentration is relatively small, the R² is a respective 0.89, 0.98, the RMSE is a respective 0.12, 0.11. The results also show that the data scatter degree is larger at the higher concentration, the R² being 0.87 and the RMSE 0.13, where the parameter deviations obtained by the two methods were the smallest under the condition of medium concentration of tracer. This is caused by the difference in the sensitivity distribution of β to the concentration. The sensitivity of β is higher at medium and low concentrations, and the fitting value is relatively accurate.

In contrast, the sensitivity to high concentration solutions is reduced, thus affecting the fitting accuracy. The dispersion D also showed similar variation characteristics, the most sensitive to concentration, with the most significant error at high concentration, with R² of 0.86 and RMSE of 3.2. The low and medium concentration error was small, the medium concentration error was the smallest, and R² was 0.92 with an RMSE of 1.72. The RMSE value of w at high concentration was 3.5, but the R2 values at three concentrations were more significant than 0.9, and the difference was slight. From the response results of the three parameters to the concentration, β and D were more sensitive to the change of concentration, and the sensitivity of medium and low concentrations was relatively high. w was less sensitive to changes in concentration, and the parameter fitting accuracy gap under different concentrations was small. According to the method proposed by [36], the sensitivity of parameters at different concentrations was calculated in this study, and the calculation results are shown in Figure 7. This result is also consistent with the results in Figure 6. D was the most sensitive to concentration, and the sensitivity value was the highest when the concentration was about 55 mg/L and 95 mg/L; β is the second; the most heightened sensitivity happens at concentrations around 35 mg/L and 70 mg/L; while the sensitivity of w to concentration ws poor which shows a slight difference with the change of concentration. This provides data support for selecting a reasonable concentration of tracer in the breakthrough curve acquisition experiment.

4. Discussion

Time moment analysis [37] is usually used to quantitatively describe the shape of BTCs and estimate solute transport parameters. This study used time moments to assess the gap between traditional and new methods. The time moment can be described as:

$$M_p={{\displaystyle \int }}_0^{\infty }{T}^{p}c\left(z,T\right)/c_{0}dT, \mathrm{P}=0, 1, 2,\dots (8)$$

(8)

where Z is a dimensionless spatial coordinate; T is dimensionless time; C and C₀ are solute concentration and input solute concentration respectively (M/L³); P is the order of the moment. Thus, the first-order standard moment $u_1^{\prime }$ represents the breakthrough time, the second-order center moment $u_2$ represents the average extension, while S represents the degree of symmetry of the breakthrough curve. In this study, time moments were used to calculate BTCs obtained by traditional methods [21] and ERT respectively, so as to obtain $u_1^{\prime }$, $u_2$, S. The soil breakthrough curves of different sites, obtained by using the traditional method and ERT, are shown in Figure 8.

Table 3. Calculation results of time moments for BTCs obtained by different methods.

Sites		${\mathit{u}}_1^{\prime }$	${\mathit{u}}_2$	S
#1	ERT	2.52	1.37	1.18
	Traditional way	2.27	1.53	1.21
#2	ERT	1.88	2.97	1.86
	Traditional way	1.96	2.54	1.97
#3	ERT	2.98	1.32	1.25
	Traditional way	3.62	1.17	1.13
#4	ERT	1.51	3.83	1.67
	Traditional way	1.43	3.25	1.72

presents the calculation results of the time moment of breakthrough curves for two different methods. From the value of S, which represents the degree of symmetry in Table 3, the variability of soil pores in No. 2 and No. 4 sites is higher than that in No. 1 and No. 3 sites.

It can be seen that the relative concentrations obtained in this series of experiments were all less than 1, which is more consistent with the results obtained by other scholars [38], mainly due to the presence of water in the soil, which has a specific dilution effect on the solution. Meanwhile, the measurement sensitivity is another contributor to the occurrence of low concentrations. Binley et al. [39] have pointed out that the low sensitivity of ERT away from the electrode led to the underestimation of the change of water, which they believed was the main reason for the mass balance error in the tracer experiment. However, in our experimental device design, we refer to the scheme of koestel et al. [40], so that the concentration deviation can be reduced to an acceptable range.

Compared with the traditional method, the first-order standard moment $u_1^{\prime }$ had an average error of 9.25%, a maximum error of 17%, and a minimum of 5%. The average error of the second-order central moment $u_2$ was 14.35%, the maximum error is 17.8%, and the minimum error was 9%; The average mistake of S was 5.4%, the maximum error was 10.6%, and the minimum error was 2.4%. The results show that the average errors were all less than 10%. The error of the coefficient S was the smallest, and the error of the second-order central moment was the largest, but the overall difference is negligible. According to previous studies [31], the more developed the macropores, the smaller the value of $u_1^{\prime }$. Therefore, it can be inferred that the macropores in the No. 4 area were the most developed, and the breakthrough time was the shortest; the macropores in the No. 1 area were relatively less developed, and the breakthrough time was longer. $u_2$ represents the average extension of the curve; the more significant the value, the more advanced the macropores. The calculation results show that the value of $u_2$ in the No. 4 area was the largest. The preferential flow had a significant effect, which is highly consistent with the calculated value of $u_1^{\prime }$.

In our experiment, it is found that the outflow concentration does not show a systematic change with the change of water supply head. This is similar to the conclusion obtained by Wehrer and Slater, [38]. Gomez et al. [41] has also pointed out that with the increase of water head, the development of finger flow may be increased, and the finger flow mainly moves laterally, which weakens the longitudinal flow and makes the resulting flow concentration show a decreasing trend. In addition, the relationship between resistivity and ion concentration established in this study adopts different fitting methods, so it is concluded that the fitting effect of power function is better, which is also confirmed by the research of Sheng et al. [42]. From the similar research carried out by relevant scholars [43,44,45], the fitting relationship between the two does not consider the physical mechanism in detail, which has a certain empirical nature. Therefore, more tests need to be carried out for different types of soil (water content, texture, and structure) in order to attempt different fitting schemes to establish a more accurate relationship between the two variables.

5. Conclusions

In this study, by establishing a fitting relationship between resistivity and tracer concentration, the concentration changes of soil tracer were captured, and the solute breakthrough curve was obtained. Results have shown that the relationship between resistivity and solute concentration is sensitive to concentration changes, and the appropriate concentration range should be 50–120 mg/L. The concentration fitting error outside this range is significant. Compared with the relevant parameters of the time moment, the solute breakthrough curve obtained by the ERT method has a small gap from the traditional method. There are differences in the sensitivity of the parameters of the “mobile–immobile model” to the concentration of the tracer, among which β and D are more sensitive. In contrast, the sensitivity of w is weaker, mainly caused by the difference in the sensitivity of different parameters to the concentration. A reasonable concentration range should be paid attention to.