Experimental and Numerical Simulation Study of Flow and Solute Transport in Pore-Fractured Media Based on High-Density Resistivity Method

Xiaosan Yan; Xiaoyong Wang; Lei Ma; Can Ge; Jiazhong Qian; Songjian Ao

doi:10.2113/2021/7010329

Abstract

Pore-fractured media is a ubiquitous phenomenon throughout the world and is a high degree of heterogeneity. The mechanism of water flow and solute transport in the media is still not fully clear. In this study, experiments were conducted in pore-fracture media with various hydraulic gradients to establish the relationship between the hydraulic gradient and specific discharge by using a side guide and to test the resistivity of tracer (sodium chloride) by using a high-density resistivity monitoring system. Quantifying models were, respectively, set up, and results were compared with the traditional methods. Main results are obtained that (1) the phenomenon from Darcian flow to non-Darcian flow was found in pore-fractured media. Both the Forchheimer and Izbash equations can well describe the flow process in the back part of the curve; (2) the phenomenon from non-Fickian to Fickian then to non-Fickian was found in fractured media; (3) good performance of CTRW-TPL has been obtained for both the larger $R^{2}$ and smaller RMSE values, their counterparts resulting from the ADE model.

1. Introduction

With rapid development in industrial and agricultural production, environmental pollution, especially groundwater pollution, becomes increasingly more serious. Some pollutants enter the fracture media through the spread of solutes in porous media. Therefore, controlling contamination migration in such a fracture media is very important for both pollution control and aquifer remediation.

Extensive research studies were conducted in the past to investigate flow and solute transport in porous media or fractured media. For example, a noninvasive image processing method was developed to map the spatiotemporal evolution of solute concentration in two-dimensional porous media [1]. Zaheer et al. [2] carried out a series of solute transport experiments in one-dimensional clay soil columns, to identify key parameters controlling the transport process. Chen et al. [3] focused on the dispersion process in rough single fractures under non-Darcian flow conditions. Huang et al. [4] designed an experiment to perform tracer tests in an unsaturated fractured rock with a horizontal fracture with a negligible matrix.

In the past decades, many analytical solutions were developed to model solute transport or contaminant transport processes in a single fracture-matrix system [5–8]. However, much literature related to solute transport in such fractured porous media often involved the so-called “numerical experiments” or “numerical simulations” where an imaginary rather than an actual experiment was conducted. It is a little hard to understand the flow and solute transport mechanism of fracture in fractured porous media. This is partially due to the little research about flow and solute transport in fractured-matrix media through actual experiments and partially due to the challenges of the measurement and quantitative analysis of flow and solute transport in such a fracture-matrix media. The first and widely applied analytical solution for solute transport in a simplified single fracture was put forward by Tang et al. [6]. Zou et al. [9] analyzed and discussed the impacts of the special assumptions for deriving the analytical solution in Ref. [5] through numerical modeling. The basic physical and mathematical assumptions are as follows: (1) the flow in the fracture was assumed as laminar; (2) no mechanical and thermal processes were considered. For simplicity, only advection and dispersion in the fracture, and molecular diffusion in the matrix, were considered in the study [9]. Zhu et al. [10] deal with reactive solute transport in a fracture-matrix system through analytical and numerical modeling methods. They also assume that the flow is driven by a simple uniform pressure gradient aligned with the fracture axis. Zhu and Zhan [11] developed closed-form solutions, together with semianalytical solutions to describe the solute penetration process in an asymmetric fracture-matrix system. To develop closed-form analytical solutions, the flow velocity in the fracture is assumed to be a constant. Zhou and Zhan [12] developed a new analytical model for reactive solute transport in a fracture-rock matrix system with asymmetric distribution of transport parameters of the rock matrix. The basic assumptions are as follows: (1) the groundwater flow only occurs in the fracture from left to right and has a constant velocity; (2) groundwater in the upper and lower rock matrixes is immobile. Chen and Zhan [13] propose Green’s function-based solution in a real-time domain for a parallel fracture-matrix system. They supposed that the $x$ and $z$ axes are the longitudinal and transverse dispersion directions, respectively. They also supposed that steady-state saturated flow with a constant velocity occurs in fractures, but not in matrixes.

