New insights into tracer propagation in partially saturated porous media

This work deals with the influence of partial saturation on the transport process of a passive tracer. Transport experiments were done in a water-wet glass micromodel combined with specific optical techniques. Full water saturation was achieved by injecting initially the background solution and then the tracer, whereas for the partial saturation conditions, the micromodel was initially saturated with oil, and then sequential the background solution and the tracer were injected at the same flow rate. We have shown that in the investigated range of water saturations it exists a transition in the oil ganglia structure and size. For high water saturations oil ganglia have one or two pores in size, however for lower water saturations they comprise an important number of pores. Transport strongly depends on the size distribution of the oil ganglia as they create large percolating paths and stagnant zones. We also showed the existence of two different types of stagnant zones: zones accessible by diffusion into pores and zones only accessible by spatially limited diffusion in films. The major advantage of using glass micromodels lies in the fact that dispersion coefficients can be computed from concentrations averaged over the pore space or from concentrations at the outlet and simultaneously from spatial concentration profiles. Curves were fitted using the Advection–Dispersion Equation (ADE) with adequate boundary conditions. The fitting quality of the temporal evolution of the average and outlet concentration was very good. However, fitting of the concentration profiles could only be done for the higher water saturations. This is due to the fact that the Representative Elementary Volume (REV) of lower water saturations is larger than the micromodel. The results show that fitting the breakthrough curve in order to determine the dispersion coefficient in a partially saturated porous medium might be misleading. Indeed, when fitting the breakthrough curves we were able to compute a dispersion coefficient even in the case where the REV of the water saturation is larger than the micromodel. Consequently, the knowledge of the local concentration profiles as a function of time is necessary as it provides an additional information on the spatio-temporal behavior of the transport process and therefor a supplementary constraint of the fitting procedure. Finally, we observed a time dependent dispersion coefficient in the regime where oil ganglia comprise several pores. This fact might be attributed to the non-Gaussian nature of the transport.


