Solute Transport in the Element of Fractured Porous Medium with an Inhomogeneous Porous Block

: In this paper, the problem of solute transport in a fractured-porous medium taking into account the non-equilibrium adsorption kinetic is studied. The solute transport in fractured-porous medium consisting of two fractures and a porous block between them located in a symmetric form is considered. The problem is then solved numerically by using the ﬁnite difference method. Based on the numerical results, the solute concentration and adsorption ﬁelds in the fractures and porous blocks are shown in graphical form. The effect of adsorption on the solute transport in a fractured-porous medium is then analyzed. In the case of different parameters in two zones, asymmetric distribution of the solute concentration and adsorption is obtained. The nonlinear kinetics of adsorption leads to an increase in the adsorption effects, conversely slowing down the rate of the distribution of concentration of the solute in the ﬂuid.


Introduction
The problem of solute transport and the fluid flow in porous [1][2][3] and fractured-porous media (FPM) [4][5][6] in recent years has received great attention. This is due to various applications where the processes of solute transport and fluid flow in the FPM are the basis of industrial, pilot industrial works on the utilization of various wastes in underground reservoirs [7,8] and the intensification of oil production by water flooding with various dissolved substances into reservoirs with fractured porous collectors [9]. One of the rational ways to analyze these processes is to compile and study the hydrodynamic models of the process.
A large number of works [10][11][12] are devoted to the problems of hydrodynamic modeling of the processes of solute transport in porous media. However, the issues of solute transport in FPM have been studied relatively less. It should be noted that one of the pioneer works was done in [13][14][15], where an approach was taken to compile and analyze the hydrodynamic models of the solute transport in FPM. In the aforementioned works, the process of solute transport in FPM is described by a combination of convective-diffusion transport, which is dominant in fractures, and diffusion transport, which is dominant in porous blocks.
The problem of the transport of pollutants in porous, fractured media has attracted considerable attention by many authors, especially in the case of radioactive waste disposal in underground storage facilities [16][17][18]. Using the Laplace transform method [19], a number of solutions to the equations for the migration of radionuclides in fractures, as well as in surrounding porous media, were obtained. Most models suggest that radionuclides are transported through fractures mainly due to advection and dispersion, while they move into the surrounding porous matrix due to molecular diffusion. It is now generally accepted that fractures can play an important role in the transport of contamination in the groundwater system, because the permeability of the fractured system is often much greater than the permeability of the porous matrix. For such a transport case, the authors of [20][21][22] obtained analytical stationary state solutions for various boundary conditions. Complete solution of transport processes can be found in [23].
In many studies of groundwater flow in fractured media, it was believed that the effective permeability of the medium is dominated by fractures [24,25]. Although fractures are the main pathways for water flow and solute transport, diffusion in neighboring porous blocks can play a significant role in the overall solute transport through the medium. One-dimensional motion through a single fracture can be described generally by the Navier-Stokes equations for non-turbulent flows of a viscous incompressible fluid in the space between two parallel planes, neglecting the inertial terms.
In [26], an exact solution to the problem of solute transport in two vertical fractures and a porous block between them is presented. In fractures, convective transport and hydrodynamic dispersion are taken into account; however, in the porous block, only molecular diffusion is taken into account. On the common surface of the porous blocks and fractures, as well as inside the block, the substance is adsorbed.
The inhomogeneity of the porous medium can have a significant effect on the solute transport, both in porous and fractured-porous media. Numerous works have been performed to study the relationship between heterogeneity and the solute transport in such media [27][28][29]. The inhomogeneity of the porous block in terms of porosity and permeability leads to heterogeneity of the medium and in the sense of diffusion solute transport in it. In [30,31], the behavior of solute transport in inhomogeneous media was numerically studied with emphasis on the influence of permeability inhomogeneities. There are two main mechanisms that control the solute transport process-diffusion and advection. The dispersion effect promotes diffusion and mechanical mixing during the solute transport due to advection of the liquid. The effects of heterogeneity on the processes of solute transport were studied using various permeability distributions, which are characterized as continuous and discontinuous models [31]. For continuous distribution models, numerical simulations show that the tracer distribution is distorted by a local change in permeability, but the global distribution behavior still resembles a uniform distribution behavior.
The process of solute transport in a system of parallel fractures located in a porous matrix was studied in [32], which takes into account advection transport, molecular diffusion and mechanical dispersion in a fracture, molecular diffusion from fractures into the matrix, adsorption on the surface and inside the matrix, and radioactive decay of the substance. The migration of radionuclides in the FPM, where porous blocks were modeled as spheres with a certain radius, is considered in [33]. The resulting solution was compared with a similar solution obtained for a system of parallel fractures. A model of parallel fractures with equal openness and width was used to analyze the labeled solute transport in the FPM [34]. In the region of small times, a single-fracture model gives satisfactory results, since the solute does not penetrate deep into the matrix and the effect of a neighboring fracture is not felt. Diffusion is dominant in matrix with high porosity. In the case of adsorption transport, the micropores of the matrix act as large sinks due to the large adsorption surface.
In [35], a solution to the problem of solute transport in a fractured medium is presented, where the solute diffuses from the fractures into the porous matrix. The solute is inert and does not react with the matrix skeleton and the fracture. The equation in a fracture is the equation of convective (advective) diffusion in the one-dimensional case, and in porous blocks it is the one-dimensional diffusion type equation. An analytical solution to the problem used to interpret laboratory data obtained in [36].
In this paper, we study the solute transport in an element of the FPM, a porous block surrounded by two fractures, taking into account the layered inhomogeneity of the porous block and non-equilibrium adsorption of the solute using schematization [35][36][37]. For the fractures and the porous block, the equations of solute transport are written separately, taking into account the mass transfer between them. Adsorption of solute is considered nonlinearly non-equilibrium both in fractures and in porous block. Based on the numerical solutions of the problem, the concentration fields of the solute and adsorbed mass are obtained. The role of adsorption and heterogeneity of the porous block on the characteristics of solute transport is estimated.