Up to now, relatively little studies have been devoted to the analysis of flow and solute transport of fractured-matrix media using actual experiments. Faulkner et al. [14] made use of a laboratory analog experiment to simulate groundwater flow and solute transport in a karst aquifer with one conduit buried in a matrix. The Darcy law is used to describe groundwater flow in the matrix, but the Stokes equation is adopted to describe flow in conduits. Chen [15] discussed flow and solute transport in a fracture with the different flow velocity, geometry, and matrix attributes. It involves matrix properties, but the contents of study are relatively simple. Li et al. [16] designed a testing system of double medium seepage characteristics which is formed of a double medium, water cycle system control module, and data collection module. The test system solves some key issues of porous media and fractured media, for example, boundary condition simulation and water exchange information. But there are still many problems to be solved.

Hence, the mechanism of flow and solute transport is still not completely understood in fractured-matrix media. How to quantify the flow and solute transport is also one of the most challenging problems.

Electrical property quantification is a useful method for better understanding of contaminant movement and monitoring in recent years [17]. Electrical resistance measurement for different geotechnical parameters of the subsurface characteristics has been reported [18–20]. Owing to accuracy and no disturbance to water and concentration fields, electrical conductivity is helpful for measurement and quantitative analysis of the solute transport behavior.

For many decades, Darcy’s law has been widely used to describe the flow in the subsurface. Also, Fickian transport is believed to be the “right” form of governing law for transport in the subsurface [21]. However, non-Darcian flow is particularly easy to occurring in heterogeneous geological formations such as in a single fracture [22, 23] or a fractured network [24, 25]. Similar to non-Darcian flow, the solute transport is found to be non-Fickian on many occasions [26–29].

The classical advection-dispersion equation (ADE) model has been widely applied to describe the tracer transport in fractures. However, many studies demonstrated the significant deviations between observed and computed breakthrough curves (BTCs) in fracture reports [30–32], which represent typical non-Fickian dynamics and cannot be efficiently described by the ADE model. There have been several attempts to seek an explanation for non-Fickian transport in fracture. These include the continuous-time random walk (CTRW), the mobile-immobile approach (MIN), and the fractional advection-dispersion equation (FADE), among others.

The objectives of this paper are to (1) experimentally disclose the phenomena from Darcian flow to non-Darcian flow in pore-fractured media, (2) set up models to describe the flow process using Forchheimer and Izbash equation, and (3) quantify the Fickian and non-Fickian behavior by comparing the CTRW-TPL model with the ADE model.

2. Theoretical Background

2.1. Model of Flow in a Single Fracture

Darcy’s law [33] for flow in a single fracture is

\begin{matrix} (1) & J = \frac{v}{K}, \end{matrix}

where

J

is the hydraulic gradient,

K

is the hydraulic conductivity (m/s), and

v

is the flow velocity (m/s). A commonly used non-Darcian flow equation is the Forchheimer law [34]:

\begin{matrix} (2) & J = - (a v + b v^{2}), \end{matrix}

where

a

and

b

are two constant coefficients.

Another commonly used non-Darcian flow equation is the Izbash equation [35]:

\begin{matrix} (3) & J = λ v^{n}, \end{matrix}

where

λ

and

n

are two constant coefficients, and

1 \leq n \leq 2

⁠.

2.2. Model of Solute Transport in a Single Fracture

2.2.1. ADE Model

The dispersive mass flux of the ADE model is proportional to the first-order spatial derivative of concentration. The one-dimensional governing equation of ADE is

\begin{matrix} (4) & \frac{\partial c}{\partial t} = D \frac{\partial^{2} c}{\partial x^{2}} - v \frac{\partial c}{\partial x}, \end{matrix}

where

c

⁠,

D

⁠, and

t

denote the concentration of the solute, dispersion coefficient, and time, respectively.

2.2.2. TPL Model

The CTRW transport equation in Laplace is usually formulated to represent the time derivative in an algebraic expression. In one dimension, the Laplace transformed concentration

\tilde{c} (x, u)

is given by

\begin{matrix} (5) & u \tilde{c} (x, u) - c_{0} (x) = - \tilde{M} (u) [υ_{ψ} \frac{\partial}{\partial x} \tilde{c} (x, u) - D_{ψ} \frac{\partial^{2}}{\partial x^{2}} \tilde{c} (x, u)], \\ \tilde{M} (u) \equiv \bar{t} u \frac{\tilde{ψ} (u)}{1 - \tilde{ψ} (u)}, \end{matrix}

where

u

is the Laplace variable, the “~” symbol represents the Laplace transformed variable,

\bar{t}

is a characteristic time, and

v_{ψ}

and

D_{ψ}

are the average transport velocity and generalized dispersion coefficient, respectively. The dispersivity can then be defined as

α_{ψ} = D_{ψ} / v_{ψ}

⁠. It is important to note that the transport velocity can be different from the mean flow velocity. Also, the generalized dispersion coefficient is distinct from that in the ADE. It is noted also that the average mass flux of the solute,

