Capillary trapping mechanism in strongly water wet systems: Comparison between experiment and percolation theory

Advances in Water Resources 79 (2015) 35–50

Contents lists available at ScienceDirect

Advances in Water Resources journal homepage: www.elsevier.com/locate/advwatres

Capillary trapping mechanism in strongly water wet systems: Comparison between experiment and percolation theory Helmut Geistlinger ⇑, Sadjad Mohammadian Helmholtz-Centre for Environmental Research – UFZ Leipzig-Halle, Department of Soil Physics, Theodor-Lieser-Str. 4, 06120 Halle/Saale, Germany

a r t i c l e

i n f o

Article history: Received 6 August 2014 Received in revised form 13 February 2015 Accepted 16 February 2015 Available online 25 February 2015 Keywords: Capillary trapping Thin-film flow Trapped gas cluster X-ray micro computer tomography Monte-Carlo experiment Percolation theory

a b s t r a c t To understand capillary trapping mechanism, we conduct a real Monte-Carlo experiment by using packed glass beads with nearly the same pore size distribution, but different stochastic realizations. We study gas phase trapping during imbibition for capillary numbers from 2  107 to 106 using X-ray micro tomography and compare the experimental results with predictions from percolation theory. We found excellent agreement. Percolation theory explains (i) that the capillary desaturation curves are not dependent on flow rate, (ii) the linear dependence of the total gas surface on gas saturation that is a direct consequence of the linear relationship between cluster surface area and cluster volume, which is a prediction from percolation theory for large finite clusters, (iii) the power-like cluster size distribution with an exponent sexp = 2.15 that only deviates by 2% from the theoretical one (stheor = 2.19), and (iv) that the maximal z-extension of trapped large gas cluster is described by the cut-off correlation length nB (B – bond number). Ó 2015 Elsevier Ltd. All rights reserved.

1. Introduction Capillary trapping of gas bubbles, residual NAPL, and oil blobs within water-saturated porous media and mobilization of such trapped, isolated phases is of central importance in many processes in oil recovery, hydrogeology and soil physics. For CCS-technology (CCS – carbon capture and storage) capillary trapping is one of the relevant storage processes [1,2]. To achieve homogeneous distribution of trapped oxygen gas bubbles and/or mobilization of trapped residual NAPL-blobs is still a challenge in remediation of contaminated groundwater [3–6]. The effectiveness of oil recovery will depend on how much isolated oil blobs remain in the porous matrix after water flooding [7,8]. Air entrapment during rain infiltration and groundwater level fluctuation is the key process in the unsaturated zone that determines the gas exchange, redox status, and the aerobic microbial activity of the unsaturated soil zone, especially in the highly transient zone – the capillary fringe [9]. Recently, we studied the capillary trapping mechanism of gas bubbles within 1 mm-glass beads by X-ray micro-computer-tomography experiments (l-CT) and found no systematic dependency of trapping efficiency on capillary number Ca (=vlw/r, v – mean velocity of the injected fluid, lw – viscosity of water, r – interfacial tension) between 2  107 and 106 [10]. To understand this behavior, we reviewed the relevant literature on imbibition ⇑ Corresponding author. http://dx.doi.org/10.1016/j.advwatres.2015.02.010 0309-1708/Ó 2015 Elsevier Ltd. All rights reserved.

process for strongly-water wet systems and for small flow rates, i.e. for capillary numbers <106. Since the pioneering work by Lenormand and Zarcone [11,12] the pore-level mechanism of capillary trapping for imbibition was studied intensively by Vizika et al. [13], by Blunt and Scher [14], by Hashemi et al. [15,16], by Constantinides and Payatakes [17], and by Joekar-Niasar et al. [18,19], and by Hilpert et al. [20]; for reviews, see [21–23]. According to experimental observations [12,17] capillary trapping can be caused by different displacement mechanisms: (i) pistonlike displacement with cooperative pore-filling and (ii) by precursor thin-film and corner flow. Sometimes the corner flow is also called thin-film flow in the literature [15,16,23] or thin-film flow and corner flow are summarized as flow in crevices [14]. If not explicitly distinguished, we will use the term thin film flow both for flow along the grain surface and along the corners of the pore space. Which displacement type dominates is strongly dependent on wettability, on flow rate, on pore-throat-geometry, and on pore-network topology (connectivity) [11,12,14]. The temporal sequence of displacement types in a strong-water wet system ‘‘glass (SiO2)–water–air’’ is governed by the order of the corresponding capillary pressures (see Table 1 in [11]). First, a precursor thin-film in front of the mean wetting front (bulk advance) will ‘‘spontaneously’’ invade the porous media and cover some region. The thickness of such thin water films is often estimated by a typical surface roughness of the order of about

36

H. Geistlinger, S. Mohammadian / Advances in Water Resources 79 (2015) 35–50

Table 1 Different imbibition mechanisms and its relation to percolation theory proposed by Lenormand and Zarcone [11]. Imbibition mechanism

Percolation type Ca = 109 large pores

109 < Ca < 106 medium-sized pores

Thin-film flow/flow by surface roughness

Bond percolation (ordinary percolation)

Flow by corners/crevices

Invasion percolation

Bond percolation + cooperative I1  pore filling Invasion percolation + cooperative I1  pore filling

10 lm [11,14,17]. However, a more detailed study by Kibbey [24] shows that depending on the surface roughness (sand grains and glass beads were compared) and the surrounding pressure of the nonwetting phase the film thickness can vary from nm- to lmscale. The spreading of this precursor thin film is always faster than the subsequent bulk advance caused by pistonlike displacement process. Under certain unstable conditions this film can swell and snap-off in the smallest throats in front of the mean displacement front without any ‘‘visible’’ continuity to the bulk wetting phase. This results in randomly filled throats (bonds) in order of increasing throat diameter and corresponds to classical bond percolation (ordinary percolation). In case where corner flow is the dominant mechanism, the interface moves by a succession of jumps. All the menisci in pores and along corners are at the same pressure which progressively decreases as the wetting fluid is injected (the volume of fluid increases in the corners, therefore the radius of curvature increases). The jump occurs when the pressure reaches the critical snap-off threshold in the smallest duct at the front. The radius of the swelling corner film starts to decrease once the critical capillary pressure for snap-off (Eq. (3d)) is reached and the gas–water interface becomes unstable. If the interface moves by such snap-off-filling events, this mechanism is similar to an invasion percolation mechanism. As the flow rate increases, but Ca is still smaller than 106, cooperative I1-pore filling mechanism (I1: cooperative pore filling with wetting fluid by Z  1 throats, Z – coordination number) also occurs, but the snap-off trapping mechanism remains the dominant mechanism. Possible imbibition mechanisms for the relevant range of capillary numbers are shown in Table 1 that is a section of Table 3 from [11]. Based on the work by Lenormand and Zarcone [11], Blunt and Scher [14] have developed a 3D-pore-level model of wetting that represents these effects (thin-film- and corner flow, cooperative pore filling mechanisms) and discussed possible trapping mechanisms and their dependence on flow rate. We will compare our experimental results with their theoretical results. With the noninvasive technique of l-CT and new sophisticated algorithms of image processing [25] the data sets are becoming highly reliable and amenable to cluster analysis at the pore scale. In order to make clear what are the new contributions of our paper, we give a short review of l-CT studies on non-wetting phase trapping during imbibition, which focus on cluster analysis. Iglauer et al. [26] analyzed the residual oil cluster size distribution and find consistency with percolation theory that predicts a power-law cluster size distribution, n(s)  ss, with s = 2.189 [27], where s denotes the number of sites (=pores) of a trapped cluster. The authors obtained a power-law exponent for the oil-wet core s = 2.12 (oil–brine–sandstone system, Ca = 107)