Introduction
Transport in partially saturated porous media is a process that has gained attention in the last decades because it is encountered in diverse energy and environmental applications and processes, such as oil recovery, CO 2 storage, or transport of pollutants in the vadose zone. The distribution of the two immiscible phases in the porous medium strongly affects transport properties of molecules miscible in one of the phases. A slight modification in the partial saturation affects the relative organization of the immiscible phases and therefore transport, making the comprehension of the latter an important issue.
Dispersion in porous media is often studied by means of tracer tests and the analysis of the temporal evolution of the tracer concentration at the outlet of a core sample (the so-called breakthrough curve). In a homogeneous porous medium the resulting curve can be analyzed using the Advection-Dispersion Equation (ADE) to obtain the dispersion coefficient. The latter equation is given by: where C is the concentration of the solute, u is the average pore velocity and D is the dispersion coefficient. x and t correspond to the spatial and temporal coordinates. The ADE is based on solute mass conservation, assuming that Fick's first law is valid. In the case of finite volume injection, the displacement distribution of the tracer corresponds to a Gaussian distribution and the evolution of its second moment r 2 is proportional to time. Thus, D can be estimated as D ¼ lim t!1 1 2 d dt r 2 . However, experimental data have shown that tracer transport may not always follow the ADE (for example [1,2]). In some types of porous media, the second moment evolves more rapidly as if some part of the tracer was transported faster than according to the ADE. In other media, important tailing in the breakthrough curves is observed, which is also not compatible with the above mentioned law. In these cases dispersion is considered to be anomalous and Fick's law fails to describe the related transport behavior.
Studying transport in partially saturated porous media, where tracer molecules are transported in one of the two phases, remains a highly challenging task and results on the dependence of the dispersion coefficient on saturation found in literature [3][4][5][6][7][8] are not always consistent. The major difficulty lies in the fact that the pore volume accessible by the tracer strongly depends on the physics of the underlying immiscible two-phase flow. Indeed, the actual spatial distribution of the two immiscible phases strongly depends on the contact angle between the rock and the wetting fluid, the interfacial tension between the two phases but also on the mean interstitial fluid velocity. Also, due to the wettability of the porous medium wetting films can be observed. Additionally, the pore size distribution as well as the drainage and imbibition cycle strongly affects the immiscible phase distribution. In some situations, phases might become badly connected leading to displacement distributions of the tracer being very different from the classical model. In this case dispersion coefficients may become time and/or space dependent. This indicates, that in the asymptotic regime, normal transport is not an appropriate model and anomalous transport should be considered. As a consequence the simple acquisition of the breakthrough curve classically applied is often not sufficient for a correct prediction of the transport and the temporal evolution of the spatial concentration profiles is desirable. An ensemble of high-level non-intrusive methods is now available addressing different variables. Nuclear Magnetic Resonance Pulsed Field Gradient (NMR-PFG) yields distributions of displacements and residence times for fluid particles [9], whereas tracer experiments followed by e.g. X-Ray (spectrometry, tomography, etc.) [10,11] and NMR [12] provide time-dependent spatial concentration profiles.
Recently, several research groups started using 2D micromodels in combination with specific optical methods to study tracer propagation in porous media at full [13] and partial saturation [14][15][16][17][18]. All research groups observed non-Gaussian transport and related them to the presence of stagnant zones in the partially saturated porous medium. However, they didn't consider delayed transport due to diffusion in films. Glass micromodels allow the direct visualization of the spatial distribution of the immiscible phases and of the tracer. Based on the time evolution of the tracer concentration, the temporal variation of its displacement variance can be determined in order to characterize dispersion [15]. Micromodels also allow the distinction between stagnant and mobile zones [16], giving further information on the coupling between saturation and dispersion. A review on micromodels, their fabrication, visualization techniques, and applications can be found in [19].
The objective of the present article is to characterize transport phenomena in fully and partially saturated porous media and particularly in the case where water films line the solid grains. We will therefore further investigate tracer transport into stagnant zones in comparison with transport into zones accessible only via films. We will also focus on the distribution of the non-wetting fluid depending on the saturation and its influence on the transport. To this goal, we are comparing dispersion coefficients obtained from the BreakThrough Curves (BTC) and those obtained from the temporal evolution of the spatial concentration profiles, as micromodel experiments allow the acquisition of both. The spatial profiles naturally provide more information as they contain the entire spatio-temporal behavior of the transport. Consequently, they constrain the fitting procedure better than BTCs.

Experimental setup
In order to visualize the tracer distribution in the saturated or partially saturated micromodel the experimental setup shown in Figure 1a is used. A detailed description of the individual parts is given below.

Porous media
All tests were performed on a commercial 2D microfluidic glass (water-wet) chip (Micronit, Fig. 2). The geometry of the 2D porous medium was designed based on micro-CT images of a sandstone. Channels were created by isotropic etching of a glass plate. One inlet and one outlet exist on the chip. The length and the width of the porous area are L = 20 mm, W = 10 mm, with the depth h = 20 lm. The mean pore size is d p = 250 lm, the porosity / = 0.57, the pore volume V p = 2.3 lL, and the permeability k = 2.5 D. Figure 2c shows the longitudinal porosity profile.

Fluids
Milli-Q water and n-dodecane (VWR Chemical) were used as immiscible fluids. The water was dyed with ink (Parker Quink Blue/black) at two concentrations: a low concentration, C b* , at 3% w/w was used as background solution in order to enhance the contrast between the fluid and the glass structure. A high concentration, C o* , at 40% w/w was used as tracer. Density and viscosity of the aqueous solutions at low ink concentration and of the oil measured: q w,b* = 1 g/cm 3 , l w,b* = 0.98 mPa s, q o = 0.75 g/cm 3 and l o = 1.42 mPa s, respectively. The interfacial tension between the oil and the background solution found c = 31 mN/m. The higher ink concentration C o* was chosen in a way that q w,o* % q w,b* , l w,o* % l w,o* , while simultaneously providing sufficient optical contrast.