Problem Formulation
Consider a FPM element consisting of a porous block (matrix), which is surrounded by both sides with fractures ( Figure 1). Fractures and porous block (and its parts R 1 and R 2 ) are considered one-dimensional and semi-infinite. The field of study of the problem consists of two parts: where θ 1 and θ 2 are the porosity coefficients in R 1 and R 2 (dimensionless quantities); D * 1 , D * 2 (m 2 /s) are the coefficients of effective diffusion in the zones R 1 and R 2 , respectively, which characterize diffusion properties of the parts of the porous block. In the porous block between the two fractures, diffusion mass transfer occurs in y direction. The equations of convective-diffusive solute transport, taking into account the adsorption and mass transfer of the solute, are written separately for fractures and a porous block in the following form [38][39][40]: ∂c m ∂t ∂c m ∂t where c m = c m (t, x, y) is the concentration of the solute in the matrix, fractures (m/s); ρ 1 , ρ 2 are the density of the medium in the R 1 and R 2 (kg/m 3 ); b 1 , b 2 are the width of fractures (m); t is the time (s).
We use the following kinetic equations in equilibrium state corresponding to the Freundlich isotherms [41,42] ∂s where k f (n) , k m(n) are the equilibrium constants of the Freundlich equation in the fracture and in the matrix (m 3 /kg; α f (n) , α m(n) are the adsorption rate coefficients characterizing the intensities of the adsorption processes in the fracture and in the matrix, respectively (s −1 ); N is constant (0 < N < 1); n is the index coefficient of the equations (n = 1 for R 1 , n = 2 for R 2 ). From Equations (5) and (6) at t → ∞, we obtain the equilibrium equations (Freundlich isotherm): It is worth mentioning that, when N = 1 at equilibrium state we obtain the linear law of adsorption isotherms (Henry Isotherm).
Let us assume that, at the initial moment the fractures and the porous block filled with pure fluid. At t > 0, the fluid with a dissolved substance is fed into the fractures, for general case with different concentrations of c 01 and c 02 , respectively. We assume that, for the time ranges used in this work, the concentration fields do not reach x = ∞. At the boundary of the porous block y = b 3 2 , the concentration of the solute remains continuous. With this formulation, the initial and boundary conditions of the problem have the following form:

Numerical Solution Algorithm
The system of Equations (1)-(4) taking into account (5) and (6) under the conditions (7)- (14) is solved by using the finite difference method [43,44]. To do this, in the region On this grid, Equations (1)-(4) are approximated as where c f k i , cm k i,j are the net values of concentrations c f (t, x) and c m (t, x, y) at net points (t k , x i ) and t k , x i , y j , respectively.
Equations (5) and (6) are approximated as follows where s f k i , sm k i,j are net values of concentrations s f (t, x) and s m (t, x, y) at net points (t k , x i ) and t k , x i , y j , respectively.
The initial and boundary conditions are approximated as follows c f k+1 where I is taken large enough so that the concentration profiles do not reach the selected conditional boundaries of the regions. After approximation, Equations (15) and (17) are reduced to where The Equations (16) and (18) are also reduced to Similarly, after approximation, Equations (19)- (22) are respectively reduced to where The systems of Equations (31) and (32) are solved by the sweep method in the i and j directions, respectively [44,45]. We use the following relations where α j+1 are indeterminate coefficients. To solve (33) and (34), the sweep method in the i and j directions is also used, respectively. For this, we use the following relations j+1 , β (4) j+1 are indeterminate coefficients. When using (39) from (31), we obtain the following recurrence formulas for determining the coefficients α Similarly, we obtain the following recurrence formulas for other sweep coefficients in (40)-(42) Starting values of coefficients α i+1 are determined based on the condition (24): α are determined based on the conditions (26) and (27): α According to recurrent formulas (44) and (46) we determine the values of α (2) j+1 , β (2) j+1 (j = 1, L − 1 ) and α (4) j−1 , β (4) j−1 (j = J − 1, L + 1 ). From condition (29), we find cm k+1 i,L in following form . From Equations (35)- (38), we find s f k+1 i and sm k+1 i,j in the zones R 1 and R 2 , respectively. We note, that sufficient stability conditions to apply Thomas' algorithm to Equations (31)- (34) are satisfied.
Knowing the concentration field, we can determine the relative current flow rate of the solute through the common boundary of the fracture and the porous block in the zones R 1 and R 2 and in the net approximation The total relative flow rates across the common boundaries (y = 0 and y = b 3 ) are defined as and the relative summary flow rates as: Integrals (50), (51) are calculated numerically based on numerical values of Q k+1 i . In this case, the discretization of the integration region is carried out in accordance with the introduced grid.

Results and Discussion
The presented problem is numerically solved using the finite difference method. Some calculation results are shown in Figures 2-8 for various model parameters. The following values of parameters are used: [35][36][37]. We used the following mesh parameters h 1 = h 2 = 0.1, τ = 1. Figure 2 shows surfaces c m /c 0 and s m for linear kinetic adsorption (in the equilibrium case, the corresponding Henry isotherm) with the same characteristics of fractures and porous block. In Figure 3, the same surfaces given for nonlinear kinetic adsorption (in equilibrium case, corresponding to the Freundlich isotherm) of the solute. Comparing the data presented, it can be seen that the surfaces of the distribution of concentration and adsorption of the solute from the two fractures are symmetric, therefore, the relative current solute flow rate Q i , i = 1, 2 from the first (i = 1 in the zone R 1 ) and second (i = 2 in the zone R 2 ) fractures in the porous block will be the same. This conclusions are similar with the results obtained by Grisak and Pickens [35]. Figures 4 and 5 show the distribution fields of the concentration of the solute and the adsorption field in the medium for two cases of nonequilibrium adsorption. In this case, it was posed that in the region R 2 solute transport without taking into account adsorption, as well as in the R 1 zone, the calculation results are presented with an increase in the adsorption coefficient k f 1 = k m1 = k and α f 1 = α m1 = α. As can be seen from these figures, that changing the values of parameters k f 1 , k m1 from k f 1 = k m1 = 3 × 10 −5 m 3 /kg to k f 1 = k m1 = 3 × 10 −4 m 3 /kg, α f 1 , α m1 from α f 1 = α m1 = 4 × 10 −5 s −1 , to α f 1 = α m1 = 4 × 10 −4 s −1 , and N from N = 0.88 to N = 0.83 (increasing the adsorption effect), leads to a reduction in the distribution of the concentration field in the zone R 1 (Figure 4). The maximum values at t = 6000 s for non-equilibrium adsorption corresponding at equilibrium conditions to the Freundlich isotherm are greater than for the Henry isotherms (Figure 5a,b).
The dynamics of the relative solute flow rate from the first (Q 1 ) and second (Q 2 ) fractures into the porous block for two cases of kinetic adsorption is shown in Figure 6. With an increase in adsorption, the distribution zone Q 1 decreases at small x, then the opposite tendency is observed. An increase in adsorption in the zone R 1 leads to an increase in the total and summary solute flow rate in the zone R 1 (Figures 7 and 8).
As can be seen from Figure 7, both cases of linear and nonlinear adsorption kinetics, a non-monotonic kinetics of the total relative flow rate of the solute from the fractures to the porous block are obtained. Such non-monotonicity is characterized by the fact that in the initial stages of the transport process in a zone close to the input section of the medium, which is in the initial sections of the medium, relatively large concentration gradients are formed. As the solute transport continues, the porous block begins to saturated with the solute and concentration gradients begin to decrease. This leads to non-monotonic dynamics of the relative flow rate. Such non-monotonicity is observed at local time scales. In the summary relative solute flow rates non-monotonicity is not observed, since they represent an integral temporal characteristic (Figure 8). The lack of non-monotonicity in Q 1sum , Q 2sum can also been explained by the fact that the solute mass transfer through the boundaries of fractures and the porous block over the entire length of the boundaries, for all times occurs in one direction-from fractures to the porous block.

Conclusions
The problem of solute transport in a FPM element, consisting of two fractures and a porous block between them, with the effect of the non-equilibrium kinetic adsorption of the solute on the transfer characteristics is considered. Based on numerical calculations, it was found that the adsorption in the R 1 zone leads to a decrease in the concentration field in the zone.
In this study, two cases of non-equilibrium adsorption are analyzed and its effect on the solute transport in the element are shown. For different effective diffusion coefficients in the porous block, an asymmetric distributions of the concentration and adsorption of the solute are obtained.
It is found that the layered inhomogeneity of the porous block leads to an asymmetric distribution of the concentration and adsorbed solute mass fields. The nonlinear kinetics of adsorption, ceteris paribus, leads to increase the adsorption effects. Thus, it slows down the rate of the distribution and the concentration of the solute in the fluid. A non-monotonic dependence is found for the current and total relative solute flow rates through the common boundary of the fractures and the porous block. Total relative solute flow rates sharply increases in small values of time, then over time it decreases slowly. This is explained by the formation of large concentration gradients in the initial sections of the medium in the initial stages of the process and their subsequent decrease. Non-monotonicity is not observed for the summary relative flow rate.
Thus, this approach together with the schematization of the porous block as a layered inhomogeneous medium by diffusion properties allows a qualitative and quantitative analysis of transport processes in inhomogeneous FPMs.