which is higher than s for the water-wet system (s = 2.05, oil– brine–sandstone system, Ca = 2  106 [28]) or a CO2/brine system (s = 2.01 scCO2–brine–sandstone system, Ca = 106 [2]). The derived s-values are in reasonable agreement with the prediction from percolation theory. However, fitting the power law to the data, the voxel number was used instead of the site number, i.e. no renormalization of the n(s)-distribution from voxels to pores was carried out. Furthermore, Iglauer et al. [8] analyzed the relationship between the surface area Ai and the volume Vi of individual clusters. They found a power law relationship, Ai = const (Vi)p with p  0.8 based on a best fit using a log–log-plot. An overview of derived s- and p-values for different multiphase systems is given in [8]. Landry et al. [29] (Ca  106) studied the microscopic Ai–Virelationship and found a similar p-value of about 0.8 based on a log–log plot. Furthermore, they studied the macroscopic relationship between the total surface area (Anw) and the total volume (Vnw) of non-wetting phase and found a strong linear relationship (regression coefficient = 0.98). Georgiadis et al. [30] conducted imbibition and drainage experiments with brine–oil fluid pairs varying interface tension from 9 to 51 mN/m. The capillary numbers were in the range between 105 and 104. Assuming that percolation theory is still valid for these high capillary numbers (compare with Table 1), they fitted their experimental data to the universal power law and could not achieve reasonable agreement, especially in the range of large cluster sizes. Since the universal power law is valid in the limit of large cluster sizes [31], the strong deviations indicate that percolation theory is not applicable and that the data set may be incomplete in this range, because of the too small sample dimensions (10 mm  20 mm). Al-Raoush and Willson [32] studied the imbibition process in 0.4–0.6 mm oil-wet glass beads for a brine–oil fluid pair for Ca = 2  106. The focus of this study was to derive a pore-throat network model with realistic pore- and throat distributions. They obtained a mean pore radius of 58 lm and mean throat radius of 39 lm and a mean coordination number of about 4. They conducted a cluster analysis for the residual brine clusters and found similar results as Geistlinger et al. [10], i.e. the majority of trapped clusters (75%) had a volume less than the average grain volume. About 8% of the blobs are branched, complex in shape, and extend over more than 10 pores. Interestingly is that the cumulative cluster size distribution shows a natural cutoff-length at about five mean grain diameters (=d50). For the individual clusters the authors derived also a shape factor (deviation from spherical shape) and the cluster orientation along the major axis. The motivation for this paper was initiated from a recently conducted experimental l-CT-study on capillary trapping mechanisms of gas bubbles within 1 mm-glass beads. We found the following interesting results [10]: (i) no systematic dependency of trapping efficiency on capillary numbers between 2  107 and 106, (ii) the majority of trapped gas bubbles (about 85%) are multipore trapped gas clusters, and (iii) a significant (R2 = 0.98) linear relationship between the total gas surface area and gas saturation. The objective of this paper is to understand these results based on percolation theory, because this is ’’. . .the natural language for describing connectivity effects. . .’’ of fluid phases during displacement of one fluid by another fluid [33]. We focus on the imbibition displacement process for strongly-water wet systems. Our paper is organized as follows: In Section 2 we give a short overview of the l-CT experiments. In Section 3 an analytical study of the fluid–fluid displacement process is carried out, and in Section 4 we compare the experimental results with predictions from percolation theory.

H. Geistlinger, S. Mohammadian / Advances in Water Resources 79 (2015) 35–50

2. Experiments 2.1. Experimental methodology A column of 30 cm length and an inner diameter of 3.18 cm was filled with 1 mm-glass beads (d50 = 0.85 mm, d10 = 0.8 mm, d95 = 1 mm). These experiments are called first series of experiments. In a second series of experiments we used a thinner column with inner diameter of 2.16 cm and total length of 20.0 cm, in order to achieve a higher image resolution. Four classes of glass beads (0.5 mm, 1 mm, 2 mm and 0.5 mm–1 mm-mixed) were used in these experiments. Details of grain size distribution and packing parameters are presented in [9]. Air was used as non-wetting fluid and air-saturated water as wetting fluid. Exemplarily, we will discuss here the experimental results of the first series of experiments. Similar results were obtained for the second series of experiments [9]. Prior to packing, the glass beads were cleaned with acetone and water and dried at 105 °C. Iglauer et al. [34] studied the impact of chemical cleaning agents on surface wettability and contact angle. Because of the high wettability of these cleaning agents, they have to be used carefully, in order to prevent surface contamination and to obtain reproducible initial conditions. To minimize surface contamination by acetone cleaning, we subsequently used an ultrasonic water bath for 2 h. The good reproducibility of our experiments indicates that the contact angle varies only in a small value range. To achieve a close, reproducible packing, small vibrations to the column were applied during the continuous filling process. All packed columns exhibit close random packing resulting in a porosity of 0.332 ± 0.005 (packing density = 1.61 kg/L, particle density = 2.41 kg/L). Before saturating the column with water the column was flushed with CO2–gas to ensure that any trapped gas phase will dissolve rapidly. The column was then saturated with water to a certain position z0 that was measured with a by-pass capillary (diameter = 5 mm). Since we have increased the water table within the column nearly quasi-statically, we assume that no gas bubbles are trapped and the initial water saturation below z0 is zero. To prove this, we evaluate the volumetric gas content by l-CT and found a value smaller than 1%. The initial state was prepared prior each new experiment. Starting from the initial state we increased the water table to z1 = z0 + 5 cm with a certain water table (WT) rising velocity of 0.1, 0.2, 0.3, 0.4, 0.5, and 0.6 cm/min. These rising velocities corresponds to capillary numbers between 2  107 and 106. We note that the capillary numbers are comparable to scCO2-injection into brine saturated sandstone cores (106) [2] and to brine flooding of oil–gas-saturated sandstone cores (106) [8]. To estimate the influence of a varying local pore structure and connectivity we conducted three experiments. Each was repacked prior to the experiment and then the same WT rising velocity of 0.5 cm/min was applied. Subsequently, the static trapped gas phase was recorded by l-CT to measure the gas phase distribution. To ensure that the trapped gas phase is static, i.e. no bubble movement and no dissolution into the air-saturated water phase occur; we conduct a long-term experiment and measured the volumetric gas content and the total gas surface area with l-CT after 20 min, 3 h, and 24 h. Both the total gas surface area and the volumetric gas constant are constant during the experiment and the relative standard deviation (SDV) < 2% and <1%, respectively. The constant total gas surface area indicates no gas bubble movement, since a new bubble configuration will exhibit in general another total gas surface area [10]. The device used in this study is an industrial X-ray scanner with a focal spot of 5 lm (X-Tek HMX225). Objects are rotated in a cone-beam (135 keV, 633 lA), and projections are recorded on a

37

CCD-Detector Panel with a resolution of 1024  512 pixels. Reconstruction of the data was done with the X-Tek CTPro software package. Important for our experiments is the spatial resolution, since it determines to which lower spatial threshold we can detect trapped gas bubbles. From simple regular cubic and face-centered packing one can estimate the following mean pore diameter dk,max = 0.56 mm [35]. The resolution of the CT-image, i.e. the voxel size has to be somewhat smaller than this diameter that is assumed to reflect the lower threshold of stable gas bubbles. In fact, the nominal resolution of our CT-images was 0.084 mm. In this case a trapped gas bubble with a diameter of 0.56 mm can be represented by about 157 voxels, i.e. 6–7 voxel along the diameter. 2.2. Experimental results 2.2.1. Pore-size distribution Fig. 1 shows the pore-size distribution (PSD) of nine different re-packed 1 mm-glass beads sediments (GBS) used for the first series of experiments. The continuous PSDs are best fits to the experimental data. The experimental PSD was obtained from the segmented CT-image using a standard method based on mathematical morphology, i.e. the opening size distribution [36]. This algorithm determines the maximum diameter of a sphere, which can be placed inside the pore space at any pore voxel position. This leads to the volume Vi of different pore size classes di and the corresponding probability distribution of sphere diameters Vi/ porosity representing the pore size distribution of the structure. The pore structure was analyzed for 3–5 different pore size classes and their relative frequencies listed in Table S1 (Supplementary material S1). We fitted two different continuous probability density functions to the relative frequencies: (i) a Normal distribution and (ii) a LogNormal distribution. Since both functions are 2parameter functions (p1, p2) determined by its mean and variance, we have 3–5 constraints for two equations. Therefore, one has to solve an optimization problem [10]. Table 2 compares the relative errors (difference between theoretical and experimental relative frequency divided by experimental frequency; columns 4 and 5) for both continuous PSDs. The respective mean relative errors are 0.37% for the Normal distribution and 14.41% for the LogNormal distribution. Hence, the Normal distribution yields a quasi-exact fit, i.e. the resulting probabilities of Normal distribution are equal to the experimental relative frequencies. We note that all PSDs are nearly identical with a mean pore size of 0.17 ± 0.005 mm, i.e. the SDV is about 3% of the mean value (see Table 2, columns 2 and 3). All three repeated experiments (7, 8, 9) exhibit the same PSD representing the mean values: mean pore radius rp = 0.17 mm and mean variance r2 = 0.36 mm2. 2.2.2. CT-images and image processing Fig. 2 shows two different sections of a reconstructed 3D-image after a WT-rising at a rate of 0.5 cm/min (experiment 7). In both images the trapped gas phase is colored red. Fig. 2(A) is a 2D-cross section showing all 3 phases (grey: glass beads, black: water), whereas Fig. 2(B) is a 3D-slice of the reconstructed image showing only the gas phase and the solid phase . The gas bubbles or gas clusters (both terms are used synonymously) are isotropic homogenously distributed, i.e. no z-correlation can be seen. The trapped gas bubbles are of different sizes; e.g. small bubbles that are trapped within single macro-pores and large gas bubbles or gas clusters that span over few grain diameters. The existence of such large gas clusters (like oil- or NAPL-ganglions) was surprising and an unexpected result, because the standard assumption in literature is that large gas clusters become unstable and split into smaller ones, due to their high surface tension, strong water