Fluid displacements
Flow is established by means of a low pressure syringe pump neMESYS 290N with two syringes filled with the dyed water at C b* and C o* . As the microfluidic chip is equipped with only one inlet, a four-way valve was placed as close as possible to the inlet of the chip to ensure that the two fluids with different ink concentrations are in contact only over a short distance and time before entering the porous medium, minimizing the pre-mixing and dispersion in the inlet tubing. Oil is injected using a plastic syringe connected to the chip through the valve. The direction of flow in the micromodel is opposite to gravity.
For full saturation, S w = 1, the micromodel was initially saturated with the background solution (Fig. 1c). Then, the tracer (at C o* ) was injected at the desired flow rate, Q (Fig. 1d). For tests under partially saturated conditions, the micromodel was initially saturated with oil ( Fig. 1b), and partial saturation was obtained by injecting the background solution at a specific flow rate, Q, until the water saturation remained constant (Fig. 1c). This corresponds to an imbibition process. Then, the tracer was injected at the same flow rate, in order to preserve the balance between the two immiscible phases (Fig. 1d). After each tracer experiment the background solution was injected at the same flowrate to clean the tracer from the micromodel. Then the experiment was performed at the next higher flow rate. Flow rates in the range of {0.1 lL/min; 1.5 lL/min} were used resulting in water saturations S w between {0.73; 1}.

Visualization technique
The experimental setup was built on an optical table. The chip holder was placed between a red LED source (Phlox) and a high speed camera (maximum frame rate: 160 fps, image size: 2000 Â 1086 pixels). A 90 mm macro Olympus objective was used resulting in a spatial resolution of approximately 11.4 lm/pixel corresponding to an effective image size of the porous medium containing 1764 Â 882 pixels. The micromodel is illuminated from behind with a spatially homogenous light intensity emitting at 625 nm (±15). The visualization method is based on the absorption of the emitted light by a sample. The absorbency depends both on the solution concentration and the sample thickness (Beer-Lambert law). The solutions inside the glass micromodel are excited by the light source and emit at the same wavelength without the need for optical filters.

Calibration of the visualization technique
In the case of the dyed water, it is necessary to correlate the emitted light not absorbed by the fluid and captured by the high speed camera with the ink concentration. The calibration method consists in saturating the porous medium with solutions at different ink concentrations in the working range (0-40% w/w) and taking an image under full saturation. From the images acquired at every concentration, the average grey level, I, in the pore space was calculated and associated with the ink concentration (Fig. 3). The calibration data were fitted by: where a, b, and c are the calibration constants depending on the specific solutions and the optical setup (light intensity, objective focus, distance, and position of the elements).

Image acquisition and processing
Based on the applied flow rates and the micromodel dimensions it was sufficient to acquire images at a frame rate of 5 fps. All raw images (8-bit RGB color) were initially converted into an 8-bit greyscale using the ImageJ software, resulting in 256 maximum different concentrations. The grey level of each pixel was then converted into an ink concentration value. Data analysis was performed at 1 fps.

Data analysis
Average concentrations over the entire pore space were calculated before the tracer injection (C b ) and at the end of the experiment (C o ). In the case of partial saturation, C o calculated from the image processing is lower than the injected ink concentration C o* . This is due to the fact that there are always pores which are difficult to be accessed by the tracer. Consequently, reaching maximum concentration in these zones takes more time than the experiment lasts. Thus, these zones reduce the final concentration C o . For equation (1) several one-dimensional analytical solutions exist [20], depending on the initial and boundary conditions. We consider the following system, corresponding to continuous injection: The latter equation stands for an infinitely long system where the solute is far from reaching the outlet. However, this approximation is often used for the interpretation of tracer experiments. The analytical solution of equation (1) considering these boundary conditions is then given by: 3 Results and discussion