j

⁠, is defined in Laplace space through

\begin{matrix} (6) & \tilde{j} \equiv \tilde{M} v_{ψ} (\tilde{c} - \frac{α_{ψ} \partial \tilde{c}}{\partial x}) . \end{matrix}

The function $\tilde{ψ} (u)$ is the heart of the CTRW formulation and characterizes the nature of the solute movement. The function $\tilde{M} (u)$ can take on several expressions, depending on the functional form of $\tilde{ψ} (u)$ ⁠. General forms that have been described in the literature are as follows:

The truncated power-law model is

\begin{matrix} (7) & \tilde{ψ} (u) \equiv {(1 + τ_{2} {u t}_{1})}^{β} \exp (t_{1} u) Γ (- β, {τ_{2}}^{- 1} + t_{1} u) / Γ (- β, {τ_{2}}^{- 1}) 0 < β < 2, (TPL), \end{matrix}

where

Γ

is the incomplete Gamma function,

t_{1}

is the lower limit time when the power-law behavior begins, and

t_{2}

is the cut-off time describing the Fickian behavior. The parameter

β

indicates different types of anomalous transport.

2.3. Model Evaluation

This study uses the CTRW calculation programs [36] and software of CXTFIT2.1 [37] to analyze and evaluate its fitting effect of BTCs in pore-fractured media using the decision coefficient

r^{2}

and the root mean square error RMSE. These two indices are expressed by

\begin{matrix} (8) & r^{2} = 1 - \frac{\sum_{i = 1}^{N} {(C_{i o} - C_{i e})}^{2}}{\sum_{i = 1}^{N} {(C_{i o} - {\bar{C}}_{i o})}^{2}}, \\ RMSE = \sqrt{\frac{1}{N} \sum_{i = 1}^{N} {(C_{i o} - C_{i e})}^{2}}, \end{matrix}

where

C_{i 0}

and

C_{i e}

are observed and computed concentrations, respectively.

N

is the number of measured values;

{\bar{C}}_{i o}

is the mean of measured values. The bigger of determination

r^{2}

and the smaller of RMSE, the curve is better fitting.

3. Method

3.1. Experimental Setup

Figure 1 illustrates the schematic diagram of the experimental setup depicting a simplified single horizontal fracture-porous system. The physical system has two domains. The fracture domain occupies the bottom of the model and the interior measures $60 cm \times 15 cm \times 6 cm$ ⁠. The spacing between the two plates was controlled through the placement of Plexiglas strips with a height of 2.0 mm along the side edges within the fracture. The porous domain occupies the rest of the space, and the interior of porous media measures $60 cm \times 15 cm \times 19 cm$ ⁠. Both sides of the fracture domain and porous domain have an inflow and outflow reservoir. The type 304 stainless steel mesh used in the laboratory analog to separate the two domains has the following properties: mesh size 1 mm. Thirty-two electrodes with a length of 1.5 cm and a diameter of 0.2 cm were arranged on the surface of the media (in order to facilitate the discussion below, the locations of the observation point are denoted as point 1, point 2, point3, and point 4, respectively). These electrodes form a grid of $2 \times 16$ ⁠, with a horizontal spacing of 3.2 cm and vertical spacing of 5.5 cm.

3.2. Experimental Method

3.2.1. The Measurement of Hydraulic Head Difference

The inflow and outflow reservoirs are connected with the fracture and porous media through a latex tube. After the flow is stable, the water level can be read out and the volume of water is collected by the graduate. The hydraulic head difference can be calculated.

3.2.2. The Measurement of the Tracer Concentration

The transport of the chloride is monitored by the high-density resistivity method. There are 32 electrodes available for measurements of solute transport. The fracture has sixteen electrodes, which are used for the measurement of resistance in fractured media. The porous domain has also sixteen electrodes, which are used for the measurement of resistance in porous media. Three sets of experiments were conducted: (1) on-line monitoring of solute transport with a hydraulic head difference (3.8 cm); (2) on-line monitoring of solute transport with a hydraulic head difference (6 cm); (3) on-line monitoring of solute transport with a hydraulic head difference (10 cm).

4. Results and Discussion

4.1. Flow Characteristics

It was found that the linear relationship between hydraulic gradient and specific discharge for Darcy’s law is no longer valid when the specific discharge is gradually increasing (the velocity of boundary point is $15.5 \times 10^{- 3}$ cm/s and $25.0 \times 10^{- 3}$ cm/s, respectively). There is a Darcy linear flow stage in front of the curved part, and there is the non-Darcy flow region after the back part curve. In this paper, linear fitting is adopted for the experimental data in the Darcy region and nonlinear fitting using equations (2) and (3) for the experimental data in the non-Darcy flow region. The fitting results are shown in Figure 2.