38

H. Geistlinger, S. Mohammadian / Advances in Water Resources 79 (2015) 35–50

Fig. 1. Pore-size distributions of nine different re-packed 1 mm-glass beads. Shown are the best fit curves (Normal distribution) to the experimental data.

Table 2 Best-fit parameters of continuous PSDs for nine different glass-beads packings. Experiment

Normal distribution

Relative errors in%

Pore radius [mm]

Sigma [mm]

Normal

LogNormal

1 2 3 4 5 6 7 8 9 Mean SDV

0.17 0.17 0.18 0.18 0.18 0.18 0.17 0.17 0.17 0.17 0.005

0.06 0.06 0.06 0.06 0.07 0.06 0.06 0.06 0.06 0.06 0.002

0.23 0.32 0.21 0.18 0.11 0.20 0.74 0.77 0.59 0.37 0.26

20.73 12.12 21.98 9.13 15.44 8.97 14.33 14.44 12.53 14.41 4.54

Fig. 6 in [17]). We present in Fig. 3 the normalized CDCs for both series of experiments and found no systematic dependence. The plateau values (=mean gas saturations) for the different glass beads are: 0.1 for the 0.5 mm-GBS; 0.056 for the 1 mm-GBS-old; 0.041 for the 1 mm-GBS-new; 0.09 for the 2 mm-GBS-new; and 0.087 for the mixed GBS [9]. The normalization is just done by dividing the actual gas saturation by its mean value. For the 1 mm-GBS we show the results of both series of experiments (first series: 1 mm-old, filled rhombuses; second series: 1 mm-new, empty squares). 2.2.4. Linear relationship between nonwetting phase surface area and phase saturation The values of the specific gas phase surface range from 0.149 to 0.296 1/mm for gas saturations ranging from 3.8% to 7.9% with a mean value of 5.6% (first series). In order to compare our data with the data presented by Landry et al. [29] for trapped oil phases in brine saturated glass beads (0.918–1.226 1/mm for oil saturations ranging from 16.7% to 21.9%), we transformed the saturation of the trapped nonwetting (nw) phase to a reduced nw-phase saturation: Snw,red = (Snw  Snw,min)/(Snw,max  Snw,min) and normalized the surface area of the nw-phase by its maximum value. Fig. 4 shows the dependence of the total gas surface area on gas saturation. Both data sets show a strong linear dependence (R2 = 0.98) of the nw-phase surface area on gas saturation (diamonds: first series of experiments; circles: Landry et al. [29]-data, first imbibition of brine into oil-saturated glass beads). This linear dependence was also found by Pan et al. [41] for residual NAPL (non aqueous phase liquid; Toluene).