Saturation profiles
The Representative Elementary Volume (REV) of the saturation is defined as the minimal volume above which the saturation S w (Dx, Dy) computed inside the box of size (DxDy) becomes independent of the box size. Computing the REV based on the saturation allows the distinction between a porous medium whose structure is sufficiently homogeneous for a potential upscaling of its transport properties and a porous medium whose heterogeneities are too important in comparison to the size of the micromodel to allow upscaling. Figure 4 shows S w (Dx = l, Dy = l) for increasing values of the box of size (DxDy = l 2 ), where l is the length scale of the each side of the box. For higher water saturations, S w (Dx, Dy) reaches a relatively constant value for large box sizes, that is close to the saturation value S w of the entire micromodel. For lower saturations, S w (Dx, Dy) doesn't reach S w . The temporal and spatial evolution of the local water saturation can reveal whether the oil distribution in the  Thus, the initially created spatial oil/water distribution is not altered. In all cases the temporal variance of the saturation, defined as r S w ¼ 1 S w,t is the instantaneous spatial average saturation over the whole micromodel and S w is the temporal and spatial mean saturation over the time interval (N t ), is on the order of 10 À4 which supports the assumption of the stable oil clusters.
Oil cluster size, expressed by the 2D surface of the ganglia A clusters , and their size distribution depend on the imposed flow rate. Figure 6 shows the oil cluster size distribution for different water saturations. As can be seen it exists a transition in the oil ganglia structure. Indeed, for high water saturations a large number of small ganglia of the size of one or two pores were observed. For lower saturations there are less oil ganglia comprising a large number of pores. This is also obvious on the field images on Figure 7. The cluster size distribution follows a power law with a slight deviation for larger ganglia sizes for all saturations. Similar results have been observed in [21]. with t h being the time for each experiment that the tracer requires to propagate in half of the pore space containing only water before the tracer injection.

Concentration fields at different saturations
At full saturation, the solute propagates gradually and almost uniformly in the porous medium. Only very few regions, where the tracer propagates with a lower rate can be seen. Also, structural heterogeneities of the micromodel lead to an inhomogeneous repartition of the tracer at the entrance, that affects the tracer front. In the case of complete saturation this effect is compensated by the homogeneous distribution of the water phase, however it is amplified in the case of partial saturation. Here, tracer propagation is more heterogeneous and the presence of oil ganglia creates additional preferential paths and stagnant zones. For S w = 0.96 some stagnant zones exist, but their number is very limited as oil ganglia have the size of the pores. However, the rise in ganglia size with the decreasing water saturation leads to a larger number of stagnant zones, but also to zones that are only accessible by water films. Figure 8 shows both specific cases where tracer propagation is slowed down. Both, preferential paths and stagnant zones create a widespread velocity distribution that is different from the one observed at full saturation. Thus tracer propagation becomes very heterogeneous.
Considering the dead-end zones of Figure 8, the volume of the tracer diffusing into the stagnant pore or invading by film at each time step can be calculated from the average tracer concentration over the total volume of the dead-end zone as V tr;dz ðtÞ ¼ where N is the total number of pixels and V dz is the volume of the dead-end zones. Figure 9 presents the evolution in time of the tracer volume entering in the zones, which is equal to DV t,dz = V(t)ÀV (t = 0), normalized by the corresponding pore volume to take into account the difference in the two pores area.