Figure 2 shows that both hydraulic gradient and specific discharge obey Darcy’s law when the specific discharge is small. With the increase of the specific discharge, the non-Darcy phenomenon appears. Equations (2) and (3) have achieved good fitting results for the fitting of non-Darcy flow in the back part of the curve. The fitting values of the coefficients in the two equations are shown in Table 1.

4.2. Changes of Conductivity under Three Different Hydraulic Head Differences

In order to investigate the change of conductivity in the pore-fractured media under the hydraulic head difference of 3.8 cm, 6 cm, and 10 cm, respectively, we analyzed changes of conductivity in the pore-fractured media. Three sets of head difference were considered in the experiments. It may need more sets to support this conclusion to avoid any abnormal data. We will take into account these questions in future experiments. And this may help to highlight the new understanding of this study. The results are shown in Figure 3.

From Figure 3, we can make several interesting observations. First, the value of conductivity increases and then decreases. But the peak value of conductivity becomes larger and larger as hydraulic head difference increases in the pore media. However, the peak value of conductivity decreases first and then increases as hydraulic head difference increases in fractured media.

The concentration of tracer was directly proportional to the conductivity (the reciprocal of the resistance). The value of resistance can be measured using electrical conductivity. In order to investigate the solute transport in the pore-fractured media under the hydraulic head differences of 3.8 cm, 6 cm, and 10 cm, both ADE and CTRW models are employed.

4.3. Simulation Results of BTCs under Hydraulic Head Difference of 3.8 cm

The experiment is finished under the hydraulic head difference of 5.3 cm in the pore and that of 1.5 cm in the fracture. Hence, the hydraulic head difference of pore-fracture media is 3.8 cm (⁠ $5.3 cm - 1.5 cm = 3.8 cm$ ⁠). The specific data of the inlet and outlet water level of pores and fractures are shown in Table 2.

In order to investigate the solute transport under the hydraulic head difference of 3.8 cm, the breakthrough curves (BTCs) under the hydraulic head difference of 3.8 cm in pore-fractured media were illustrated and fitted by ADE and CTRW-TPL models. The results are shown in Figure 4 and Table 3.

It was found that the good performance of CTRW-TPL as shown in Figure 4 is also reflected by the larger $R^{2}$ and smaller RMSE values than their counterparts resulting from the ADE model. Furthermore, the estimated $β$ values in CTRW-TPL were in the range of 0.90~1.20, which was less than 2, a value representing non-Fickian transport.

4.4. Simulation Results of BTCs under the Hydraulic Head Difference of 6 cm

The experiment is finished under the hydraulic head difference of 7.5 cm in the pore and that of 1.5 cm in the fracture. Hence, the hydraulic head difference of pore-fracture is 6 cm. The specific data of the inlet and outlet water level of pores and fractures are shown in Table 4.

The breakthrough curves (BTCs) under hydraulic head difference of 6 cm in pore-fractured media were illustrated and fitted by ADE and CTRW-TPL models. The results are shown in Figure 5 and Table 5.

It was found that the good performance of CTRW-TPL as shown in Figure 5 is also reflected by the larger $R^{2}$ and smaller RMSE values than their counterparts resulting from the ADE model. Furthermore, the estimated $β$ values in CTRW-TPL in the fracture media followed an increasing pattern towards 2.00. This implied that the transport would eventually evolve into a Fickian form under the hydraulic head difference of 6 cm.

4.5. Simulation Results of BTCs under the Hydraulic Head Difference of 10 cm

The experiment was finished done under the hydraulic head difference of 12.5 cm in the pore and that of 2.5 cm in the fracture. Hence, the hydraulic head difference of pore fracture is 10 cm. The specific data of the inlet and outlet water level of pores and fractures are shown in Table 6.

The breakthrough curves (BTCs) under the hydraulic head differences of 10 cm in pore-fractured media were illustrated and fitted by ADE and CTRW-TPL models. The results are shown in Figure 6 and Table 7.

It was found that the good performance of CTRW-TPL as shown in Figure 6 is also reflected by the larger $R^{2}$ and smaller RMSE values than their counterparts resulting from the ADE model. Furthermore, the estimated $β$ values in CTRW-TPL were in the range of 0.80~0.99, which was $0 \leq β \leq 1$ ⁠, a value representing non-Fickian transport.

5. Conclusion

