Regular Article
New insights into tracer propagation in partially saturated porous media
IFP Energies nouvelles, 1 et 4 avenue de BoisPréau, 92852 RueilMalmaison Cedex, France
^{*} Corresponding authors: vlasios.leontidis@ifpen.fr; daniela.bauer@ifpen.fr
Received:
9
July
2019
Accepted:
26
March
2020
This work deals with the influence of partial saturation on the transport process of a passive tracer. Transport experiments were done in a waterwet 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 spatiotemporal 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 nonGaussian nature of the transport.
© V. Leontidis et al., published by IFP Energies nouvelles, 2020
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 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 socalled 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:(1)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 σ^{2} is proportional to time. Thus, D can be estimated as .
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–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 twophase 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 highlevel nonintrusive methods is now available addressing different variables. Nuclear Magnetic Resonance Pulsed Field Gradient (NMRPFG) yields distributions of displacements and residence times for fluid particles [9], whereas tracer experiments followed by e.g. XRay (spectrometry, tomography, etc.) [10, 11] and NMR [12] provide timedependent 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–18]. All research groups observed nonGaussian 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 nonwetting 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 spatiotemporal behavior of the transport. Consequently, they constrain the fitting procedure better than BTCs.
2 Methodology
2.1 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.
Fig. 1 (a) Experimental setup, (b–d) injection sequence for the imbibition process. 
2.2 Porous media
All tests were performed on a commercial 2D microfluidic glass (waterwet) chip (Micronit, Fig. 2). The geometry of the 2D porous medium was designed based on microCT 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 μm. The mean pore size is d_{p} = 250 μm, the porosity ϕ = 0.57, the pore volume V_{p} = 2.3 μL, and the permeability k = 2.5 D. Figure 2c shows the longitudinal porosity profile.
Fig. 2 (a) Commercial microfluidic chip, (b) representation of the porous media, (c) porosity profile. 
2.3 Fluids
MilliQ water and ndodecane (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: ρ_{w,b*} = 1 g/cm^{3}, μ_{w,b*} = 0.98 mPa s, ρ_{o} = 0.75 g/cm^{3} and μ_{o} = 1.42 mPa s, respectively. The interfacial tension between the oil and the background solution found γ = 31 mN/m. The higher ink concentration C_{o*} was chosen in a way that ρ_{w,o*} ≈ ρ_{w,b*}, μ_{w,o*} ≈ μ_{w,o*}, while simultaneously providing sufficient optical contrast.
2.4 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 fourway 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 premixing 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 μL/min; 1.5 μL/min} were used resulting in water saturations S_{w} between {0.73; 1}.
2.5 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 μm/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 (BeerLambert 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.
2.6 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:(2)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).
Fig. 3 Calibration curve for the image technique. 
2.7 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 (8bit RGB color) were initially converted into an 8bit 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.
2.8 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 onedimensional analytical solutions exist [20], depending on the initial and boundary conditions. We consider the following system, corresponding to continuous injection:(3)
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:(4)
3 Results and discussion
3.1 Saturation profiles
The Representative Elementary Volume (REV) of the saturation is defined as the minimal volume above which the saturation S_{w}(Δx, Δy) computed inside the box of size (ΔxΔy) 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} (Δx = l, Δy = l) for increasing values of the box of size (ΔxΔy = l^{2}), where l is the length scale of the each side of the box. For higher water saturations, S_{w}(Δx, Δy) 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}(Δx, Δy) doesn’t reach S_{w}.
Fig. 4 REV saturation calculation (the dashed lines are the corresponding average saturation in the whole micromodel). 
The temporal and spatial evolution of the local water saturation can reveal whether the oil distribution in the micromodel is modified by the tracer injection. The longitudinal (S_{wx}) and transversal (S_{wy}) saturation profiles (Figs. 5a and 5b for S_{w} = 0.85) remain relatively constant over the whole duration of the experiment. Slight fluctuations can be attributed to the compression or extension of oil ganglia due to small variations in the pressure field. Thus, the initially created spatial oil/water distribution is not altered. In all cases the temporal variance of the saturation, defined as , where 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.
Fig. 5 (a) and (b) Water saturation profiles for partially saturated conditions (S_{w} = 0.85) for different time steps and comparison with average saturation. 
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].
Fig. 6 Oil ganglia size distribution for different water saturations. 
Fig. 7 Tracer concentration fields in full and partial saturation. 
3.2 Concentration fields at different saturations
Figure 7 shows the tracer concentration distribution in the micromodel at full (S_{w} = 1) and partial saturation (S_{w} = 0.96, 0.85, 0.82, 0.73) at τ = 1, τ = 1.3 and at the end of the experiment, where the normalized time is 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.
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.
Fig. 8 Concentration field at S_{w} = 0.73 showing a stagnant zone where the tracer only propagates by diffusion (black square) and a stagnant zone that is only accessible by water films (red square and arrow). 
Considering the deadend 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 deadend zone as , where N is the total number of pixels and V_{dz} is the volume of the deadend zones. Figure 9 presents the evolution in time of the tracer volume entering in the zones, which is equal to ΔV_{t,dz} = V(t)−V(t = 0), normalized by the corresponding pore volume to take into account the difference in the two pores area.
Fig. 9 Temporal evolution of the tracer volume entering in the deadzones for S_{w} = 0.77. 
3.3 Tracer propagation zones
Zones containing tracer are defined as the areas where the local normalized concentration is [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 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 nonaccessible during the experimental time.
Fig. 10 Area of zones containing tracer as a function of τ. 
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 nonnegligible.
Fig. 11 as a function of the water saturation: (a) S_{1}, (b) S_{2}. 
3.4 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}(τ), to the breakthrough curves, C_{BTC}(τ), (corresponding to the average concentration in the outlet pixels) and to the spatial concentration profiles, C(x, τ), 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}(τ) and C_{BTC}(τ) and the corresponding fit with equation (4) of C_{avg}(τ) and C_{BTC}(τ). It can be seen that C_{avg}(τ) and C_{BTC}(τ) 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).
Fig. 12 Black line: temporal evolution of the average concentration in the micromodel C_{avg}(τ) and the breakthrough curves C_{BTC}(τ); red line: corresponding fit using equation (4). 
However, when fitting equation (4) to the spatial concentration profiles, C(x, τ), (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).
Fig. 13 Spatial concentration profiles at different times, black line: corresponding fit using 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}(τ), C_{BTC}(τ) (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}(τ) and C_{BTC}(τ) slightly overestimate those computed from the spatial concentration profiles. Also, dispersion coefficients obtained from the spatial concentration profiles are not constant and may show an evolution in time.
Fig. 14 Normalized dispersion coefficients as a function of the saturation. 
Figure 15 shows for S_{w} = 1 and S_{w} = 0.85. Whereas for S_{w} = 1, remains relatively constant, increases in time for S_{w} = 0.85. The fact that for S_{w} = 0.85 is time dependent indicates a nonFickian behavior and hence anomalous transport in some partially saturated porous media. From the fact that is time dependent it can be stated that the evolution of the tracer variance is not proportional to time but varies with σ ≈ 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 nonGaussian 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, nonGaussian dispersion is favored by lower water saturations (S_{w} ≤ 0.85).
Fig. 15 Temporal evolution of the dispersion coefficients obtained from the spatial concentration profiles for partial and full saturation. 
4 Conclusion
This work deals with the influence of partial saturation on the transport of a passive tracer. Experiments were done using a waterwet glass micromodel allowing the simultaneous observation of the immiscible phases and the tracer propagation. Using a waterwet 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 socalled Mobile–ImmobileModels [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 NMRPFG experiments allowing a direct correlation between transport and immiscible twophase 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 spatiotemporal 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 nonGaussian nature of the transport.
References
 Gouze P., Le Borgne T., Leprovost R., Lods G., Poidras T., Pezard P. (2008) Nonfickian dispersion in porous media: 1. Multiscale measurements using singlewell injection withdrawal tracer tests, Water Resour. Res. 44, 6. [Google Scholar]
 Charlaix E., Hulin J.P., Plona T.J. (1987) Experimental study of tracer dispersion in sintered glass porous materials of variable compaction, Phys. Fluids 300, 6, 1690–1698. [CrossRef] [Google Scholar]
 Nützmann G., Maciejewski S., Joswig K. (2002) Estimation of water saturation dependence of dispersion in unsaturated porous media: Experiments and modelling analysis, Adv. Water Resour. 250, 5, 565–576. [Google Scholar]
 Maraqa M.A., Wallace R.B., Voice Th.C. (1997) Effects of degree of water saturation on dispersivity and immobile water in sandy soil columns, J. Contam. Hydrol. 250, 3–4, 199–218. [Google Scholar]
 Haga D., Niibori Y., Chida T. (1999) Hydrodynamic dispersion and mass transfer in unsaturated flow, Water Resour. Res. 350, 4, 1065–1077. [Google Scholar]
 Padilla I.Y., Yeh T.C.J., Conklin M.H. (1999) The effect of water content on solute transport in unsaturated porous media, Water Resour. Res. 350, 11, 3303–3313. [Google Scholar]
 Vanderborght J., Vereecken H. (2007) Review of dispersivities for transport modeling in soils, Vadose Zone J. 60, 1, 29–52. [Google Scholar]
 Birkholzer J., Tsang C.F. (1997) Solute channeling in unsaturated heterogeneous porous media, Water Resour. Res. 330, 10, 2221–2238. [Google Scholar]
 Guillon V., Fleury M., Bauer D., Néel M.C. (2013) Superdispersion in homogeneous unsaturated porous media using nmr propagators, Phys. Rev. E 870, 4, 043007. [Google Scholar]
 Latrille C., Néel M.C. (2013) Transport study in unsaturated porous media by tracer experiment in a dichromatic Xray experimental device, in: EPJ Web of Conferences, EDP Sciences, Vol. 50, p. 04002. [CrossRef] [EDP Sciences] [Google Scholar]
 Fourar M., Radilla G. (2009) Nonfickian description of tracer transport through heterogeneous porous media, Transp. Porous Media 800, 3, 561–579. [Google Scholar]
 Guillot G., Kassab G., Hulin J.P., Rigord P. (1991) Monitoring of tracer dispersion in porous media by nmr imaging, J. Phys. D: Appl. Phys. 240, 5, 763. [CrossRef] [Google Scholar]
 de Anna P., JiménezMartnez J., Tabuteau H., Turuban R., Le Borgne T., Derrien M., Méheust Y. (2013) Mixing and reaction kinetics in porous media: An experimental pore scale quantification, Environ. Sci. Technol. 480, 1, 508–516. [Google Scholar]
 JiménezMartnez J., de Anna P., Tabuteau H., Turuban R., Le Borgne T., Méheust Y. (2015) Porescale mechanisms for the enhancement of mixing in unsaturated porous media and implications for chemical reactions, Geophys. Res. Lett. 420, 13, 5316–5324. [Google Scholar]
 JiménezMartnez J., Le Borgne T., Tabuteau H., Méheust Y. (2017) Impact of saturation on dispersion and mixing in porous media: Photobleaching pulse injection experiments and shearenhanced mixing model, Water Resour. Res. 530, 2, 1457–1472. [Google Scholar]
 Karadimitriou N.K., JoekarNiasar V., Babaei M., Shore C.A. (2016) Critical role of the immobile zone in nonfickian twophase transport: A new paradigm, Environ. Sci. Technol. 500, 8, 4384–4392. [Google Scholar]
 Karadimitriou N.K., JoekarNiasar V., Brizuela O.G. (2017) Hydrodynamic solute transport under twophase flow conditions, Sci. Rep. 70, 1, 6624. [Google Scholar]
 Hasan S., JoekarNiasar V., Karadimitriou N.K., Sahimi M. (2019) Saturation dependence of nonfickian transport in porous media, Water Resour. Res. 550, 2, 1153–1166. [Google Scholar]
 Karadimitriou N.K., Hassanizadeh S.M. (2012) A review of micromodels and their use in twophase flow studies, Vadose Zone J. 110, 3. [Google Scholar]
 Van Genuchten M.Th. (1980) Determining transport parameters from solute displacement experiments, U.S. Dept. of Agriculture, Science and Education Administration, U.S. Salinity Laboratory, Riverside, CA. [Google Scholar]
 Oughanem R., Youssef S., Bauer D., Peysson Y., Maire E., Vizika O. (2015) A multiscale investigation of pore structure impact on the mobilization of trapped oil by surfactant injection, Transp. Porous Media 109, 673–692. [Google Scholar]
 Van Genuchten M.Th., Wierenga P.J. (1976) Mass transfer studies in sorbing porous media i. analytical solutions 1, Soil Sci. Soc. Am. J. 400, 4, 473–480. [Google Scholar]
 Haggerty R., Gorelick S.M. (1995) Multiplerate mass transfer for modeling diffusion and surface reactions in media with porescale heterogeneity, Water Resour. Res. 310, 10, 2383–2400. [Google Scholar]
 Néel M.C., Bauer D., Fleury M. (2014) Model to interpret pulsedfieldgradient nmr data including memory and superdispersion effects, Phys. Rev. E 890, 6, 062121. [Google Scholar]
 Fleury M., Bauer D., Néel M.C. (2015) Modeling of superdispersion in unsaturated porous media using nmr propagators, Micropor. Mesopor. Mater. 205, 75–78. [CrossRef] [Google Scholar]
All Figures
Fig. 1 (a) Experimental setup, (b–d) injection sequence for the imbibition process. 

In the text 
Fig. 2 (a) Commercial microfluidic chip, (b) representation of the porous media, (c) porosity profile. 

In the text 
Fig. 3 Calibration curve for the image technique. 

In the text 
Fig. 4 REV saturation calculation (the dashed lines are the corresponding average saturation in the whole micromodel). 

In the text 
Fig. 5 (a) and (b) Water saturation profiles for partially saturated conditions (S_{w} = 0.85) for different time steps and comparison with average saturation. 

In the text 
Fig. 6 Oil ganglia size distribution for different water saturations. 

In the text 
Fig. 7 Tracer concentration fields in full and partial saturation. 

In the text 
Fig. 8 Concentration field at S_{w} = 0.73 showing a stagnant zone where the tracer only propagates by diffusion (black square) and a stagnant zone that is only accessible by water films (red square and arrow). 

In the text 
Fig. 9 Temporal evolution of the tracer volume entering in the deadzones for S_{w} = 0.77. 

In the text 
Fig. 10 Area of zones containing tracer as a function of τ. 

In the text 
Fig. 11 as a function of the water saturation: (a) S_{1}, (b) S_{2}. 

In the text 
Fig. 12 Black line: temporal evolution of the average concentration in the micromodel C_{avg}(τ) and the breakthrough curves C_{BTC}(τ); red line: corresponding fit using equation (4). 

In the text 
Fig. 13 Spatial concentration profiles at different times, black line: corresponding fit using equation (4). 

In the text 
Fig. 14 Normalized dispersion coefficients as a function of the saturation. 

In the text 
Fig. 15 Temporal evolution of the dispersion coefficients obtained from the spatial concentration profiles for partial and full saturation. 

In the text 