The Effects of Pore Geometry on Late Time Solute Transport with the Presence of Recirculation Zone

Abstract

The solute transport process in porous media is central to understanding many geophysical processes and determines the success of engineered applications. However, fundamental understanding of solute transport in heterogeneous porous media remains challenging especially when inertial effects are significant. To address this challenge, we employed direct numerical simulations in a variety of intrapore geometries at a high Reynolds number (Re = 10) flow regime, where recirculation zones (RZs) are present with significant inertial effects. We find that the volume of RZs depends on pore geometries. Moreover, RZs serve as an immobile domain that can trap and release solutes that lead to non-Fickian transport, characterized by the early arrival and heavy tailing of breakthrough curves and bimodal residence time distributions (RTDs). Lastly, the late time portion of RTDs is fitted to the power law function with determined exponent n, where n depends on the pore geometries and consequently the volume of RZs. Our study sheds light on the mechanisms of an immobile zone on the solute transport, especially improving our understanding of late time transport tailing in pressurized heterogeneous porous media.

1. Introduction

Fluid flow and transport processes in subsurface porous media assume the applicability of linear rate laws, i.e., Darcy’s law and Fick’s law, respectively [1]. This assumption normally holds in low pressure gradient environments where inertial effects are negligible [2]. However, in some pressurized or engineered geological settings, e.g., the subsurface sequestration of CO₂, the applicability of linear rate laws is questionable due to non-negligible inertial effects [1,3,4,5].

Indeed, the transition from Darcy to Forchheimer flow with increasing inertial effects has been widely acknowledged based on numerous laboratory and field experiments and numerical simulations [4]. This suggests the invalidity of linear rate law (i.e., Darcy’s law) for describing fluid flow in porous media. With the development of computational ability and high-resolution devices, researchers found that the growth of recirculation zones (RZs) in porous media due to increasing inertial effects can shrink the main flow channel, resulting in a reduction of effective permeability for porous [6,7] and fractured [8,9,10] media.

Furthermore, the development and growth of RZs have profound effects on solute transport by trapping and afterward releasing solutes back to the main flow channel [11], where RZs can be regarded as one type of immobile zones that retard solutes in the flowing systems. Additionally, the mass transfer process between the immobile zone (i.e., RZs) and mobile zone (i.e., main flow channel) can lead to a late time transport tailing process [12]. Such tailing process has been widely recognized in many geological environments, including the hyporheic zones [13,14], the nested groundwater systems [15], and the fractures [16,17,18].

To better characterize transport tailing, many macroscopic models have been proposed, such as multi-rate mass transfer [12], continuous time random walk [19], and mobile-immobile model [20]. These models are mostly applicable when single-modal residence time distributions (RTDs) are present. However, when inertial effects are profound, the growth of RZs might lead to bimodal or multi-modal RTDs, which has been demonstrated in fractured media [10]. In such circumstances, how to understand and accurately capture the tailing process is largely unknown, especially for heterogeneous porous media in pressurized environments [21,22].

To address the above-mentioned challenge, we used direct numerical simulations to solve fluid flow and solute transport in microscopic intrapore geometries with varying shapes, where RZs were present due to significant inertial effects when Reynolds number (Re) = 10; this helps us to exclusively diagnose the role of geometry-dependent RZs on solute transport in porous media, particularly on the late time tailing process.

2. Materials and Methods

2.1. Diverse Intrapore Geometries for Numerical Simulations

The first order heterogeneity of porous media can be represented by the converging-diverging pore geometries [23,24]. To fully scrutinize the effects of pore geometry on the solute transport process, a two-dimensional (2-D) axis-symmetric pore geometry was used to represent flow and transport in a three-dimensional pore (Figure 1a–c) [6]. Ideally, the pore geometries can be characterized by a non-dimensional hydraulic shape factor β [6,7]:

$$\beta =\frac{S_A\times L}{V}$$