In this paper, series of experiments in pore-fractured media were carried out for flow and solute transport with different hydraulic gradient conditions. The results reveal the following:

It was found that with the increase of specific discharge, the linear relationship between the hydraulic gradient and the specific discharge will deviate. Namely, Darcy flows stage in front of the curved part. However, there is the non-Darcy flow region after the curved part when the specific discharge is large enough.

The relationship between hydraulic gradient and specific discharge in the pore media can be fitted by either the Forchheimer or the Izbash equation in the back part of the curve. The Forchheimer and Izbash equations can be expressed by $J = - 1.46 \times 10^{- 4} v^{2} + 0.0177 v$ and $J = 0.0278 v^{0.796}$ ⁠, respectively. The Darcian equation can be expressed by $J = 0.0138 v$ in the front part of the curve. The same relationship in the fracture media can be fitted by either the Forchheimer or the Izbash equation in the back part of the curve. The Forchheimer and Izbash equations can be expressed by $J = - 1.06 \times 10^{- 5} v^{2} + 0.0018 v$ and $J = 0.0072 v^{0.546}$ ⁠, respectively. The Darcian equation can be expressed by $J = 0.00163 v$ in the front part of the curve.

The phenomenon from non-Fickian to Fickian then to non-Fickian was found in fractured media when hydraulic head difference ranges from 3.8 cm to 10 cm. The classical ADE model was incapable of capturing the long-tailing of BTCs. The fitting results of the CTRW-TPL model are far better than ADE in fitting the long tailing. Besides, the good performance of CTRW-TPL was also reflected by the larger $R^{2}$ and smaller RMSE values than their counterparts resulting from the ADE model.

Data Availability

The data supporting the conclusions of the study can be obtained through asking for the author.

Conflicts of Interest

This manuscript presents no conflicts of interest.

Acknowledgments

This research was partially supported by the National Natural Science Foundation of China (Grant nos. 41831289, 41877191, and 42072276).

Water level	Pore (cm)	Fracture (cm)
Inlet of water level	23.8	20.0
Outlet of water level	18.5	18.5
Hydraulic head difference	5.3	1.5

Water level	Pore (cm)	Fracture (cm)
Inlet of water level	23.8	20.0
Outlet of water level	18.5	18.5
Hydraulic head difference	5.3	1.5

Media	ADE				TPL
Media	$V_{A}$ mm/s	$D_{A} \times 10^{- 4}$ m²/s	$r^{2}$	RMSE	$V_{ψ}$ mm/s	$D_{Ψ} \times 10^{- 4}$ m²/s	$β$	$Lg (t_{1})$	$Lg (t_{2})$	$r^{2}$	RMSE
Left pore	0.5	0.6	0.729	0.152	0.3	0.6	1.20	$2.0 e + 01$	$5.0 e + 01$	0.942	0.032
Right pore	0.3	0.6	0.830	0.140	0.4	0.2	0.92	$1.0 e + 01$	$2.0 e + 01$	0.982	0.021
Left fracture	0.8	0.65	0.847	0.167	0.8	0.8	0.90	$3.0 e + 01$	$8.0 e + 01$	0.974	0.031
Right fracture	0.48	0.3	0.810	0.105	0.38	0.22	1.18	$2.0 e + 01$	$4.0 e + 01$	0.920	0.038

Media	ADE				TPL
Media	$V_{A}$ mm/s	$D_{A} \times 10^{- 4}$ m²/s	$r^{2}$	RMSE	$V_{ψ}$ mm/s	$D_{Ψ} \times 10^{- 4}$ m²/s	$β$	$Lg (t_{1})$	$Lg (t_{2})$	$r^{2}$	RMSE
Left pore	0.5	0.6	0.729	0.152	0.3	0.6	1.20	$2.0 e + 01$	$5.0 e + 01$	0.942	0.032
Right pore	0.3	0.6	0.830	0.140	0.4	0.2	0.92	$1.0 e + 01$	$2.0 e + 01$	0.982	0.021
Left fracture	0.8	0.65	0.847	0.167	0.8	0.8	0.90	$3.0 e + 01$	$8.0 e + 01$	0.974	0.031
Right fracture	0.48	0.3	0.810	0.105	0.38	0.22	1.18	$2.0 e + 01$	$4.0 e + 01$	0.920	0.038

Water level	Pore (cm)	Fracture (cm)
Inlet of water level	26.0	20.0
Outlet of water level	18.5	18.5
Hydraulic head difference	7.5	1.5

Experimental and Numerical Simulation Study of Flow and Solute Transport in Pore-Fractured Media Based on High-Density Resistivity Method

Abstract

1. Introduction