1 Introduction
Two-phase flow is encountered widely in the field of hydrology because soil and groundwater pollution by dense non-aqueous phase liquids (DNAPLs), such as chlorinated solvents (e.g., trichloroethylene (TCE)), constitute a large and serious environmental problem (Bano et al., 2009; Bohy et al., 2004). When infiltrating through porous media, migration of DNAPLs can result in highly fingered fluid distributions. The occurrence of fingering is caused by flow instabilities due to differences in viscosity and density between DNAPL and water (Khataniar and Peters, 1992; Lindner and Wagner, 2009; Riaz and Tchelepia, 2006; Tehelepi et al., 1994). These fingers propagate rapidly, resulting in an early breakthrough compared to the displacement of a stable front. Generally, two approaches are used to study the displacement process through porous media. The first approach uses the continuum description of the porous medium associated with macroscopic laws (e.g., Muskat's law) to represent the volumetric flow vector for each fluid phase on the REV scale (Barenblatt and Patzek, 2003; Bear, 1972; Galusinski and Saad, 2009). The second approach is based on the microscopic description of the pore geometry and on the physical laws of flow and transport within the pores (Blunt and King, 1990; Laroche and Vizika, 2005; Lowry and Miller, 1995). One of the common tools that has been used to study the displacement process is pore-scale network modelling, which uses simple geometric shapes to represent the pore space. Many network models have been developed to study a wide range of displacement processes, including drainage and imbibitions (Hilpert and Miller, 2001; Joekar-Niasar et al., 2009; Nordhaug et al., 2003; Patzek, 2000; Singh and Mohenty, 2003). The most critical aspect of network model construction is the appropriate choice of the pore structure and its capacity to represent the real physical system (Celia et al., 1995; Raoof and Hassanizadeh, 2009). Generally, an equivalent network to the real porous medium is used. This network can be based on statistical distributions and correlations of basic morphologic parameters, such as pore bodies and pore throats. Typical parameters used for building an equivalent network include pore body size distribution, pore throat size distribution and throat length distribution. In the majority of models, the parameters are fictive and do not represent a realistic porous medium. Consequently, these models cannot effectively analyse and reproduce physical experiments.
Different approaches exist for identifying pore structures (Saad et al., 2010). One of the common approaches is the inverse method. Pressure-saturation curves measured in the laboratory for a given two-phase fluid system in a porous medium are fitted using theoretical relationships to obtain the equivalent geometric parameters of the network model. Other experimental approaches are direct methods, such as the use of mercury porosimeters to assess the pore-size distribution of the porous medium. Because direct measurements of pore space structure are generally costly and time consuming, recent conceptual approaches have been developed to derive the void size distribution (VSD) from the grain size distribution (GSD). The GSD is widely available for many soil types and can be routinely determined in the laboratory. A novel approach based on the geometric proprieties of the packing of spheres has been elaborated (Dodds, 1980; Geoffrey and Mellor, 1995). This approach computes the VSD resulting from a specific packing of spheres. The shape and size-distribution of the voids within a random packing of an array of spheres and the effect of GSD on the packing have received considerable attention (Rouault and Assouline, 1998).
In this article, we present a novel network model for simulating drainage processes in a homogeneous porous medium. The developed model is based on the same basic flow characteristics as the model introduced by Nordhaug et al. (2003). However, our model considers capillary effects at water-DNAPL interfaces. Capillary pressure is assumed to be located at the boundary between a pore body and a pore throat. To quantify the flow rate in a pore throat, which is connected to a pore body separated by a water-DNAPL interface, we used the Washburn equation. Furthermore, to obtain a predictive network model, the model parameters used to generate both a pore-body size distribution and pore-throat size distribution were quantified using a geometrical model of an array of spheres based on the grain-size distribution combined with a probabilistic approach. The numerical model was then applied to a laboratory drainage experiment (DNAPL/water) that was conducted on a sand-filled column under stable and unstable flow conditions. The pressure at the column entrance section and the arrival times distribution of the DNAPL front in a cross-section, measured during the displacement process, were compared to the numerical results. A sensitivity study of the size of pore throats and the size of the lattice was used to analyse the influence of these model parameters on the displacement of water by DNAPL.
2 Network model
Pore-scale network models typically represent the pore space of the medium using simplified geometries. Within these geometric representations, the models solve equations to explicitly track the location of all fluid–fluid interfaces within the network. We developed a pore-throat network model to simulate the drainage process in saturated porous media. Drainage occurs during an immiscible displacement of a wetting fluid (water) by a nonwetting fluid, such as trichloroethylene.
2.1 Geometry of the network
The pore network that is simulated in the model is a cubic lattice with constant node spacing. The pore bodies are represented by spheres. Similarly, pore throats are cylindrical with a circular cross-section (Fig. 1). The locations of fluid-fluid interfaces are restricted to the connections between the pore bodies and pore throats. The pore space is constructed by assigning radii to the lattice elements using a Weibull distribution (Weibull, 1951), which is frequently used to construct network analogs of porous media. The distribution is attractive due to its flexibility and to the fact that it is based on fixed parameters that describe the characteristics of the sizes of pore bodies and pore throats of the network. A minimum, mean and maximum radius (r_{min}, r_{mean} and r_{max}, respectively) specified for both pore bodies and pore throats are therefore required to randomly generate the pore body-size distribution and pore-throat size distribution in the network. The probability distribution that was used is given by:
$$g\left(x\right)=\left\{\begin{array}{l}\frac{x}{{x}_{2}^{2}}\mathrm{exp}\left(-\frac{{x}^{2}}{{x}_{2}^{2}}\right)\text{\hspace{0.28em}}\text{for}\text{\hspace{0.17em}}x\le {x}_{3}\\ 0\text{\hspace{0.28em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{0.17em}}\text{for}\text{\hspace{0.17em}}\text{x}\text{\hspace{0.17em}}>{\text{\hspace{0.17em}x}}_{3}\end{array}\right.$$ | (1) |
The porosity of the equivalent porous medium is assumed to be equal to the areal porosity of a horizontal transect of the network model. The reason for using areal porosity is that volumes do not scale as the third power of the radius, whereas areas scale as second power of the radius. Assuming that the pore throat volumes are negligible compared to the total volume of pore bodies, the porosity ø[-]of the pore-throat network is given by:
$$\varphi =\frac{{\sum}_{\text{i}}{A}_{\text{pore,}i}}{{n}_{x}{n}_{y}{n}_{z}{L}^{2}}$$ | (2) |
2.2 Governing flow equations of the pore-throat model
The developed network model makes the following assumptions: (i) the fluids used are Newtonian, incompressible and immiscible; (ii) the flow is laminar; and (iii) the local capillary pressure in the pore bodies is negligible, resulting in only one mean pressure within a pore body, independent of the phase saturation of that pore body. The mass balance equation for each pore is therefore written as follows:
$$\sum _{j}{Q}_{\mathrm{ij}}={Q}_{i}^{e}$$ | (3) |
If no meniscus is present at pore throat entry, the flow rate Q_{ij} through the pore throat is given by Poiseuille's Law
$${Q}_{ij}=\frac{\pi {r}_{ij}^{4}}{8L{\mu}_{\alpha}}\left({P}_{i}-{P}_{\text{j}}-\left({z}_{i}-{z}_{j}\right){\rho}_{\alpha}g\right)$$ | (4) |
If a meniscus is present at pore throat entry, the flow rate Q_{ij} through the pore-throat is expressed by the Washburn equation (Washburn, 1921):
$${Q}_{ij}=H\u2329{P}_{i}-{P}_{j}-\left({z}_{i}-{z}_{j}\right){\rho}_{\alpha}g-{P}_{c}\left({r}_{ij}\right)\u232a\frac{\pi {r}_{ij}^{4}}{8L\overline{\mu}}\left({P}_{i}-{P}_{j}-\left({z}_{i}-{z}_{j}\right)\overline{\rho}g\right)$$ | (5) |
P_{c}(r_{ij}) [ML^{−1}T^{−2}] used in Eq. (5) denotes the capillary entry pressure assigned to pore throat connecting pore i and j that needs to be exceeded by the DNAPL to enter the water filled pore throat. The capillary entry pressure is given by Young Laplace's Law as follows:
$${P}_{c}\left({r}_{\mathrm{ij}}\right)=\frac{2\sigma cos\theta}{{r}_{ij}}$$ | (6) |
$H\u2329X\u232a$ is a Heaviside function that stands for the displacement condition of the interface in a pore throat. The throat can be blocked by capillarity if the capillary entry pressure (P_{c} (r_{ij})) has not yet been reached. The Heaviside function is defined as:
$$H\u2329X\u232a=\left\{\begin{array}{c}\hfill 1\text{\hspace{1em}}\text{for}\text{\hspace{0.17em}}X>0\hfill \\ \hfill 0\text{\hspace{1em}}\text{for}\text{\hspace{0.17em}}X\le 0\hfill \end{array}\right.$$ | (7) |
In situations where a meniscus is present at boundaries between a pore body and a pore throat, the identity of the fluid flowing through the pore throat (DNAPL or water, or capillary blocked) can only be determined after computation of the pressure value of the whole network. Therefore, we use an “effective” viscosity $\overline{\mu}$ and an “effective” density $\overline{\rho}$ for the displacing fluid in Eq. (5), which are defined by:
$$\begin{array}{l}\overline{\mu}={\mu}_{nw}{S}_{\mathrm{ij}}+\left(1-{S}_{\mathrm{ij}}\right){\mu}_{\text{w}}\\ \overline{\rho}={\rho}_{\text{nw}}{S}_{\mathrm{ij}}+\left(1-{S}_{\mathrm{ij}}\right){\rho}_{\text{w}}\\ {S}_{\mathrm{ij}}=\left({S}_{i}+{S}_{j}\right)/2\end{array}$$ | (8) |
Introducing the appropriate expression of flow rate Q_{ij} in Eq. (3) leads to a set of algebraic equations with the pore body pressures as unknowns. This set of equations is written in matrix form as follows:
$$\text{AP}=\text{B}$$ | (9) |
The saturation field is updated after computation of the pressure terms by quantifying a characteristic minimum time (hereafter called T_{min}), which corresponds to a minimum filling time for each of the pore bodies. The calculated pressure field is used to quantify fluid fluxes through the pore throats. The ratio of free volume in a pore body (not invaded by displacing fluid) to the flow rate at a given interface location gives the filling time needed for that pore body to be invaded. The saturations in each pore body are then updated using:
$${S}_{i}^{n+1}={S}_{i}^{n}+\frac{{Q}_{i}\times {T}_{\mathrm{min}}}{{V}_{i}}$$ | (10) |
During one time step, only one pore body reaches full nonwetting phase saturation. This pore body generates the additional interfaces, which are located at all connecting pore throats that are filled with wetting fluid. Newly created interfaces are tested for stability, and appropriate expressions for flow rate Q_{ij} are then identified. This leads to a modified matrix of conductance terms that is reformulated for each time step in the set of equations (Eq. (9)). The field pressure is computed and the procedure for updating saturation in pore bodies is repeated; this provides a transient response for pressure and saturation fields.
In the case of a single-phase flow, the intrinsic permeability k [L^{2}] of the equivalent porous medium can be furthermore computed using Darcy's Law as follows:
$$k=\frac{{\mu}_{\alpha}{Q}_{i\text{,total}}^{e}h}{A\left({\Phi}_{\text{up}}-{\Phi}_{\text{down}}\right)}$$ | (11) |
3 An approach based on grain-size distribution to derive the radii of pore bodies and pore throats
3.1 Theory
To generate a distribution of pore-body sizes and pore throat-sizes for the network, we used a statistical distribution that specifies values of the minimum, mean and maximum size (Section 2.1). Assessment of these network parameters is essential to deliver a quantitative model and to obtain reliable simulations of laboratory experiments. Our modelling concept is based on geometric properties of packing spheres, which are combined with a probability approach proposed by Rouault and Assouline (1998). Models using packing spheres represent the porous medium as a geometric system that can be described by tetrahedral subunits formed by groups of four touching spheres. The probability of occurrence of each subunit is expressed in terms of the volume of the spheres, their probability of being present and the probability that specific pairs of spheres in contact may exist. The pore space in a subunit can be approximated by a narrow throat in each of the four faces of the tetrahedral subunit that meet at a pore body, which occupies the central cavity of the tetrahedron. Fig. 2 shows an example of a tetrahedral subunit formed by four touching spheres and an approximation of the pore space structure. The size of the pore throat is evaluated by a cylinder with a diameter given by the sphere in between the three spheres in the face. The pore body-size is determined by the sphere created in between the four touching spheres. Formulas given in the Rouault and Assouline (1998) work for computing pore body size (r_{p}), pore throat size (r_{ci}) and their relative probability densities were used.
The input data for sphere packing are in the form of a discrete sphere-size distribution, which can be derived from the GSD of the porous medium. We describe the GSD of a sand used in our laboratory experiment by a normal distribution, which is commonly used for an analytical description of grain-size distribution and can easily be applied in practice. The density probability function f(d) of this distribution is given by:
$$f\left(d\right)=\frac{1}{\sqrt{2\pi {\sigma}_{\text{g}}}}\mathrm{exp}\left(-\frac{1}{2}{\left(\frac{d-{d}_{50}}{{\sigma}_{\text{g}}}\right)}^{2}\right)$$ | (12) |
For each set of four radii, the probability density of the tetrahedron subunits and the corresponding radius of the “kissing” sphere are obtained. The final result is one pore-size distribution curve and four throat-size distribution curves. Assuming that the generated pore body-sizes and pore throat-sizes are normally distributed, which is generally assumed for a large sampling number, the assessment of the mean value and the standard deviation of corresponding curves allows for calculation of various ranges of pore body radii and pore throat radii. These parameters are geometric input data of the pore-throat model.
3.2 Application to a medium-size sand
In this section, we use the conceptual approach described in Section 3.1 to parameterise the equivalent network of the H2F sand that was used in the physical experiment. H2F is a medium-size quartz sand with a hydraulic conductivity of about 8 × 10^{−4} m/s (k = 8 × 10^{−11} m^{2}), a total porosity of 40% and a low fraction of organic content (foc = 0.09%, based on NFT 31–109) (Bohy et al., 2006).
Fig. 3 shows the grain-size distribution curve of the H2F sand used in the experiment and its corresponding derived sphere-size distribution considered in the model of sphere packing. The mean value and the standard deviation of the sphere-size distribution allow positive values at each distance even at three standard deviations away from the mean value. Twenty-one spheres of different sizes were considered in the packing. This relatively high number of spheres ensured the reliability of statistics on the evaluation parameters as C^{4}_{21} (5985) possible sets of four radii might be formed. The radii of chosen spheres are equally spaced; this means that the radii are arithmetically related. The sizes of the radii vary between a minimum value of 0.05 mm and a maximum value of 0.34 mm. The ratio between the maximum and minimum sizes satisfies the packing condition mentioned in the work of Dodds (1980).
The generated pore body and pore throat-size distributions are shown in Fig. 4. In contrast to the pore body-size distribution, the packing model allows us to compute four different distributions that can be used to describe the pore-throat size. The pore body-size distribution has a shape similar to the sphere-size distribution. The results agree well with the findings of Rouault and Assouline (1998). The forms of throat-size distributions vary and occupy a large range of radii. Based on the assumption that the volume of pore throats is negligible compared to that of pore bodies, we limited our analysis to the pore-throat distribution that presents the smallest mean radius (r_{c1}).
Table 1 shows various values derived for both pore-body radii and pore-throat radii using the mean value and the standard deviation of the corresponding generated size distributions. For both size distributions, Case 1, Case 2 and Case 3 refer to the minimum and maximum radii derived from a size distribution as one time, two times and three times the standard deviation of the size distribution, respectively. For fixed distribution parameters of pore body radius and by combing the pore body radii with the pore throat radii given in Case 1, Case 2 and Case 3, we obtained at least nine possible realisations of the pore space geometry of our network model. To reduce the number of possible realisations of the geometry of the network model, we used the macroscopic properties of medium-size H2F sand, such as the intrinsic permeability (k = 8 × 10^{−11} m^{2}) and porosity (ø = 40%). For each of the nine possible combinations, we fixed the distance between pore-body centres at L = 0.14 mm, which is greater than two times the maximum radii of pore bodies. A single-phase water flow was simulated using the generated network model with given pressure boundaries. The intrinsic permeability and the porosity of the network, calculated with Eqs. (2) and (11), respectively, were then compared to the measured values of the sand used in the laboratory experiment (Table 2). The combination of Case 1 of the pore-body radii with Case 2 of the pore throat radii gave the best fit for the macroscopic properties of the sand. While the calculated permeability of 8.6 × 10^{−11} m^{2} was quite close to the experimental value, the porosity of our network model (ø = 47%) overestimated the porosity of our sand. The porosity might be better estimated by increasing the distance L between pore body centres. However, this would cause a decrease in the equivalent macroscopic permeability.
Rayons minimal, moyen et maximal de pores et de capillaires, calculés par le modèle d’empilement de sphères.
Pore body radius | Pore throat radius | |||||
r_{min} (mm) | r_{mean} (mm) | r_{max} (mm) | r_{min} (mm) | r_{mean} (mm) | r_{max} (mm) | |
Case 1 | 0.042 | 0.048 | 0.057 | 0.026 | 0.03 | 0.036 |
Case 2 | 0.036 | 0.048 | 0.061 | 0.022 | 0.03 | 0.04 |
Case 3 | 0.03 | 0.048 | 0.067 | 0.017 | 0.03 | 0.046 |
Propriétés macroscopiques du réseau de pores et de capillaires (calculées pour neuf réalisations), comparées à la perméabilité intrinsèque et à la porosité du sable moyen H2F.
Pore body radius | Pore throat radius | Macroscopic proprieties | |||||
r_{min} (mm) | r_{mean} (mm) | r_{max} (mm) | r_{min} (mm) | r_{mean} (mm) | r_{max} (mm) | ø (%) | K (10^{−11} m^{2}) |
0.042 | 0.048 | 0.057 | 0.026 | 0.03 | 0.036 | 47 | 12.87 |
0.042 | 0.048 | 0.057 | 0.022 | 0.03 | 0.040 | 47 | 8.6 |
0.042 | 0.048 | 0.057 | 0.017 | 0.03 | 0.046 | 47 | 13.2 |
0.036 | 0.048 | 0.061 | 0.026 | 0.03 | 0.036 | 49 | 12.84 |
0.036 | 0.048 | 0.061 | 0.022 | 0.03 | 0.040 | 49 | 7.64 |
0.36 | 0.048 | 0.061 | 0.017 | 0.03 | 0.046 | 49 | 13.42 |
0.30 | 0.048 | 0.067 | 0.026 | 0.03 | 0.036 | 53 | 11.73 |
0.30 | 0.048 | 0.067 | 0.022 | 0.03 | 0.040 | 53 | 11.83 |
0.30 | 0.048 | 0.067 | 0.017 | 0.03 | 0.046 | 53 | 12.05 |
4 Numerical simulation of a laboratory experiment
4.1 Model parameters
One of the objectives of our study was to use experimental results (Nsir, 2009) for verification of the network model developed. Experiments were conducted in a sand-filled glass column with a length of 70 cm and inner diameter of 10 cm. The DNAPL used was TCE, which was injected at a constant flow rate of 40 mL/min. The displacement was performed from bottom to top, which reproduced a stable displacement condition, and from top to bottom, which reproduced an unstable gravity condition. The experimental device was equipped with pressure transducers for monitoring pressure at the inlet and outlet sections of the glass column. The experimental programme also allowed the measurement of the arrival time of the DNAPL front at eight points on a control section equipped with optical fibre sensors (Nsir, 2009).
Simulation results reported here were carried out on a network of 7×7×4858 pore bodies. The spatial discretisation used reproduces the total height of the physical model but only a small fraction of the column section. We assumed that the chosen network is sufficiently large to be a representative elementary volume (REV) which may reproduce the macroscopic properties of the sand column such as intrinsic permeability and porosity. The pore-body and pore-throat radii considered in the simulations correspond to Case 1 for pore bodies and Case 2 for pore throats, as given in Table 1. The distance between pore-body centres was fixed at L = 0.14 mm.
The network model is initially water-saturated, and a constant flow rate of 29.6 × 10^{−4} mL/min of TCE was applied at the inlet section of the model. The injection rate was calculated by accounting for the ratio between the modelled section and the total section of the glass column. At the outlet section of the model, the observed water head was specified as a constant pressure boundary. The simulation stopped when the duration of the DNAPL injection that occurred in the experiment was reached.
4.2 Stable case of water displacement by DNAPL
To illustrate the case of stable displacement, the displacing fluid (TCE) was injected at the bottom of the network. The boundary pressure applied to the outlet section corresponded to a water height of 18 cm. Fig. 5a shows a comparison of calculated and measured inlet pressures as a function of time. According to the experimental observations, the value of the pressure in the whole domain increased with time. The water/DNAPL interface was stabilized by the gravity effect; the potentially destabilizing influence of the viscosity ratio on the front behaviour was not significant. The increase in pressure was therefore caused by the high density of the invading fluid, which overrides the decrease of viscous pressure during the displacement of the less viscous DNAPL.
Even though the increasing pressure behaviour observed in the experiment was adequately reproduced, the measured pressure increases were stronger than what was calculated experimentally. This can be explained by the chosen arrangement of pore and throat sizes in the pore-throat model. The very narrow throats that exist in the real porous medium were not considered in the network and may require a higher displacement pressure to be invaded by the displacing fluid. Using a wider pore body-size distribution might also reduce this misfit of inlet pressure.
The curves representing the distributions of dimensionless arrival times of the water/DNAPL front in the experiment and the simulation were similar (Fig. 5b). The calculated standard deviations were very close to the measured values (14 s calculated compared to 9 s measured in the experiment). According to our experimental observations, the distribution of front arrival times at the control section is almost uniform. Due to the stabilizing effect of gravity, a flat piston-like front was observed moving through the network. Consequently, very little wetting fluid was left behind and a compact pattern of the invading fluid occurred. The network model adequately reproduced the stability mechanisms of displacement caused by the positive contribution of gravity.
4.3 Unstable case of water displacement by DNAPL
In the case of unstable water displacement, the DNAPL was injected at the top of the network using the same boundary and initial conditions as described in section 4.1. For this displacement, we expected a much more irregular front with possible gravity fingering. Fig. 6a shows the variation of displacement pressure measured and calculated by the model. Due to the constant injection rate, the pressure value in the domain decreased as the less viscous and denser fluid invaded the system. The sudden increase of pressure shown in Fig. 6a corresponds to the breakthrough of the invading fluid at the outlet section. Under these displacement conditions, gravity forces are destabilizing and dominate the movement, resulting in typical DNAPL fingers in the displaced fluid. The macroscopic water/DNAPL interface was no longer horizontal and many throats were blocked by the capillary effect. The number of active invasion paths was therefore low, and the flow section of DNAPL was reduced. The contribution of viscous pressure was also limited and was dominated by the contribution of the higher density of the invading fluid. This resulted in a decrease of the inlet pressure. The inlet pressure observed was reproduced by the numerical model, but the rate at which it decreased in the model was different from the measured one, the calculated pressure decrease was more substantial than that observed. In addition, we recorded an earlier breakthrough time than we observed. This result can be explained by the number of active invasion paths that were lower in the network model than in the real porous medium. This effect may be due to the narrow generated pore-throat size distribution that contains many large pore throats.
Fig. 6b shows a comparison between the distribution of calculated arrival times of the front and those observed at the control section. The statistical data of both distributions are similar. They both illustrate highly non-uniform distributions of arrival times with significant standard deviations. The calculated standard deviation was 75 s, compared to 70 s measured in the laboratory experiment. Thus, gravity fingering was numerically confirmed by the distribution of arrival times of the front at the chosen cross-section of the network. However, the mean arrival time of the water/DNAPL front calculated at the control section was approximately 633 s, compared to the observed time of approximately 747 s. This finding directly confirms the earlier computed breakthrough of TCE at the downstream section (Fig. 6a).
4.4 Sensitivity study
A sensitivity study of specific model parameters, such as the size of the lattice and the size of pore throats, was undertaken to analyse the influence of these parameters on the displacement of water by DNAPL. Furthermore, we assessed whether a reasonable modification of their values might explain the differences between the simulated and measured DNAPL pressures at the inlet section, as well as the simulated and measured arrival times of the DNAPL front. Specifically, we analysed the unstable case of water displacement by DNAPL.
4.4.1 Effect of the lattice size
Using a larger section of network may be more valid for the appropriate representation of the real porous medium. However, larger networks suffer from computation times that can become, in some cases, very high, even with powerful computers. Fig. 7a shows the calculated inlet pressure obtained from simulations carried out on networks of two lattice sizes (5×5×4858 and 6×6×4858 pore bodies) that were smaller than the reference case and one lattice size (9×9×4858) that was larger than the reference case discussed in Section 4.3. The size of the reference case was 7×7×4858. On a Linux-computer with 3956 MB RAM and a Celeron processor running at 1998 MHz, the corresponding Central Processing Unit (CPU) times were very different; the times were one day, four days, thirteen days and twenty days for the two small lattices, the medium lattice and the large lattice, respectively. In all cases, the modelled inlet pressure behaviour was similar; it decreased with increasing time, particularly for the reference case and the case of the bigger lattice (the pressure profiles for those cases were almost the same). However, the breakthrough times calculated for the bigger lattice and the reference lattice were shorter than the times calculated for the two smaller lattice sizes. The smaller lattices were closer to the observed breakthrough time. Thus, developed displacement instabilities are more significant in bigger lattices than in smaller lattices. The arrival time distribution of the front water/DNAPL calculated and measured (Fig. 7b) in the same control section confirms this conclusion. The corresponding curves for both bigger lattices were almost the same and were more widely distributed than those computed using the smaller lattices. The standard deviation of the distribution of arrival times was 75 s and 79 s for lattice sizes 7×7×4858 and 9×9×4858, respectively. These values are very close to the measured standard deviation of 72 s. The standard deviation arrival times for the smaller lattices sizes (5×5×4858 and 6×6×4858) were 59 s and 64 s, respectively. A large lattice size leads to a spatial pore body size distribution and pore throat size distribution that is more heterogeneous than for a small lattice. Consequently, preferential pathways are more likely to occur during displacement. As a result, the water/DNAPL front arrives earlier at the downstream section of the model in the case of bigger lattice sizes. Thus, the distribution of arrival times of the water/DNAPL front is not uniform, as shown in Fig. 7b. In contrast, based on the numerical results obtained from the 7×7×4858 and 9×9×4858 lattices, which were quite similar, we concluded that a lattice size of 7×7×4858 adequately reproduces the unstable displacement process in the given porous medium.
4.4.2 Effect of the variance in pore-throat size
In this study, we used the three distribution parameters for pore throat radius as given in Table 1 (Case 1, Case 2 and Case 3) to analyse the influence of the variance in the pore throat size on the unstable displacement of water by DNAPL. Case 1 has a pore throat size variation range smaller than Case 2, which is hereafter referred to as the reference case. Case 3 has pore throat size variation range larger than the reference case. The size of the network considered in our simulations was the same for each pore throat size distribution (with identical pore body size distribution). We chose the smallest lattice size (5×5×4858) because it required less computational time and memory.
Fig. 8a shows the calculated pressure for different network constructions that were compared to the measured pressure. We note that the inlet pressure increased with increasing variance of the pore throat distribution and was accompanied by a later breakthrough of the water/DNAPL front. Apparently, the pore-throat radius controls the dynamic displacement of the water/DNAPL interface in the network. Case 3 had the smallest pore-throat sizes. This resulted in an increase in pressure in the pore bodies at the upstream boundary of the model, which was required to maintain the prescribed constant mass flow rate for DNAPL. A higher local pressure field permitted the DNAPL to overcome the capillary entry pressure in the throats and therefore resulted in a more uniform invasion of the water-saturated pore bodies by DNAPL in the model transect. However, the pore-throat size distributions of Case 1 and Case 2 had much larger throats, resulting in preferential pathways that were more likely to occur. This made the overall invasion of the section harder. As a consequence, the model predicted earlier breakthrough times, and the pressure field decreased rapidly. This may explain why the computed breakthrough of the displacement front in Case 3 was later than in Cases 1 and 2. Furthermore, in analysing the calculated DNAPL saturation field (not presented here), the DNAPL front in Case 3 was characterised by several gravity-dominated fingers, whereas in Case 1, fewer fingers were formed.
Fig. 8b shows a comparison between distributions of dimensionless arrival times of the water/DNAPL front calculated and measured at the control section. The prediction with a pore-throat size distribution as in Case 3 was quite close to the experimental data. The corresponding curve was more widely spread than those obtained with the narrow pore-throat size distributions of Case 1 and Case 2. The calculated standard deviation relative to the pore-throat distribution of Case 3 was approximately 69 s, which is quite close to the experimentally determined value of 72 s. However, the calculated standard deviation was only 49 s and 59 s for Case 1 and the reference case, respectively. Case 3 provides a wider pore-throat radius distribution. This may allow much invasion pathways to be formed, leading to a better reproduction of the displacement front observed in the laboratory experiment. Because the pore-throat size distribution variability was lower in Case 1 and Case 2, one may think that the number of accessible pore throats for the displacing fluid was reduced and therefore only a few fingers resulted.
5 Conclusions
We developed a pore-throat network model that adequately describes the dynamic drainage process. A geometric model of packing spheres combined with a probabilistic approach was used to obtain realistic pore sizes and throat sizes based on grain-size distribution. The pore-scale network model was tested using experimental data. The prediction of the temporal evolution of pressure at the inlet section and the distribution of arrival times of the water/DNAPL front in a given cross-section matched well with laboratory measurements. An important feature of the model is its capability to reproduce the observed pressure behaviour for stable and unstable displacement regimes. The sensitivity study indicated that the lower pore-throat size strongly influenced the instability of DNAPL displacement. However, additional numerical studies may be necessary to analyse the influence of pore body size on unstable DNAPL displacement. Furthermore, application of the developed numerical model to a drainage experiment recently conducted in a homogeneous fine sand may allow an additional verification step of the model to derive pore body sizes and pore throat sizes in the case of a low permeability porous medium.
The network model used in the present work considers the converging–diverging nature of pores and is based on the pore body size distribution and pore-throat size distribution. However, no correlation exists between pore throat sizes and neighbouring pore bodies. Besides, pore space asperities and roughness were not considered. Additional quantitative information about the degree of correlation between pore-body size and pore-throat size are needed to take into account more features of the pore space geometry of the porous medium. The lattice sizes used resulted in time consuming simulations. This prohibits application of the model to a larger part of the cross-section of the porous medium. One possible way to overcome this problem might be to develop a more efficient algorithm to solve the large set of equations.
Acknowledgements
Financial support for this research was received from the programme REALISE (REseau Alsace de Laboratoires en Ingénierie et Sciences pour l’Environnement). The Région Alsace, the GDR “Hydrodynamique et Transferts dans les Hydrosystèmes Souterrains” (INSU-CNRS), and the Conseil Scientifique de l’Université Louis-Pasteur Strasbourg are gratefully acknowledged. Furthermore, the authors would like to thank Professor Insa Neuweiler and the two other anonymous reviewers for their valuable comments and suggestions, which helped to improve the article significantly.