(1)

where S_A is the total surface area of microscopic pores, L is the domain length in the x direction, and V is the pore volume (Figure 1c).

Here, a broad range of pore geometries (β = 5.5, 6.9, 8.5, 12.3, 15.5, 18.9, and 21.9) was particularly designed with L = 2.2 mm and the pore throat radius R = 0.2 mm (Figure 1c). Small β cases are cavity-like pores whereas large β cases are slit-like pores. These geometries were used to identify the geometry-dependent flow and transport behaviors, given that β is related to the growth of RZs induced by significant inertial effects within intrapore [7], which is further discussed below.

2.2. Direct Numerical Simulation of Fluid Flow via Navier-Stokes Equations

The fluid flow in porous media is fundamentally governed by the Navier-Stokes equations:

$$\{\begin{array}{c}\rho \left(\mathit{u}·\nabla \right)\mathit{u}=-\nabla p+\mu {\nabla }^2\mathit{u}\\ \nabla ·\mathit{u}=0\end{array}$$

(2)

where u = [u, w] is the velocity vector in the two-dimensional pores, p is pressure, ρ (1 × 10³ kg/m³) and μ (1 × 10⁻³ Pa·s) are, respectively, water density and viscosity.

We applied a constant flux (q = 0.05 m/s) at the inlet (right) and the Dirichlet boundary condition at the outlet (left); this drove fluid flow from right to left (Figure 1d). The symmetric boundary condition was used at the bottom of the axis-symmetric pore with no-slip walls for the rest of the boundaries (Figure 1c). Note that, the specific flux was used to ensure inertial effects were profound where Reynolds number (Re) = 10. Re is defined as:

$$Re=\frac{\rho qR}{\mu }=10$$

(3)

The pore geometries were spatially discretized in a way, such that the inner mesh size was larger than those around the boundaries. This ensured the accuracy of numerical results. In total, we had ~1.6 × 10⁵ triangular elements and numerical results were insensitive to a further refinement of the domain. The implementation of numerical simulation of fluid flow was conducted in COMSOL Multiphysics, a commercial finite-element based software. Each run took about a wall-clock of a few minutes.

2.3. Automatic Delineation of Recirculation Zone in Single Pores

The resultant flow fields can be used to further analyze the spatial distribution of RZs in diverse pore geometries. This was achieved by using the zero-flux method proposed by previous studies [10,25]. That is, we exported the velocity field from COMSOL with a regular cell size (~1 × 10⁻⁴ m), and the resolved velocity field was imported into Matlab code for accurately detecting the interface between RZs and the main flow channel (Figure 2). Moreover, the volume of RZs (V_RZ) was estimated to possibly establish a relationship between V_RZ and the degree of transport tailing as explained below.

2.4. Direct Numerical Simulation of Solute Transport via Advection-Diffusion Equation

Solute transport in microscopic pores is essentially governed by the advection-diffusion equation:

$$\frac{\partial C}{\partial t}+\nabla ·\left(\mathit{u}C\right)=\nabla ·\left(D_m\nabla C\right)$$

(4)

where C is concentration, t is time, D_m is the molecular diffusion coefficient (2.03 × 10⁻⁹ m²/s). Here, we applied a step function at inlet using the Dirichlet boundary condition, i.e., a constant concentration (C₀ = 1 kg/m³) at inlet (right in Figure 1c) when time > 0. Moreover, an open boundary condition with zero dispersive/diffusive flux was set at the outlet (left in Figure 1c). All other boundaries were set as no-flux boundary conditions. The entire domain was free of solute at time = 0.

Like fluid flow, the solute transport process was implemented in COMSOL Multiphysics using the same meshing scheme (Figure 3). We used a small time step (~1 × 10⁻³ s) to ensure the accuracy of numerical results considering a relatively high velocity field with Re = 10. Each run for transient solute transport took about a wall-clock of a few hours depending on the pore geometry, i.e., β.

2.5. Breakthrough Curves and Bimodal Residence Time Distributions

Based on solute transport simulation results, the breakthrough curves (BTCs) at the outlet (left in Figure 1c) was estimated by quantifying the dimensionless flux-weighted concentration (C′) over time (Figure 4a):

$$C^{\prime }=\frac{uC}{Q{C}_0}$$

(5)

where C is numerically-derived concentration over time at outlet solved by the Equation (4), u is computed velocity at outlet solved by the Equation (2), Q is specified discharge (Q = 1 × 10⁻⁵ m²/s), C₀ is the specified concentration at inlet.

Moreover, to better demonstrate the late time transport tailing process, the BTCs were converted to residence time distributions (RTDs); this was achieved by taking the time derivative of C′ and plotting it against time (Figure 4b). Note that the time can be normalized by using pore volume (PV), which is defined as:

$$PV=\frac{V}{Q}$$

(6)

where V is volume of microscopic pores.

3. Results and Discussion on Solute Transport in Diverse Intrapore Geometries

3.1. Development of Recirculation Zones in Microscopic Intrapore

The development of RZs can be formed in heterogeneous porous media [6,11]. However, the size of RZs would remain unchanged without considering inertial effects. By including inertial effects, RZs would gradually grow, where more energy could be dissipated [11]. The growth of RZs might reduce the apparent permeability of porous media and result in a nonlinear relationship between flow rate and pressure gradient, i.e., nonlinear fluid flow [6,7]. For instance, a previous study recently found that the threshold for identifying the transition from Darcy flow to nonlinear (Forchheimer) flow depends on intrinsic permeability and slip boundary condition [4], where intrinsic permeability basically depends on pore geometry, i.e., β.

In this study, Re was set to be constant (Re = 10) such that we can exclusively identify the effects of pore geometry on fluid flow and solute transport, whereas RZs were remarkably present (Figure 2). Regarding the flow field, we find that V_RZ decreases with an increase of β (Figure 2 and Figure 5). This is essentially attributed to a fact that a small β indicates ample room for the growth of RZs, while a large β suggests much less room for the development of RZs (Figure 1a). Another typical feature of a flow field with RZs is that the magnitude of velocity in the main flow channel is much larger than that in RZs, which is congruent with previous studies [11,24,26].

Interestingly, the automatic delineation of RZs is demonstrated to be a viable tool for identifying the interface between RZs and the main flow channel considering complex flow behavior with more than one RZs directly connecting each other. For example, there were two RZs rotating in exactly opposite directions, i.e., clockwise and anti-clockwise, in a microscopic pore (Figure 2c,d).

3.2. Trapping and Releasing of Solutes in Recirculation Zones

Solute travels along with fluid following advection and diffuses induced by the concentration gradient [1]. Due to the distinct velocity field between RZs and the main flow channel (Figure 2), solutes migrate orders of magnitude faster in the main flow channel than those trapped in the RZs (Figure 3). That is, solute transport is mainly advection-dominated in the main flow channel, whereas the transport process tends to be diffusion-dominated or slightly advection-dominated in the RZs depending on pore geometries, i.e., β. Consequently, for all cases with a broad range of β, we observe that solutes promptly sweep over the main flow channel. Afterward, solutes slowly migrate from the main flow channel, across the interface, and diffuse into RZs, i.e., solutes are temporarily trapped in the RZs.

The transport process within RZs depends on β. Specifically, when β is small, the volume of RZ and the magnitude of velocity is relatively larger than in the cases with large β. As a result, solutes would slightly advect along the boundary of RZs and then diffuse into RZs’ center for the small β cases (Figure 3). This is consistent with previous studies on fractures [25]. In contrast, solutes are mostly diffusing into the RZs, especially for the large β cases with two hydraulically-connecting RZs (Figure 3).

The trapped solutes in RZs would eventually release back to the main flow channel given a long time period [8]. The trapping and releasing of solutes can affect the late time solute transport, especially the tailing process as suggested by many previous studies [8,16]. Consequently, the contrasting solute transport behavior due to different RZs in diverse pore geometries is expected, i.e., RZs impose appreciable impacts on the tailing process. Take CO₂ subsurface sequestration, for example, the dissolved CO₂ might be trapped by the RZs. In some cases, the trapping process can even create the possibility of two dissolution stages due to this difficulty in solute transport as illustrated in the previous study [5].

3.3. Non-Fickian Transport with Bimodal Tailing Process Depending on Intrapore Geometry

The information of the transport process is collectively encrypted in BTCs and RTDs [12,27]. For the generalized mobile-immobile systems, the non-Fickian transport is expected, because the transport process cannot be accurately described by the Fickian theory in such systems [12,20]. Our numerical simulations with the presence of RZs support this. That is, two typical transport behaviors that directly suggest non-Fickian transport are observed in this study: (a) The early arrival is clearly demonstrated in the BTCs and RTDs (Figure 4). This is because Fickian transport normally converges to the point (1, 0.5) in the space of dimensionless concentration and pore volume of BTCs, and it reaches its peak value when time = 1 pore volume in RTDs. However, none of the cases shows the Fickian behavior demonstrating early arrival behavior. (b) The late time transport tailing is also demonstrated in BTCs and RTDs as further discussed below.

Due to the presence of RZs (i.e., a kind of immobile domain) at a high Re flow regime for all cases with diverse pore geometries, BTCs all exhibit tailing behavior where concentration takes a relatively long time to reach its asymptotic value (Figure 4a). In fact, the pattern of tailing depends on β (Figure 4). That is, the cavity-like pores with a small β demonstrate a heavier tailing than the slit-like pores with a large β as shown in BTCs.

The RTDs is an alternative method to better demonstrate the tailing behavior [28]. Clearly, the transport tailing is obviously observed in RTDs (Figure 4b), which is exemplified by the bimodal RTDs. The bimodal RTDs suggest the impacts of RZs on the solute transport process, where the secondary peak is primarily attributed to the releasing of solutes in RZs. This phenomenon is consistent with previous findings for fractures using the particle tracking method [10].

We further quantified the degree of tailing by fitting a power law function (y = x⁻ⁿ) to the late time portion of RTDs after the secondary peak (black line in Figure 4b). Specifically, we selected the last 40 points excluding the very end of 30 points. The fitting power law exponent n indicates the degree of tailing, i.e., a large n suggests a lighter tailing, whereas a small n suggests a heavier tailing. By plotting n against β, we find that n increases with β (Figure 5), which is exactly opposite to the relationship between the volume of RZs (V_RZ) and β (Figure 5). The results suggest that V_RZ largely controls the tailing behavior, where large RZs would have more capacity to capture solutes and thus lead to a heavier tailing. Our new finding of relationship between n and β can be possibly used for upscaling transport process via deterministic method [29] for complex porous media in nature [30].

4. Conclusions

To better understand solute transport in pressurized heterogeneous porous media, we used direct numerical simulations for simulating fluid flow and solute transport in microscopic intrapore geometries when inertial effects are significant, where the pore shape geometry is represented by the non-dimensional hydraulic shape factor β. The recirculation zones (RZs) are present at a high Reynolds number flow regime for all β cases. We find that the volume of RZs depends on β. Moreover, the presence of RZs significantly affects solute transport by trapping and releasing solutes, resulting in non-Fickian transport with prominent early arrival and long tailing of solute breakthrough curves and bimodal residence time distributions. The degree of tailing characterized by the power-law exponent increases with β, which can be fundamentally explained by the volume of RZs that decreases with β. Our study helps to explain the heavy tailing of solute transport observed in heterogeneous porous media, especially in pressurized environments with bimodal residence time distributions.