wettability, and the large pore-throat aspect ratio of glass beads packings [10]. Image processing and phase segmentation was conducted for quantitative evaluation of the trapped gas phase. The experimental porosity was used as a constraint for testing the segmentation algorithm. A short overview of the procedure is presented here (for details see [10]). First, images were de-noised using Total Variation method introduced by Rudin et al. [37]. The power of this method is that the intensity variations within phases are smoothed while the contrast at the boundaries between different phases is preserved. The smoothed images then were segmented into three phases (i.e. solid, water, and gas) by maximizing image entropy [38]. Local information of voxels was also included in the thresholds by coupling Kapur et al. [38] method with conditional dilation [37]. Segmented images were quantified using Minkowski functionals [39], and gas clusters were analyzed by modifying the labeling method of Hoshen and Kopelman [40]. All steps of image processing were performed using QuantIm-library (Available online at http://quantim.ufz.de). Segmented images were imported to VGStudioMax software for visualization.

2.2.5. Cluster statistics: cluster size distribution Fig. 5(A) shows the relative frequencies fr(s) of trapped gas clusters of size s (s – voxel number) for experiment 7. The correspondP ing cumulative cluster size distribution CDFðsÞ ¼ s1 frðsÞ is shown in Fig. 5(B). We found an approximately power-law size distribution for all experiments by fitting a continuous power func-

2.2.3. Flow rate dependence of trapped gas phase The most common way of correlating experimental data for trapped nonwetting phase is in terms of a relationship between the residual nonwetting phase saturation Srnw and the capillary number. Such a relationship is usually called capillary desaturation curve (CDC) [23]. In most cases the CDC is normalized by its corresponding plateau value, i.e. its limit for small capillary numbers (compare with Fig. 14.6 in [23]; with Fig. 13 in [14]; and with

to the experimental data (k – constant). tion CDFðsÞ ¼ 1  ðk=sÞ For experiment 7 one obtains an exponent s = 1.96 with a relative error of 12%. The least-square fit was conducted for the complete data set (4206 clusters, 774 cluster size groups, largest cluster size = 8442 voxel). Note that Fig. 5(A) and (B) show only a section of the whole data set. We also used a Gamma-distribution proposed by Stauffer and Aharoni [31] (Eq. (31)) and obtained a better fit with a relative error of 5%. Analyzing all data of the first series of experiments (1 mm-GBS) we found an averaged s-value of 1.9 with

s1

H. Geistlinger, S. Mohammadian / Advances in Water Resources 79 (2015) 35–50

39

Fig. 2. Images showing the residual trapped air (red colored) after water imbibition for experiment 7: (A) 2D-cut through a 3D-reconstructed image, (B) 3D-slice of a 3Dreconstructed image. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 3. Capillary desaturation curves for four different classes of glass beads.

an averaged relative error of 9%. Comparing our s-values with data in the literature, we found good agreement. For water-wet multiphase systems Iglauer et al. [8] found also a power-law distribution with an exponent s = 2.02–2.03 for residual oil clusters and s = 2.04 for residual gas clusters (for a review of s-values from literature, see Table 3 in [8]).

2.2.6. Surface area-volume relationship for trapped gas clusters We calculate for each trapped gas cluster i both its surface area Ai and its volume Vi. Fig. 6(A) shows a log–log plot and Fig. 6(B) a linear–linear plot of Ai(Vi) for experiment 7. We fit the data to a power-law relationship Ai / V pi and found for the log–log plot

p = 0.833 (relative error = 1=N

PN qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi2 ðAi;fit  Ai Þ =Ai = 7%; regression i

coefficient R2 = 0.98) and for the linear–linear plot p = 1 (relative error = 8%; R2 = 0.99). The difference in both fits is that for the log–log-plot small clusters have a larger weight compared to the linear–linear plot. Based on a log–log fit Iglauer et al. [8] obtained similar p-values of about 0.75. 2.2.7. Cluster shape and maximal z-extension The total number of trapped gas clusters within the column volume (=34 mL) for experiment 7 is 3695. First we consider gas bubbles with a radius smaller than 0.28 mm ( dk,max/2). Such gas bubbles are trapped within a single-pore. The number of

40

H. Geistlinger, S. Mohammadian / Advances in Water Resources 79 (2015) 35–50

Fig. 4. Total nw-phase surface area (gas: diamonds; oil: circles) versus nw-phase saturation. The data for the trapped oil phase are taken from Landry et al. [29].

Fig. 5. (A) Relative frequency distribution of cluster sizes s in voxel (blue circles) for experiment 7. The red curve represents the best fit to a power-law distribution with an exponent of 1.96 and the gray curve to a Gamma-distribution. (B) Normalized cumulative cluster sizes distribution CDF(s) (blue circles) and the best fits to a power-law- and Gamma distribution. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

single-pore-trapped gas bubbles is 1958 that is about 53% of the total number of trapped gas bubbles. However, they only contain 16% of the trapped gas volume. If we consider multi-pore-trapped

gas bubbles with a radius between dk,max/2 and the grain radius rk = dk/2 = 0.44 mm, then the number of such bubbles is 1367, i.e. 37% of all trapped gas bubbles. The corresponding gas volume

H. Geistlinger, S. Mohammadian / Advances in Water Resources 79 (2015) 35–50

41

Fig. 6. (A) Log–log-plot of the cluster surface area Ai versus cluster volume Vi. (B) Linear–linear plot.

amounts 40%. Although 90% of all trapped gas bubbles are smaller than a grain, they contain only 56% of the total gas volume. The remaining half of the gas volume is trapped by large network-like gas clusters (ganglia) shown in Fig. 7. Typically, these large clusters consist of a network of connected pores. The physical reason for such unexpected, ganglia-like clusters in strongly wet gas–water systems is the failure of thick continuous water films within the domain of these large clusters. Therefore, unstable interfaces cannot occur and snap-off into smaller clusters is suppressed [10]. Furthermore, we analyze the cluster height (=maximal z-extension) of all trapped large gas clusters shown in Fig. 8. It is obvious that no experimental cluster height exceeds a certain threshold value of about 4 mm. This demonstrates the influence of buoyancy introducing a natural cut-off length, which is one key prediction from percolation theory as will be discussed in Section 4.4. Furthermore, large cluster are of central importance for deriving the microscopic Ai–Vi-relationship, because they contain about half of the trapped gas volume.

3. Theoretical considerations of the fluid–fluid displacement process 3.1. Pore-throat network, geometry of the glass beads packing We can idealize our 1 mm-glass bead packing as a simple cubic pore-throat network of spherical pores with the radius rp

connected by cylindrical throats with a radius rt shown in Fig. 9. Fig. 9(A) shows the 3D-pore-throat network (blue-shadowed spheres and red cylinders) and grains in between (yellow spheres) that are shown separately in Fig. 9(B). In a more realistic packing the coordination number Z of adjacent grains or pores will change between simple cubic (=sc, i.e. loosest packing with Z = 6) and facecentered-cubic (=fcc) packing (densest packing with Z = 12). We can estimate the experimental coordination number Zexp using a /(Z)-relationship [23] for close random packing of equal-sized spheres:

/ ¼ 0:0043  Z 2  0:1193  Z þ 1:072

ð1Þ

(/  porosity). For the experimental porosity of 0.332 one obtains Zexp = 9.36. Important packing parameters are the minimal and maximal pore diameter dk,min and dk,max. We estimate the throat and pore radius by the minimal and maximal diameter:

rt ¼ dk;min =2 ¼ nmin ðZ exp Þ  dk =2

ð2aÞ

rp ¼ dk;max =2 ¼ nmax ðZ exp Þ  dk =2

ð2bÞ

where dk denotes the grain diameter, and nmin, nmax interpolating function between sc- and fcc-packing with nmin = 0.27 and nmax = 0.58. This gives a throat radius rt = 0.11 mm and a pore radius rp = 0.25 mm. If one uses the experimental PSD and estimates the throat radius by rp,mean  r = 0.11 mm and the pore radius by rp,mean + r = 0.23 mm (see Table 2), one realizes that the interpolation algorithm yields reasonable values (for details see [35]).

42

H. Geistlinger, S. Mohammadian / Advances in Water Resources 79 (2015) 35–50

Fig. 7. Different types of network-like trapped gas clusters. The maximal z-extension reaches 50–60 voxels, i.e. 4–5 mm.

Fig. 8. Cluster height (maximal z-extension) of all trapped large gas clusters for experiment 7.

Obviously, cylindrical pore throats shown in Fig. 9 are a poor approximation of the star-shaped cross section shown in Fig. 10(A). Therefore we approximate the star-shaped cross section by a rectangular quadratic pore throat with a throat width x = 2rt = 0.22 mm and a throat length lt = dk/2 = 0.425 mm shown in Fig. 10(B). Following Lenormand and Zarcone [11] we approximate the spatial fluid–fluid distribution by a core-annular flow shown in Fig. 10(C), where the core fluid is the nonwetting fluid – air in our case – and the annular fluid is the wetting fluid – water in our case. Furthermore, we assume that the whole packing of glass spheres is covered by a thin water film of approximately 300 nm [24,10]. This thin water film is shown in Fig. 10 as a blue line around all grains. 3.1.1. Time scales for thin film flow, bulk advance and throat filling In this Section we will show that thin film flow and throat flow can occur at times scales that are comparable with the duration of the experiments. For the lowest flow rate (WT-rising velocity = 0.1 cm/min, Ca = 2  107) a precursor thin film is therefore

likely to cover the whole sediment within the column. As a second consistency test we will study the cooperative I1-pore filling mechanism, which is the basic trapping mechanism in invasion percolation and demonstrate that I1-displacement can lead to an unstable interface configuration of the thin film that can cause film swelling and snap-off. In this case the time scale for throat filling is about 1.5 min and therefore much smaller than the duration of the imbibition experiment (see Appendix A). As pointed out above the interaction between precursor thin film flow, corner flow and bulk advance lead to different cooperative pore and throat filling mechanisms. These are determined by the corresponding capillary pressures and unstable configurations of the gas–water interface. Based on a realistic geometry, topology and a square throat cross section we will discuss the relevant flow types that can occur in close random glass beads packings. The relevant local capillary pressures that drive the flow of the wetting fluid (complete wetting), cooperative pore filling and interface instability in a square throat of width x are given by the following expressions [11]:

43

H. Geistlinger, S. Mohammadian / Advances in Water Resources 79 (2015) 35–50

PI1 ¼

2r R1

ð3cÞ

and for unstable snap-off configuration within the throat induced by swelling water films

Ps ¼

Fig. 9. Pore-throat-network for simple cubic packing of 1 mm-glass beads Fig. (A): red – throats that are connected to pores, yellow – grains. Fig. (B) shows only the grains, i.e. the equal-sized 1 mm-glass beads. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

2r r ¼ : x rt

ð3dÞ

Choosing as reference pressure the pressure of the nonwetting phase (pnw = 0), the pressures at the front of the wetting phase are given by the negative values of the corresponding capillary pressures. As Lenormand and Zarcone [11] has emphasized there is always a sequential order of the relevant displacement processes. Since pistonlike displacement needs a frontal meniscus over the whole throat, the flow rate has to be sufficiently high (Ca > 104) to provide enough wetting fluid for the filling process. As the capillary number decreases, thin film flow and corner flow ‘‘have enough time’’ to occur. In this case pistonlike displacement is only possible after snap-off or during the I1-displacement process. Comparing the order of the absolute values of the pressure gradient thin film flow has the highest value followed by corner flow, and I1-displacement, i.e. the smaller the invaded cross section (equivalent radius requ.) the higher the driving pressure gradient. On the other hand the flow rate is proportional to the fourth power of requ. and therefore the flow resistance in smaller ducts is significantly higher. In order to estimate times scales for typical displacement types, one has to take into account both effects. To illustrate the hierarchical order of the flow processes, we calculated the time behavior of the precursor thin film flow followed by pistonlike throat flow based on the analytical solution proposed by Lenormand and Zarcone [11], modified by Blunt and Scher [14], and Constantinides and Payatakes [17] (for details see Appendix A). The model assumes that thin-film flow is caused by surface roughness, i.e. by small crevices or surface channels along the grains, where these channels are tiny cylindrical tubes of radius rc, spaced a distance 4rc apart. The effective radius rc of such channels is determined by the thin film thickness d. Using d = 300 nm as a realistic estimate for glass beads [24], one obtains rc = 380 nm. Applying a constant flow rate q0 the onset of pistonlike flow in a throat needs a threshold time t0 to establish a pressure Pp (Eq. (3a)) at the entrance of the throat. During this time the thin film advances a threshold length l0 driven by the viscous pressure difference DPw = Pthin film  Pp. The onset of throat flow is shown in the upper Fig. of Fig. 11(A). Calculating the threshold values for our 1 mm-glass beads for a capillary number of 106, i.e. a mean WT-rising velocity of 0.5 cm/min, one obtains t0 = 2 s and l0 = 3 mm. That is the thin film covers only 3 layers (3 grain diameter) ahead of the throat meniscus. For times larger than t0 the constant flow rate has to ‘‘feed’’ both throat and thin film flow (see Fig. 11(A), lower figure). The time development of the precursor thin film region l(t) in dimensionless variables (dimensionless time s = t/tscale, k = l/lscale) is described by its inverse function:



sðkÞ ¼ s0  ðk  k0 ÞLn

k1 k0  1

 ð4Þ

For pistonlike displacement within a throat

Pp ¼

4 2r ¼ x rt

ð3aÞ

for corner and thin-film flow (flow in crevices [14]; rc – radius of the curvature of the interface within the crevice, r – interface tension)

Pthin film ¼

2r rc

ð3bÞ

for cooperative I1-displacement (R1 – curvature of the meniscus within the pore)

(see Appendix A). We evaluate l(t) for two times: tmin = 10 min and tmax = 50 min, i.e. the corresponding times for water table rising at v = 0.5 cm/min and 0.1 cm/min, respectively. During tmin the precursor film region has extended from 3 to 111 mm, i.e. covers about one third of the column length. For smaller WT-rising velocity, i.e. smaller injection rate, one obtains for the threshold values: t0 = 50 s and l0 = 15 mm. During tmax the throat meniscus has moved 50 mm and the precursor film has extended from 15 to 240 mm. Hence, the thin film front has nearly reached the end of the column (length = 300 mm) and wets the whole sediment.

44

H. Geistlinger, S. Mohammadian / Advances in Water Resources 79 (2015) 35–50

Fig. 10. Pore throat geometry for simple cubic packing of glass spheres after Lenormand and Zarcone [11]. (A) is a cross section along the dashed line of (B) and (C) is an rectangular approximation of the star-shaped throat shown in (A). The wetting fluid (blue) flows along the corners, whereas the nonwetting fluid (white) flows through the throat centre/bulk (core-annular flow). Note that the grains are covered by a thin wetting film (blue solid line). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 11. Precursor thin film flow (light gray) and throat flow (dark gray) according to [11]: (A) time dependence of the precursor thin film region l(t): upper figure: threshold (t0, l0) at which throat flow starts, lower figure: cooperative thin film and throat flow. (B) Flow of a wetting fluid through a porous media with capillary trapping of the nonwetting phase (white) and some throat filling within the precursor film region by snap-off instability. (A) and (B) are modified figures taken from [11,14].

An interesting question is, how much time is needed for throat filling in front of the mean bulk front. Choosing a throat position of 5 cm (=height of the WT-rise) ahead of the bulk front, one obtains a typical time scale of about 1.5 min (see Appendix A). Since this time is smaller than the final time of the WT-rising, tmin, it is likely that throat filling will occur; schematically shown in Fig. 11(B). 3.1.2. Occurrence of unstable thin water films during cooperative I1pore filling It is instructive to study, how instability and snap-off can occur during cooperative I1-pore filling. Since the capillary pressure for pistonlike advance in a quadratic throat (Eq. (3a)) is always larger than the snap-off capillary pressure (Eq. (3d)), instability at the fluid–fluid-interface and hence snap-off can never arise. However, if the front invades the pore, the curvature of the front meniscus decreases as described by Eq. (3c). If the capillary pressure at the front becomes smaller than a critical value, snap-off can occur in the adjacent throat (throat 1 in Fig. 12). This is shown schematically in Fig. 12 for the I1-mechanisms. As discussed above and shown in Fig. 12 the grains are covered by thin wetting films. These films can swell because of an unstable fluid–fluid interface indicated by the dashed line. The menisci in throats 2 to 4 are shown at time t1 (thick solid line, higher capillary pressure) and at time t2 (dashed line, lower capillary pressure). Since the pressure of the nonwetting phase (gas phase) is

connected to the outlet and is therefore constant (viscous pressure drop is neglected), the wetting phase pressure at the front, pw = pnw  pc, has to increase as the menisci move pistonlike in throats 2–4 from t1- to t2-position. The critical value for snap-off is given by the capillary pressure within the throat (Eq. (3d)). Using the derived values for throat radius rt = 0.11 mm and pore radius rp = 0.25 mm, the critical capillary pressure is about 639 Pa. Since this pressure is larger than PI1 = 594 Pa (R1 = rp) instability can occur, i.e. the thin wetting film swells in the throat as long as the throat is filled.

4. Comparison of experimental results with predictions from percolation theory If capillary forces are the dominant forces at pore scale and viscous forces can be neglected, the flow patterns can be described by percolation theory. Let us estimate the viscous pressure drop across the throat length lt using the generalized Darcy equation for the wetting fluid

DPv isc 

lw q w lt k

ð5aÞ

where qw is the flow rate of the wetting fluid, lw its viscosity, and k its absolute permeability, and the relative permeability is of order 1.

H. Geistlinger, S. Mohammadian / Advances in Water Resources 79 (2015) 35–50

45

Fig. 12. Fluid–fluid-displacement in pore throats for the cooperative I1-pore filling mechanism. Pistonlike displacement occurs in throats 2, 3, and 4. Snap-off displacement by swelling thin wetting films occurs in throat 1.

Comparing this to a typical capillary pressure within a throat (Eq. (3a)), the ratio of viscous to capillary forces is

DPv isc r t lt n  Ca  / ¼ Ca  / min ¼ 17  Ca Pc 2k 2d

ð5bÞ

where we have used an empirical relationship between permeability and grain size for random close packing of spheres of radius rk, k ¼ d  r2k , d  2.7  103 [42,43]. For the range of experimental capillary numbers 2  107 . . . 106 capillary forces are by 5–6 orders of magnitude larger than viscous forces. Consequently, we expect that percolation theory will describe our experimental results. 4.1. Flow rate dependence of trapped gas phase Blunt and Scher [14] developed a 3D network model that describes thin film flow and bulk advance as well as cooperative pore filling mechanisms and discussed the different flow types that can occur at different flow rates. Of special interest are the investigations of the change in trapped nonwetting phase saturation with flow rate. They found that the residual nonwetting phase saturation is nearly constant, if aCa is in the range between 104 . . . 103 (see Fig. 13 in [14]), where a = b(lt/rt)3 and b = 109 is a geometrical factor of the Hagen-Poiseuille flow in a square throat for a wetting fluid with a contact angle of zero and no slip at the nonwetting/wetting interface. The gradual change in residual saturation at very low flow rates is caused by the transition from a flow regime dominated by snap-off to a connected, invasion percolationlike advance. Both flow types are realistic in the considered range of experimental capillary numbers 107 . . . 106 (see Table 1). For this range we calculate the parameter aCa = 1.4  104. . .7.0  104. Consequently, one expects no flow rate dependence in the considered range of experimental water

table rising velocities. This percolationlike behavior of the CDCs is shown in Fig. 3. We note that the percolation limit (solid line) is a mean value over an ensemble of realizations [14], i.e. the individual realizations lead to some fluctuations around this mean value. This is demonstrated for the first experimental series. As shown in Fig. 1 all the PSDs are nearly identical. Hence this experimental series represents a real Monte-Carlo experiment drawing different stochastic realizations from the same probability density function shown in Fig. 1. 4.2. Linear relationship between nonwetting phase surface area and nonwetting phase saturation To understand the linear relationship between the surface area and the trapped gas volume of the nonwetting phase (see Fig. 4), we first study the relationship between the surface area of a single trapped gas cluster and its volume. Consider a cluster of size s that forms a spherical bubble, then one would expect that its surface area is proportional to the 1/D-root of its volume

Ai  s11=D

ð6Þ

i.e. for dimension D = 3 an exponent p = 2/3. Obviously, this exponent does not agree with our experimental findings of the first experimental series, namely where p has values between 0.833 and 1 (see Fig. 6). Following Stauffer and Aharony (see Ch. 3.1. in [31]) one has to imagine that large clusters are always penetrated by holes like a ‘‘swiss cheese’’ and therefore their total surface area is a sum of an outer and inner surface area. Thus, one expects that for large, but finite clusters

Ai  s 1 :

ð7Þ

In network-like clusters as shown in Fig. 7 one can see these holes and obviously the number of connected pores Np gives both the

46

H. Geistlinger, S. Mohammadian / Advances in Water Resources 79 (2015) 35–50

cluster surface area, Npa, and the cluster volume, Npv, where a is the surface area and v the volume of some averaged pore. We note that Stauffer and Aharony [31] give a strict proof that for sufficiently large clusters the surface area is always proportional to its volume (see Eq. (44b) in [31]). Since the large clusters contain nearly 50% of the trapped gas volume the linear–linear plot shown in Fig. 6(B) is the more realistic plot, because the large clusters have a higher weight compared to the log–log plot shown in Fig. 6(A). Therefore, the linear relationship of the total gas surface area on gas saturation shown in Fig. 4 is a direct consequence of the linear dependence of large network-like clusters. Consequently, our experimental results are in excellent agreement with percolation theory. We analyzed the data from Al-Raoush and Willson (Fig. 10 in [32]) and found a strong linear Ai–Vi-relationship with a regression coefficient of 0.96. This result indicates that percolation theory holds also for an oil-wet multiphase system (oil-wet glass beads, oil – wetting phase, brine – nonwetting phase). 4.3. Cluster statistics Following Stauffer and Aharony (see Ch. 2.2. in [31]) then the probability that the cluster contains exactly s sites (nodes, pores) P is PDFðsÞ ¼ ns  s= ns  s, where ns defines the normalized cluster number, i.e. the cluster number per lattice site. The cumulative cluster size distribution is then given by s X CDFðsÞ ¼ PDFðsÞ:

ð8Þ

1

At the percolation threshold pc (or the terminal point, where the non-wetting fluid becomes disconnected; see [14]) and for large clusters the normalized cluster number shows universal power law decay with the Fisher exponent s [44]:

ns ðpc Þ / ss

ð9Þ

which gives

CDFðs; pc Þ / s1s :

ð10Þ

The critical exponents for 3D are listed in Table 2 in [31]: s = 2.18 and m = 0.88 (see Eq. (12)). These critical exponents are the same both for ordinary percolation and invasion percolation [45,46]. Hence, analyzing the critical exponent s is not a unique criterion to distinguish between both percolation processes, i.e. between bond percolation due to precursor thin film flow and throat filling by snap-offs and invasion percolation due to pistonlike bulk advance and cooperative I1-pore filling. Note that s is the number of sites (bonds) that corresponds to pores (throats) and not to the number of voxels. In order to compare the experimental distribution shown in Fig. 5 with the universal power law behavior Eq. (10), one has to renormalize the experimental distribution by the number of voxels that one mean pore contains. The renormalization yields CDFð~sÞ, where ~s now describes the number of pores (=sites) spanning a cluster. For regular packings of glass beads the use of a mean pore size should be reasonable and was taken from the corresponding PSD of each experiment (about 35 voxels, see Table 1). The best fit to the data for experiment 7 is shown in Fig. 13(A). Note that the cluster size s is now given in sites. The renormalized exponent s = 2.16 for experiment 7 and the averaged s-value of 2.15 show now excellent agreement with universal power law behavior, i.e. the deviation is smaller than 2%. To our knowledge this renormalization is not taken into account in the literature. For example Iglauer et al. [2] compare their voxelbased cluster-size distribution to percolation theory and contend a good agreement to the prediction from percolation theory: . . .’’In a

strongly water-wet system (oil-brine), we find an approximately power-law size distribution with s = 2.05, close to the prediction from percolation theory.’’. . . We note that this renormalization introduces a lower cutoff for clusters that contain less than 35 voxels (see similar discussion in [28]). Since the universal power law holds only for large clusters this cutoff should have no significant impact on s. Instructive is the comparison with the oil-wet multiphase system. Since Al-Raoush and Willson [32] derived a mean pore size of 58 lm, we were able to renormalize their cumulative cluster size distribution shown in their Fig. 6. The data show a universal power-law behavior with s = 1.91 (see Fig. 13(B)). The deviation from the exact curve (dashed curve in Fig. 13(B)) may be attributed to the sparse data set (only 75 clusters).

4.4. Maximal z-extension of trapped gas clusters Now we study the effect of buoyancy on the size of trapped gas clusters following Wilkinson [45] and Blunt et al. [46]. The relative importance of capillary to buoyancy forces can be estimated by the bond number

B¼

Dqgr t lt : 2r

ð11Þ

B represents the ratio of the buoyancy (hydrostatic) pressure in a single throat of height lt to the throat capillary pressure (Eq. (3a)). Consider the imbibition of the wetting fluid water from below, the occupation probability is reduced by the hydrostatic gradient Dqgz (Dq – density difference of fluids, g – earth acceleration) and hence, it varies linearly with the cluster height z. Therefore, the buoyancy force will affect the most extended clusters and truncate ns(pc). Physically, this leads to a maximum correlation length nB in the system according to

nB ¼ lt

 a 1 B

ð12Þ

with the critical exponent a = m/(1 + m) = 0.468. The interpretation of the length nB in the imbibition case is that it is an upper limit of the largest trapped clusters of the nonwetting fluid and amounts to 6 mm for trapped gas clusters within the 1 mm-glass beads. If one compares this threshold value with the experimental cluster height (maximal z-extension) of the largest trapped gas clusters shown in Fig. 8, it is obvious that no experimental cluster height exceeds this value, even though some of the clusters come near to this value, i.e. have a height of about 4 mm. We also calculate the cutoff-length for the oil-wet multiphase system [32] estimating the throat length to 100 lm by visual inspection of Fig. 2 in [32] and obtained a cutoff-length of about 5 mm. Comparing this value with the maximal extension of the largest clusters (about 3 mm; see Fig. 8 in [32]) suggests that no cluster exceeds this natural cutoff-length. Consequently, the prediction from percolation theory holds for the oil-wet multiphase system, too. We note that in l-CT studies where the natural cutoff-length has the same order as the sample dimensions this could cause a cut-off for large trapped clusters and could lead to incomplete data sets. As an instructive example we calculated the cutoff-length for the l-CT study conducted by Georgiadis et al. [30] and obtained about 6 mm (mean pore radius 35 lm, interface tension 50.5 mN/m; B ¼ Dqgr 2p =r [45]). We believe that this is the reason for the incomplete data set measured in [30] and the high sensitivity of the power-law fit to the upper cluster size limit as discussed in [30].

H. Geistlinger, S. Mohammadian / Advances in Water Resources 79 (2015) 35–50

47

Fig. 13. Best fit of experimental data (blue circles) to a power-law distribution (red curves): (A) fit to the data of experiment 7. (B) Fit to the experimental data from Al-Raoush and Willson (Fig. 6 in [32]). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

5. Conclusions 5.1. Consistency test of proposed capillary trapping mechanisms The fluid–fluid displacement and the capillary trapping mechanism in the range of small capillary numbers 2  107. . . 106 can be described (i) by ordinary (bond) percolation or (ii) by invasion percolation. The interplay between thin film flow, pistonlike bulk advance and snap-off due to unstable interface configurations are essential to understand the experimental results on capillary trapping in 1 mm-glass beads packings. Based on experimental pore size distributions and packing parameters we derived geometrical parameters of an idealized pore-throat network, where a realistic square throat was used. We showed that thin film flow and throat flow occur in times scales that are comparable with the duration of the experiments. For the lowest flow rate (WT-rising velocity = 0.1 cm/min, Ca = 2  107) a precursor thin film is therefore likely to cover the whole sediment within the column. This can lead to some throat filling in front of the mean bulk front. Consequently, bond percolation could be responsible for the trapping mechanism. As a second test we study the cooperative I1-pore filling mechanism, which is the basic trapping mechanism in invasion percolation. We found that I1-displacement leads to an unstable interface configuration of the thin film that will cause film swelling

and snap-off. The time scale for throat filling is about 1.5 min and therefore much smaller than the duration of the imbibition experiment. These findings are consistent with the capillary trapping mechanism of invasion percolation. 5.2. Comparison with percolation theory, new contributions We compare our experimental results of a strongly-water wet multiphase system with predictions from percolation theory and found excellent agreement. Percolation theory explains our central experimental result, namely that the capillary desaturation curves are not dependent on flow rate in the range of capillary numbers 2  107–106. Compared to recent l-CT studies [8,26,28–30,32] our paper makes the following new contributions: (a) We renormalize the measured cluster size distribution and found excellent agreement with the predicted s-value from percolation theory. (b) We realized that the macroscopic linear Anw–Vnw-relationship is a direct consequence of the microscopic linear Ai– Vi-relationship. (c) Realizing that about 50% of the residual gas volume is trapped by large gas clusters (about 10% of all clusters) the log–log-plot gives the large clusters a too low statistical weight and small gas clusters a too high statistical weight.

48

H. Geistlinger, S. Mohammadian / Advances in Water Resources 79 (2015) 35–50

Consequently, we used instead of the log–log plot a linear– linear plot and obtained a strong linearity (regression coefficient = 0.99) that is predicted from percolation theory. (d) We study the influence of gravity (buoyancy) and analyzed the maximal length of large clusters in z-direction and found reasonable agreement with the cut-off length nB predicted by percolation theory. (e) Analyzing the data set for an oil-wet multiphase system (oilwet glass beads, oil – wetting phase, brine – nonwetting phase) [32] we found reasonable agreement with predictions from percolation theory.

estimate the precursor thin film flow region l(t) (see Fig. 11), we assume that near the inlet the nonwetting fluid is displaced by pistonlike movement from throats and by the I1-displacement from pores. We neglect the viscous pressure drop in the nonwetting fluid and choose the pressure of the nonwetting phase as reference pressure, i.e. the pressure in the wetting phase is negative. We assume that the thin film flow is caused by surface roughness and that it will flow along surface crevices (noted as surface channels with a radius of rc and a spacing of 4rc). Then the number of such channels in a square throat (width x) is

x : rc

Nc ¼ Since the key predictions from percolation theory (linear Ai–Virelationship, universal power-law, and existence of a cut-off length) hold for large clusters, which can span over 4–5 mm as shown in Fig. 7, the relative error of the limited resolution of l-CT (85 lm) is about 2%. We think that this is the reason for the excellent agreement between our experimental results and the predictions from percolation theory. We note that there is an inherent connection between universal power-like behavior of the cluster-size distribution and the residual gas saturation determined through the cut-off length. The summation of the cluster-size distribution yields the volume of the P residual non-wetting phase: Sðsmax Þ ¼ s1max PDFðsÞ  s, i.e. the larger the upper summation limit smax (=cut-off length nB) the higher the trapped volume. Since the cut-off length depends on the density difference, interfacial tension and pore-throat geometry, the different cut-off lengths will naturally lead to different residual nonwetting phase saturations. As discussed in Geistlinger et al. [10] for oil and CO2 there are larger clusters indicating a larger cut-off length and consequently a higher residual non-wetting phase saturation, which is indeed observed [2,26,28]. Moreover, the universal power-like behavior is independent of details of the local geometry of the pore space, like coordination number and poresize distribution. Thus it will be also valid for other types of porous media, like natural sands and sandstones. Further research is necessary to establish that all these multiphase systems belong to the same universality class [45]. 5.3. Is capillary trapping caused by invasion percolation or by bond percolation? To decide if capillary trapping is cause by invasion percolation, i.e. by unstable I1-displacement or by bond percolation, i.e. by precursor thin film flow, snap-off and throat filling, we conduct preliminary experiments in 1 mm-glass-beads monolayer. There is some experimental evidence that the trapping mechanism is caused by invasion percolation and not by ordinary bond percolation. However, to quantify the differences of capillary trapping between horizontal and vertical flow, for instance the impact of gravitational instability [47], further work is necessary. Acknowledgments The authors gratefully acknowledge funding of the project Dynamically Capillary Fringe (DYCAP) by the German Research Foundation DFG and Dr. Steffen Schlüter and Prof. Hans-Jörg Vogel for critical discussion and valuable hints. Furthermore, we thank MSc. Iman Ataei for conducting the 1-mm-monolayer experiments.

ðA1Þ

The effective channel radius is defined through the film thickness d:

4d

rc ¼

p

ðA2Þ

:

The wetting fluid is injected with constant flow rate q0 in a throat of infinite length. As long as the absolute value of the wetting phase pressure is smaller than the capillary pressure of the throat no meniscus can be established and the fluid will only flow as thin film along the surfaces up to a threshold time t0, where throat flow starts. We calculate the threshold length of the thin film region l0 assuming Hagen–Poiseuille flow in each surface channel

q¼

p 8l

r 4c

jDpj ; l0

ðA3Þ

where the absolute value of the pressure difference in the wetting fluid is given by

jDpj ¼

2r 4r  :: rc x

ðA4Þ

For t < t0 the flow rate q0 is equal to the thin film flow rate, i.e.

q0 ¼ qthin film ¼ Nc q:

ðA5Þ

Inserting (A1), (A3) and (A4) into the threshold condition (A5), one obtains the threshold length for rc  x

l0 ¼

p r2c 4 x  Ca

;

ðA6Þ

where the throat capillary number is defined as

Ca ¼

q0 l : x2 r

ðA7Þ

Filling up Nc surface channels of length l0 yields the threshold time

t0 ¼

p2

r3c : 4 Ca  q0

ðA8Þ

When the threshold pressure has been reached, a meniscus is formed within the throat and pistonlike flow through the throat starts. For t > t0 the flow rate q0 will ‘‘feed’’ both throat and thin film flow:

q0 ¼ qthroat þ qthin film ¼ qthroat þ Nc  q:

ðA9Þ

By definition the throat flow rate at position z(t) (see Fig. 11) is given by

qthroat ¼ S

dz ; dt

ðA10Þ

Appendix A

and the channel flow rate at position z(t) + l(t) by

A.1. Extent of precursor thin film flow and time scale for throat filling

q¼s

The following considerations are partly adapted from Lenormand and Zarkone [11] and Blunt et al. [14]. In order to

where S = x2 and s = pr 2c define the cross section of the throat and the channel, respectively. Inserting (A10) and (A11) into (A9) and

dðzðtÞ þ lðtÞÞ ; dt

ðA11Þ

H. Geistlinger, S. Mohammadian / Advances in Water Resources 79 (2015) 35–50

using dimensionless variables yields the differential equation for the time development of the precursor thin film region:

dk kþk ¼ 1 þ g; ds

ðA12Þ

where g = s/S is the ratio of the cross sections, k the dimensionless length and s the dimensionless time that are defined as follows

k¼

l lscale

;

s¼

[7]

[8]

[9]

t t scale

ð13abÞ

;

[10]

with the scaling factors

lscale ¼

rc ; 4  Ca

t scale ¼

r c x2 : 4Ca  q0

[11]

ð14abÞ

The analytical solution k(s) for g  1 is given as its inverse function

[12]

s(k)

sðkÞ ¼ s0  ðk  k0 ÞLn



 k1 : k0  1

[13]

ðA15Þ

Using the scaling factors (14a,b) the dimensionless threshold length and threshold time are

k0 ¼

pr c x

;

s0 ¼

p  r 2 c

2x

:

ðA16Þ

The steady-state is nearly reached at s = 3. At this time the higher capillary pressure of the surface channel is compensated by an increasing length-dependent friction force. The steady-state extent of the thin film region is

k1 ¼ 1:

[14] [15]

[16]

[17]

[18]

ðA17Þ [19]

Time scale for throat filling by snap-off within the precursor thin film region The time that is needed to fill a throat within the precursor thin film region is given by

t¼

8lxlt l ; pr3c jDpj c

ðA18Þ

where lt is the throat length, rc the channel radius, x the throat width, jDpj ¼ 2r(1/rc  1/x) the absolute value of the pressure difference, and lc the throat position within the thin film region. Appendix B. Supplementary data

[20]

[21]

[22]

[23]

[24]

[25]

Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.advwatres.2015. 02.010.

[26]

References

[27]

[1] Juanes R, Spiteri EJ, Orr Jr FM, Blunt MJ. Impact of relative permeability hysteresis on geological CO2-storage. Water Resour Res 2006;42:W12418. http://dx.doi.org/10.1029/2005WR004806. [2] Iglauer S, Paluszny A, Pentland CH, Blunt MJ. Residual CO2 imaged with X-ray microtomography. Geophys Res Lett 2011;38(L21403):2011G. http:// dx.doi.org/10.1029/2011GL049680. [3] Johnson PC, Johnson RL, Brucea CL, Leeson A. Advances in in situ air sparging/biosparging. Biorem J 2001;5:251–66. http://dx.doi.org/10.1080/ 20018891079311. [4] Jeong SW, Corapcioglu MY. Force analysis and visualization of NAPL removal during surfactant-related floods in a porous medium. J Hazard Mater 2005;126:8–13. http://dx.doi.org/10.1016/j.jhazmat.2005.06.015. [5] Kaoa CM, Chena CY, Chenb SC, Chiena HY, Chen YL. Application of in situ biosparging to remediate a petroleum-hydrocarbon spill site: field and microbial evaluation. Chemosphere 2008;70:1492–9. http://dx.doi.org/ 10.1016/j.chemosphere.2007.08.029. [6] Geistlinger H, Lazik D, Krauss G, Vogel H-J. Pore-scale and continuum modeling of gas flow pattern obtained by high-resolution optical bench-scale

[28]

[29]

[30]

[31] [32]

49

experiments. Water Resour Res 2009;45:W04423. http://dx.doi.org/10.1029/ 2007WR006548. Chatzis I, Morrow NR, Lim HT, et al. Magnitude and detailed structure of residual oil saturation. In: Paper SPE/DOE-10681, presented at the 3rd symposium on enhanced oil recovery. Tulsa; April 4–7, 1982. http://dx.doi. org/10.2118/10681-PA. Iglauer S, Paluszny A, Blunt MJ. Simultaneous oil recovery and residual gas storage: a pore-level analysis using in situ X-ray micro-tomography. Fuel 2013;103:905–14. http://dx.doi.org/10.1016/j.fuel.2012.06.094. Mohammadian S, Geistlinger H, Vogel H-J. Quantification of gas phase trapping within the capillary fringe using micro-CT, Special section: Dynamic Processes in Capillary Fringes, Vadose Zone J. http://dx.doi.org/10.2136/vzj2014.06.0063. Geistlinger H, Mohammadian S, Schlueter St, Vogel H-J. Quantification of capillary trapping of gas clusters using X-ray microtomography. Water Resour Res 2014;50:4514–29. http://dx.doi.org/10.1002/2013WR014657. Lenormand R, Zarcone C. Role of roughness and edges during imbibition in square capillaries. In: SPE-paper No. 13264, in Proceedings of the 59th ann. tech. conf. of the SPE. Houston, TX (SPE, Richardson, TX; 1984. http://dx.doi. org/10.2118/13264-MS. Lenormand R, Touboul E, Zarcone C. Numerical models and experiments on immiscible displacements in porous media. J Fluid Mech 1988;189:165–87. http://dx.doi.org/10.1017/S0022112088000953. Vizika O, Avraam DG, Payatakes AC. On the role of the viscosity ratio during low-capillary number forced imbibition in porous media. J Colloid Interface Sci 1994;165:386–401. http://dx.doi.org/10.1006/jcis.1994.1243. Blunt MJ, Scher H. Pore-level modeling of wetting. Phys Rev E 1995;52:6387–403. http://dx.doi.org/10.1103/PhysRevE.52.6387. Hashemi M, Dabir B, Sahimi M. Dynamics of two-phase flow in porous media: simultaneous invasion of two fluids. AIChE J 1999;45:1365–82. http:// dx.doi.org/10.1002/aic.690450702. Hashemi M, Sahimi M, Dabir B. Monte Carlo simulation of two-phase flow in porous media: invasion with two invaders and two defenders. Phys A 1999;267:1–33. http://dx.doi.org/10.1016/S0378-4371(98)00661-X. Constantinides GN, Payatakes AC. Effects of precursor wetting films in immiscible displacement through porous media. Transp Porous Media 2000;38:291–317. http://dx.doi.org/10.1023/A:1006557114996. Joekar-Niasar V, Hassanizadeh SM, Pyrak-Nolte LJ, Berentsen C. Simulating drainage and imbibition experiments in a high-porosity micromodel using an unstructured pore network model. Water Resour Res 2009;45(W02430):2007W. http://dx.doi.org/10.1029/2007WR006641. Joekar-Niasar V, Doster F, Armstrong RT, Wildenschild D, Celia MA. Trapping and hysteresis in two-phase flow in porous media: a pore-network study. Water Resour Res 2013;49:4244–56. http://dx.doi.org/10.1002/wrcr.20313. Hilpert M, Miller CT. Pore-morphology-based simulation of drainage in totally wetting porous media. Adv Water Resour 2001;24:157–77. http://dx.doi.org/ 10.1016/S0309-1708(00)00056-7. Payatakes AC, Dias MM. Immiscible microdisplacement and ganglion dynamics in porous media. Rev Chem Eng 1984;2:85. http://dx.doi.org/ 10.1515/REVCE.1984.2.2.85. Celia MA, Reeves PC, Ferrand LA. Recent advances in pore scale models for multiphase flow in porous media. Rev Geophys 1995;33:1049–57. http:// dx.doi.org/10.1029/95RG00248. Sahimi M. Flow and transport in porous media and fractured rock. Weinheim: VCH-Wiley; 2011. http://dx.doi.org/10.1002/ 9783527636693. Kibbey TCG. The configuration of water on rough natural surfaces: implications for understanding air–water interfacial area, film thickness, and imaging resolution. Water Resour Res 2013;49:4765–74. http://dx.doi.org/ 10.1002/wrcr.20383. Schlüter S, Sheppard A, Brown K, Wildenschild D. Image processing of multiphase images obtained via X-ray microtomography: a review. Water Resour Res 2014;50:3615–39. http://dx.doi.org/10.1002/2014WR015256. Iglauer S, Fernø MA, Shearing P, Blunt MJ. Comparison of residual oil cluster size distribution, morphology and saturation in oil-wet and water-wet sandstone. J Colloid Interface Sci 2012;375:187–92. http://dx.doi.org/ 10.1016/j.jcis.2012.02.025. Lorenz CD, Ziff RM. Precise determination of the bond percolation thresholds and finite-size scaling corrections for the sc, fcc and bcc lattices. Phys Rev E 1998;57:230–6. http://dx.doi.org/10.1103/PhysRevE.57.230. Iglauer S, Favretto S, Spinelli G, Schena G, Blunt MJ. X-ray tomography measurements of power-law cluster size distributions for the nonwetting phase in sandstones. Phys Rev E 2010;82:056315. http://dx.doi.org/10.1103/ PhysRevE.82.056315. Landry CJ, Karpyn ZT, Piri M. Pore-scale analysis of trapped immiscible fluid structures and fluid interfacial areas in oil-wet and water-wet bead packs. Geofluids 2011;11:209–27. http://dx.doi.org/10.1111/j.14688123.2011.00333.x. Georgiadis A, Berg S, Makurat A, Maitland G, Ott H. Pore-scale microcomputed-tomography imaging: nonwetting-phase cluster-size distribution during drainage and imbibition. Phys Rev E 2013;88:033002. http://dx.doi.org/ 10.1103/PhysRevE.88.033002. Stauffer D, Aharony A. Introduction to percolation theory. revised second ed. Philadelphia: Taylor and Francis; 1994. Al-Raoush RI, Willson CS. A pore-scale investigation of a multiphase porous media system. J Contam Hydrology 2005;77:67–89. http://dx.doi.org/10.1016/ j.jconhyd.2004.12.001.

50

H. Geistlinger, S. Mohammadian / Advances in Water Resources 79 (2015) 35–50

[33] Sahimi M. Flow phenomena in rocks: from continuum models to fractals, percolation, cellular automata, and simulated annealing. Rev Mod Phys 1993;65:1393–534. http://dx.doi.org/10.1103/RevModPhys.65.1393. [34] Iglauer S, Salamah A, Sarmadivaleh M, Liu K, Phan Chi. Contamination of silica surfaces: impact on water–CO2-quartz and glass contact angle measurements. Int J Greenhouse Gas Control 2014;22:325–8. http://dx.doi.org/10.1016/ j.ijggc.2014.01.006. [35] Geistlinger H, Krauss G, Lazik D, Luckner L. Direct gas injection into saturated glass beads: transition from incoherent to coherent gas flow pattern. Water Resour Res 2006;42:W07403. http://dx.doi.org/10.1029/2005WR004451. [36] Vogel H-J, Weller U, Schlueter S. Quantification of soil structure based on Minkowski functions. Comput Geosci 2010;36:1236–45. http://dx.doi.org/ 10.1016/j.cageo.2010.03.007. [37] Rudin LI, Osher S, Fatemi E. Nonlinear total variation based noise removal algorithms. Phys D 1992;60:259–68. http://dx.doi.org/10.1016/01672789(92)90242-F. [38] Kapur JN, Shao PK, Wong AKC. A new method for gray-level picture thresholding using the entropy of the histogram. Comput Vision Graphics Image Process 1985;29:273–85. http://dx.doi.org/10.1016/0734-189X(85) 90125-2. [39] Vogel HJ, Kretzschmar A. Topological characterization of pore space in soil sample preparation and digital image-processing. Geoderma 1996;73:23–38. http://dx.doi.org/10.1016/0016-7061(96)00043-2.

[40] Hoshen J, Kopelman R. Percolation and cluster distribution. I. Cluster multiple labeling technique and critical concentration algorithm. Phys Rev B 1976;14:3438–45. http://dx.doi.org/10.1103/PhysRevB.14.3438. [41] Pan C, Dalla E, Franzosi D, Miller CT. Pore-scale simulation of entrapped nonaqueous phase liquid dissolution. Adv Water Resour 2007;30:623–40. http:// dx.doi.org/10.1016/j.advwatres.2006.03.009. [42] Bryant S, King PR, Mellor DW. Network model evaluation of permeability and spatial correlation in a real random sphere packing. Transp Porous Media 1993;11:53–70. http://dx.doi.org/10.1007/BF00614635. [43] Bryant S, Blunt M. Prediction of relative permeability in simple porous media. Phys Rev A 1992;46:2004–11. http://dx.doi.org/10.1103/PhysRevA.46.2004. [44] Fisher ME. Phys B 1967; 255–61. [45] Wilkinson D. Percolation effects in immiscible displacement. Phys Rev A 1986;34:1380–91. http://dx.doi.org/10.1103/PhysRevA.34.1380. [46] Blunt MJ, King MJ, Scher H. Simulation and theory of two-phase flow in porous media. Phys Rev A 1992;46:7680. http://dx.doi.org/10.1103/PhysRevA. 46.7680. [47] Joseph DD, Renardy YY. Fundamentals of two-fluid dynamics, Part II: lubricated transport, drops and miscible liquids. Interdisciplinary applied mathematics series, vol. 4. New York: Springer; 1993.