Tracer propagation zones
Zones containing tracer are defined as the areas where the local normalized concentration is C eff ðx; y; tÞ ¼ C ðx;y;tÞÀC bÃ C oÃ ÀC bÃ > 0:05 [14]. S prop is defined as the area where the tracer actually propagates normalized by the area of the total pore space [16]. Figure 10 shows the temporal evolution of S prop for S w = 1, 0.85 and 0.73. In all cases two distinct slopes dSprop ds can be observed. The first slope (S 1 ) corresponds to the propagation of the tracer in the mobile zones, whereas the second slope (S 2 ) corresponds to the invasion of stagnant zones or zones accessible only by films. Only in the case of full saturation does the tracer propagate at the end of the experiment in the entire pore space. Even at high water saturation oil ganglia have formed stagnant zones non-accessible during the experimental time. Figure 11 gives (a) S 1 and (b) S 2 as a function of S w . Considering S 1 , two distinct regions can be seen. For high water saturations (S w > 0.85), due to the limited number of data, no specific correlation between the water saturation and the slopes can be established. However, a relatively linear relation between these parameters can be observed for lower water saturations. S 2 increases with decreasing water saturations. This increase is due to tracer propagation into stagnant zones by molecular diffusion and particularly into zones only accessible by films. As the number of stagnant zones increases with decreasing water saturation transport kinetics into these zones become non-negligible.

Determination of the dispersion coefficient
It is common in tracer experiments to analyze the concentration at the outlet and to compute the dispersion coefficient. Glass micromodels together with optical techniques provide the opportunity of real time measurements over the whole media. Comparing the breakthrough curves at the outlet of the micromodel or the temporal evolution of the average concentration in the entire micromodel to the spatial concentration profiles provides further information on the transport in partially saturated porous media.
Dispersion coefficients were determined by fitting equation (4) to the temporal evolution of the average concentration in the micromodel, C avg (s), to the breakthrough curves, C BTC (s), (corresponding to the average concentration in the outlet pixels) and to the spatial concentration profiles, C(x, s), using the velocity as additional fitting parameter. The velocity, which depends on the actual water saturation, the spatial distribution of the oil ganglia and the imposed by the pump flow rate, was chosen as a free parameter in order to account for the fluctuations of the flow velocity in the pores resulting either from the continuous change of the oil ganglia shape (due to pressure fluctuations) or from a potential error in the imposed flow rate mainly at the lower flow rates and water saturations where the size of the oil clusters is large. Unfortunately, volume changes cannot be visualized and measured since the experimental technique provides 2D images of a 3D system (depth of the micromodel). Figure 12 shows C avg (s) and C BTC (s) and the corresponding fit with equation (4) of C avg (s) and C BTC (s). It can be seen that C avg (s) and C BTC (s) are similar except from the time shift, that corresponds to the time the tracer  needs to reach the end of the micromodel. Good fitting was achieved even for low water saturations (R 2 = 0.9972, coefficient of determination).
However, when fitting equation (4) to the spatial concentration profiles, C(x, s), (Fig. 13), fitting quality is not identical. Very good fitting was obtained for full saturation (R 2 % 0.99), however particularly for lower water saturations the fitting quality is very poor (S w = 0.73, R 2 = 0.81 for the best fit). This can be explained by the size of the REV of the water saturation which is larger than the size of the micromodel, which demonstrates the strong heterogeneity of the phase distribution at this lower water saturation. Thus, the results show that fitting the breakthrough curve in order to determine the dispersion coefficient in a partially saturated porous medium might be misleading as equation (4) perfectly fits the breakthrough curve even if the size of the micromodel is smaller than the REV. However the spatial and temporal evolution of the concentration doesn't follow equation (4).
Considering S w = 0.85, 0.96, and 1 the fitting quality is sufficiently good (R 2 % 0.93, %0.96, %0.99) to determine the dispersion coefficients from the spatial concentration profiles. Values of the normalized dispersion coefficients D* are given in Figure 14 obtained from C avg (s), C BTC (s) (all saturations) and from the spatial concentration profiles (S w = 0.85, 0.96, 1). In the saturation range 0.85 < S w < 1 dispersion coefficients increase with decreasing water saturation. This is due to the heterogeneity in the velocity induced by the partial saturation. As can be seen, values obtained from C avg (s) and C BTC (s) slightly overestimate those computed from the spatial concentration profiles. Also, dispersion coefficients D Ã x ðsÞ obtained from the spatial concentration profiles are not constant and may show an evolution in time. Figure 15 shows D Ã x ðsÞ for S w = 1 and S w = 0.85. Whereas for S w = 1, D Ã x ðsÞ remains relatively constant, D Ã x ðsÞ increases in time for S w = 0.85. The fact that D Ã x ðsÞ for S w = 0.85 is time dependent indicates a non-Fickian behavior and hence anomalous transport in some partially saturated porous media. From the fact that D Ã x ðsÞ is time dependent it can be stated that the evolution of the tracer variance is not proportional to time but varies with r % t a . The variance can also be computed using the time derivative of the data (that is to subtract each image from the next one). However, this methodology enhances the noise of the raw data. Further filtering is impossible to be applied  without losing important information. The non-Gaussian nature of the transport is due to the coupling between fast flowing percolating paths, stagnant zones and stagnant zones only accessible by water films where the tracer propagation is extremely slow (see Fig. 8). Similar results have been observed in [15]. Stagnant zones as well as percolating paths result from oil ganglia that are larger than one or two pores. Thus, non-Gaussian dispersion is favored by lower water saturations (S w 0.85).

Conclusion
This work deals with the influence of partial saturation on the transport of a passive tracer. Experiments were done using a water-wet glass micromodel allowing the simultaneous observation of the immiscible phases and the tracer propagation. Using a water-wet glass micromodel as porous medium we also focused on the influence of film flow. We have shown that in the investigated range of water saturations it exists a transition in the oil ganglia structure. Indeed, for high water saturations a large number of small ganglia of the size of one or two pores were observed. For lower saturations there are less oil ganglia comprising a large number of pores. Transport strongly depends on the oil ganglia distribution as they create large percolating paths and stagnant zones. We also showed the existence   of two different types of stagnant zones: zones accessible by diffusion into pores and zones only accessible by spatially limited diffusion in films. The latter zones are characterized by a lower invasion rate. Thus, two different transport kinetics into stagnant zones exist. Generally, the so-called Mobile-Immobile-Models [22] are used to describe transport in partially saturated porous media. The model is based on one exchange term between the mobile phase and the stagnant zones. However, the existence of wetting films influences the exchange kinetics and a supplementary term might become necessary to describe this type of transport as proposed by Haggerty and Gorelick [23]. Nevertheless, these models don't take into account the topology of the underlying immiscible cluster size distribution and sometimes fail to correctly describe the exchange between stagnant and mobile zones [16]. Another possibility could be the use of a fractional version of the ADE [24]. Parameters of these models can be determined from NMR-PFG experiments allowing a direct correlation between transport and immiscible two-phase flow [25].
The major advantage of using glass micromodels lies in the fact that dispersion coefficients can be computed from the average or outlet concentration and simultaneously from spatial concentration profiles. Curves were fitted using the ADE with adequate boundary conditions. The fitting quality of the temporal evolution of the average concentration and the breakthrough curve were very good for all saturations. However, fitting of the spatial concentration profiles could only be obtained for higher water saturations. This is due to the REV size of the saturation that strongly depends on the oil ganglia size and number. The results show that fitting the breakthrough curve in order to determine the dispersion coefficient in a partially saturated porous medium might be misleading. Indeed, when fitting the breakthrough curves we were able to compute a dispersion coefficient even in the case where the REV of the water saturation is larger than the micromodel. Consequently, the knowledge of the local concentration profiles as a function of time is necessary as it provides an additional information on the spatio-temporal behavior of the transport process and therefor a supplementary constraint on the fitting procedure. We finally showed that the dispersion coefficients obtained from the spatial concentration profiles may become time dependent for lower water saturations, proving the non-Gaussian nature of the transport.