Quantitative Visualization of Two‐Phase Flow in a Fractured Porous Medium

Two‐phase fluid flow in fractured porous media impacts many natural and industrial processes but our understanding of flow dynamics in these systems is constrained by difficulties measuring the exchange of water between fracture and adjacent porous matrix. We present a novel experimental system that allows quantitative visualization of the air and water phases in a single analog fractured porous medium. The fracture system consists of a sintered‐glass porous plate in contact with an impermeable glass plate. A reservoir connected to the porous plate allows control of pore pressure within the porous medium. The fracture fills and drains through the porous matrix and flow manifolds along two edges of the fracture. The fracture is mounted in an imaging system that includes a controlled light‐emitting diode panel and a charge‐coupled‐device camera. Flow and pressure are controlled and monitored by a computer during experiments. To demonstrate this system, we carried out a series of cyclic drainage and imbibition experiments in fractures bounded by porous media with different pore‐size distributions in the porous matrix. Images of the drainage process demonstrate that the air‐water distribution within the fracture evolves differently than has been observed in non‐porous fractured systems. Specifically, we observed limited trapping of water within the fracture during drainage. Conversely, during imbibition, because air cannot exit through the porous matrix, significant regions of air became entrapped once pathways to the fracture boundaries became water filled. The differences in phase evolution led to substantial differences in the evolution of estimated relative permeability with saturation.


Introduction
Two-phase flow in fractured porous media plays an important role in natural processes such as infiltration into fractured rock and engineered processes such as enhanced oil/gas recovery (Karpyn et al., 2009;Rangel-German & Kovscek, 2002), geological CO 2 sequestration (Vafaie et al., 2023), and remediation of groundwater contaminated by non-aqueous phase liquids (Dearden et al., 2013).Related two-phase flow processes can be broadly categorized as either drainage (non-wetting phase displaces the wetting phase) or imbibition (wetting phase displaces the non-wetting phase).
Early studies of two-phase flow through fractures considered fractures in a non-porous matrix.These studies included experiments in transparent models (e.g., Nicholl et al., 1994;Su et al., 1999) or replicas (e.g., Persoff & Pruess, 1995;Wan et al., 2000) and invasion percolation simulations in variable aperture fields (e.g., Glass et al., 1998;Xu et al., 1998;Yang et al., 2012).In such fractures, flow of the two fluids occurs only through the fracture and the distribution of the phases depends on characteristics of the fracture aperture and the nature and history of the displacement processes.Furthermore, because water and air enter only through the fracture edges, both the wetting and non-wetting fluid may become entrapped and immobilized in regions that are isolated from the fracture boundaries.
When the fractured matrix contains non-negligible porosity, flow of one or both fluids can occur through the fracture and the porous matrix.Two general experimental approaches have been used to study two-phase flow in fractured porous media.The first approach uses micromodels that represent a two-dimensional (2D) cross-section through a fracture; the fracture is a 2D channel and the adjacent porous matrix is a 2D slice of porous medium (e.g., Haghighi et al., 1994;Rangel-German & Kovscek, 2006;Wan et al., 1996).Such experimental systems allow direct visualization of the interaction of multiple phases within the fracture and porous matrix, but neglect the 3D interaction of the phases within the fracture induced by aperture variability.
The second approach combines two-phase flow through cores of fractured porous rock with X-ray computed tomography to observe the distribution of phases within the fracture and porous matrix (Arshadi et al., 2018;Karpyn et al., 2009;Rangel-German & Kovscek, 2002).This has the advantage of providing measurements of the distribution of two or more phases within the pore space (both fracture and porous matrix) in fractured cores.However, the temporal and spatial resolution of CT scans constrains the scalability of these studies.For example, Arshadi et al. (2018) imaged 5-mm segments of a larger fractured core at a 2.5-μm voxel resolution.Such high spatial resolution facilitates identification of phases within pores in the matrix, but inherently limits the size of the fracture that can be imaged.
We developed a novel fractured porous media experimental test cell that consists of a translucent porous glass surface mated with a transparent non-porous glass surface.This system allows imbibition and drainage experiments in which water enters and exits the fracture through the porous matrix.Quantitative visualization techniques facilitate direct measurement of the evolving phase distribution within the 15 × 15-cm fracture at a spatial resolution of ∼75 μm and a temporal resolution of ∼1 Hz.We demonstrate the new system through a series of drainage and imbibibition experiments in fractures with two different pore-size distributions in the porous matrix and corresponding differences in the distribution of fracture apertures.We present the details or our new experimental system and fabrication of the fracture model in Section 2; Section 3 describes the experimental procedure used to demonstrate the results along with the required data analyses; Section 4 presents the results of the demonstration experiments; and Section 5 provides concluding remarks.

Experimental System
The experimental system (Figure 1a) includes a test stand that rigidly supports a 12-bit camera (Quantix KAF-6303e; 2,048 × 3,072 pixels), red backlight panel (with an emitted wavelength of ∼625 nm) and the experimental model.The test stand can rotate from 90°to 90°so gravitational forces acting on the fluid in the fracture can be varied.The spatial resolution of images of the fracture plane was ∼75 × 75 μm.Opaque fabric covers the test stand to minimize stray light in the imaging system.Similar experimental systems have been used to study a range of flow and transport processes in single variable aperture fractures (e.g., Detwiler et al., 1999Detwiler et al., , 2002;;Jones & Detwiler, 2016;Nicholl et al., 1999).
In this study, we have developed a novel fracture test cell that includes a 15 × 15-cm fracture created by mating a porous glass bottom surface with a non-porous translucent glass plate top surface (Figure 1b).A unique feature of this configuration is that the bottom porous surface is both permeable and translucent.Thus, water can flow in and out of the fracture through the fracture edges and through the porous matrix and the evolving distribution of air and water within the fracture lead to variations in light transmission that can be imaged and quantified during experiments (Section 3.2).
For demonstration of this method, we chose a smooth top fracture surface.In this system, fracture aperture variability results predominantly from roughness of the porous glass plates used for the bottom fracture surface.We used two different porous glass surfaces (Rudong Shundao Glass Instrument Factory, China) with reported pore-size distributions of 16-30 μm (MF -fracture with medium pore-size matrix) and 4-7 µm (FF -fracture with fine pore-size matrix) (Table S1 in Supporting Information S1).Thus, for FF, both the mean and variance of the resulting aperture field when the surface is mated with the smooth glass surface are smaller than for MF.The characteristic length scale of aperture variability (e.g., correlation length) is also smaller for FF than MF.
Previous studies have used light transmission to measure fracture aperture by filling fractures created with two translucent non-porous surfaces (e.g., Detwiler et al., 1999Detwiler et al., , 2002;;Jones & Detwiler, 2016).Unfortunately, because the lower surface is porous in our new fracture test cell, it is not possible to fill only the fracture with dyed fluid, which is necessary to effectively use light transmission to measure fracture aperture.Thus, in this study, we use observations of phase distribution at different capillary pressures to estimate fracture apertures within each of the fractures based on the Laplace-Young equation (Glass et al., 1998).
Two 1.9-cm-thick fused-quartz windows supported by aluminum frames enclosed the fracture surfaces.Clear polyvinyl chloride gaskets separated each fracture surface from the fused-quartz window creating empty cavities between each fracture surface and the supporting window.A needle through the gasket into the lower cavity provided an inlet/outlet for water flow in/out of the porous matrix.To prevent leakage from the edges of the porous lower surface, a rim of dyed epoxy was applied along the peripheral edges of the porous glass (see Text S1, Figure S1, and Figure S2 in Supporting Information S1 for additional details).
After placing the smooth glass and the porous glass surfaces in contact to create the fracture, normal stress was applied to the frame by tightening the connecting bolts to a uniform torque (typically 1.7 Nm). Figure 1b shows a rigid frame surrounding the fracture with bolts that apply force to the no-flow boundaries (left and right sides) and flow manifolds (top and bottom).The flow manifold has a ∼3 × 5-mm channel along the entire width of the fracture to ensure that pressure gradients along the manifold channel are negligible relative to pressure gradients within the fracture.
Fluid entered and exited the fracture through tubing connected to the cavity needle and the flow manifolds.For the experiments presented here, a Mariotte bottle connected to the cavity needle served to control the head in the permeable matrix (bottom fracture surface).The Mariotte bottle was positioned on an analytical balance (Mettler Toledo MS4002S/03) on a variable-height stage, which allowed reproducible head changes.The maximum range of the stage was ±70 cm relative to the middle of the fracture plane, but we limited the range of displacements to avoid breakthrough of air into the cavity backing the porous surface.We define the capillary head as: where p a , p w , and p c are the atmospheric, water, and capillary pressure, respectively, ρ w is the density of water and g is acceleration due to gravity.During all experiments, the air pressure is equal to the atmospheric pressure p a , because the flow manifold channels, through which air enters and exits the fracture, are connected to the atmosphere, and a pressure transducer (Validyne DP15-42) monitored p w at 0.167 Hz.Mass flow rate in and out of the fracture was recorded by the analytical balance.A computer connected to the experimental system controlled data acquisition from each of the sensors (camera, pressure transducer, and balance) using Labview (e.g., Bitter et al., 2006).Manometers adjacent to the pressure transducer facilitated periodic calibration of the pressure transducer but were isolated from the fracture during drainage and imbibition experiments.Because Mariotte bottles lead to small pressure oscillations when bubbles release from the vent tube, we terminated the vent tube with a 16-gauge needle and applied a 0.4 atm vacuum to the head space in the bottle.This caused a steady stream of bubbles from the vent tube and negligible pressure oscillations.To minimize evaporation losses from the Mariotte bottle, we humidified the vent air entering the bottle (Figure 1).

Cyclic Drainage and Imbibition Experiments
To demonstrate the capabilities of this new experimental system, we carried out several cyclic drainage and imbibition experiments.Horizontal experiments were carried out in the MF and FF models (Experiments MFH and FFH, respectively) to investigate the effect of matrix pore geometry and aperture variability, and one vertical experiment was conducted in FF model (Experiment FFV) to explore the added effect of gravity.

Experimental Procedure
Each experimental sequence involved: primary drainage → primary imbibition → secondary drainage → secondary imbibition (DR1 → IMB1 → DR2 → IMB2).Before each experiment, carbon dioxide was injected through the cavity, porous matrix and dry fracture to displace air from the test cell.Then deionized, de-aired water was injected to saturate the fracture.Prior to initializing the first drainage sequence, the boundary conditions for the fracture were established.For the horizontal experiments, the flow manifolds and all connected tubing were drained so the manifolds were filled with air at atmospheric pressure.For the vertical experiment, the top flow manifold and connected tubing were drained, establishing a zero-pressure, air boundary at the top of the fracture, while the bottom manifold and connected tubing were valved off, establishing a no-flow boundary at the bottom of the fracture.We note these experiments were designed as a controlled test of the experimental system; the system can be adapted to allow a wide variety of initial and boundary conditions representing different two-phase flow processes of interest.
The drainage-imbibition cycles were conducted by sequentially varying the capillary head in the cavity through a set of static displacements of the Mariotte bottle.Image acquisition began before initially changing the position of the Mariotte bottle and continued until the final drainage or imibibition step.During each step, the valve between the cavity and Mariotte bottle (E) was closed as the height of the bottle was adjusted.The valve was then opened and the fracture was allowed to drain or fill until subsequent images indicated negligible change in phase distribution (approximately 10 min per step).Each drainage-imbibition cycle was completed during a single day followed by an approximately 12-hr pause before completing the secondary drainage-imbibition cycle.

Measurement of Evolving Phase Distribution
To aid interpretation of the images acquired during experiments, we developed an image processing script in MATLAB to convert raw images to binary images that distinguished the two phases (air/water).Figure 2 shows the steps of the image processing procedure.At the start of each experiment, we took a sequence of 100 reference images of the saturated matrix and fracture, which we averaged to yield a single, low-noise reference image (Figure 2a).To account for nonuniformities in light transmission through the porous matrix and fracture, we normalized each experimental image (Figure 2b) by the reference image.The natural logarithm of the resulting normalized field quantifies light absorbance at each pixel (Figure 2c).
Light scattering at the interface between the porous matrix and the fracture complicates differentiating air and water within the fracture.Rather than a sharp edge, the air-water interface appears in the absorbance field as a diffuse zone where the values transition from near zero for water to approximately 0.2 for air (Figure 2c).To quantify the location of the interface, we sought a global threshold that minimized the number of air clusters during the primary drainage cycle.The rationale for this approach is that, during primary drainage, we expect the formation of a connected air cluster originating from the fracture inlet with minimal fragmenting of this cluster until the beginning of the subsequent imbibition cycle.Due to noise in the images, smaller threshold values cause localized water-filled regions to be misidentified as air resulting in an increase in the number of air-filled clusters.
Larger threshold values cause some thin air-filled channels to be misidentified, separating the large invading cluster into multiple clusters.
To determine the global absorbance threshold, we developed an algorithm that sequentially incremented the threshold over a range that included the likely global threshold value, binarized the field according to each threshold value and counted the resulting number of discrete air clusters.We repeated this process for each image during the primary drainage cycle.Plotting the average number of segmented air clusters, N ave , versus the threshold values for experiments FFH and FFV, T FFH and T FFV , reveals distinct minima for these curves.The respective optimal threshold, T ⋆ , for these two experiments were T ⋆ FFH = 0.083 and T ⋆ FFV = 0.105.We selected a global value of T ⋆ = 0.094 ± 0.014 (average from the two experiments ±15%, Figure 2d) as the optimal threshold.Figures 2e and 2f show the results of binarizing using the upper and lower bounds for T ⋆ (T ⋆ lb = 0.08 and T ⋆ ub = 0.108) and demonstrate that the most significant discrepancies occur where thin tendrils of water (black phase) connect two larger regions of water.We consider the range of possible interface locations as a source of uncertainty in calculations of saturation presented in Section 4. For Experiment MFH, the optimal threshold was T ⋆ MFH = 0.042 ± 0.006 (Figure S3 in Supporting Information S1).

Experimental Results
Images acquired during each experiment allow us to quantify the evolving phase distribution within the fracture.Figure 3 compares the primary and secondary drainage and imbibition cycles for horizontal experiments in the MF and FF (MFH and FFH, respectively) and a vertical experiment in FF (FFV).The colors reflect air occupancy at sequential values of Ψ during each step of the experiment, with warm colors indicating smaller Ψ and cool colors indicating larger Ψ; black regions remained water-filled for all Ψ.The gray regions in the secondary drainage figures are regions that remained air-filled at the end of primary imbibition.
For all experiments, air entered the fracture only after Ψ exceeded the air entry pressure of the fracture (Ψ a,e ).Initial air entry occurred after the step from Ψ = 174-227 mm for MFH, Ψ = 273-325 mm for FFH, and Ψ = 195-302 mm for FFV according to the primary drainage (DR1) process (Figure 3).The different values of Ψ a,e for MF and FF reflect differences in the fracture aperture along the air inlet boundaries for the two fractures (MF and FF).Though the nonporous glass surface and porous matrix were placed in contact to create the fractures, the size of the sintered beads used to create MF were larger than those used for FF (Figure S4 in Supporting Information S1), resulting in a larger aperture and lower Ψ a,e for MF.The Laplace-Young relationship relates Ψ a,e to the corresponding fracture apertures (see Text S2 and Figure S5 in Supporting Information S1 for details) and suggests that the apertures along the flow boundaries are between 43 and 50 μm for FF and between 71 and 77 μm for MF.Note, differences in Ψ a,e for the two FF experiments (FFH and FFV) are due to using the fracture midline as the datum; the result is that the inlet in FFV is 76 mm higher than the inlet for FFH.The binarized distributions of air and water within the fracture at each increment of Ψ (Figure 3) allow us to quantify the areal fraction of the fracture occupied by water, S A w .This serves as a surrogate for volumetric water saturation, which we cannot precisely quantify because we have only capillary-pressure-derived estimates of the fracture aperture and not direct measurements of aperture variability within the fracture.Plotting areal air saturations S A a (1 S A w ) versus time for each drainage and imbibition step (Figure S6 in Supporting Information S1) demonstrates that the fracture approached equilibrium after each drainage and imbibition step, but for many steps, saturation was still changing slowly when the next pressure step was initiated.
Figures 4a-4c show Ψ plotted against S A w for each cycle of Experiments MFH, FFH, and FFV, respectively.The Ψ versus S A w plots exhibit significant hysteresis, which can be understood by comparing the phase distributions during different cycles of each experiment.As primary drainage proceeded through sequential steps of Ψ, air entered regions of progressively smaller aperture.For all experiments, air first filled regions near the air-filled flow manifolds and then advanced into the center of the fracture.For the horizontal fractures, air eventually formed a connected path between the two manifolds.Then, with further decreases in Ψ, the air-occupied region expanded toward the no-flow boundaries.Similar large-scale displacement patterns are observed in all experiments; this is likely due to the clamping pressure applied to the aluminum frames during fracture assembly (Section 2) causing smaller apertures around the perimeter of the fracture.The small-scale features of the displacement patterns differ for the and MF experiments, likely due to the difference in the porous matrix, which influences the magnitude and variability of fracture aperture.
A common feature of all experiments was the relative absence of regions of trapped water within the drained region of the fractures (i.e., isolated black regions surrounded by colored regions in Figure 3).This differs from experimental observations in fractures bounded by non-porous, impermeable surfaces, where regions of the draining phase become disconnected from the fracture edges and entrapped within the fracture (e.g., Chen et al., 2017;Detwiler et al., 2002).Here, water regions that become isolated during drainage eventually drain through the porous matrix if the aperture is large enough that Ψ exceeds the local air entry pressure in the fracture.
During primary imbibition (IMB1), water fills the smallest aperture regions along the no-flow edges of the fractures first and then gradually advances toward the center of the fracture with each increment of Ψ.After water filled the fracture along the two flow manifolds, the remaining air became entrapped and immobilized (dark red regions in second row from top in Figure 3).In contrast to the drainage process, the trapping observed during imbibition is similar to that observed in fractures bounded by non-porous, impermeable surfaces.As a result, potentially large regions of air may become entrapped regardless of the presence of the bounding porous matrix.
During secondary drainage (DR2) in FFH, air spreads more readily through the fracture than during DR1 due to the regions of trapped air remaining after IMB1.The result is that, for the horizontal experiments, similar distributions of air and water within the fracture occur at lower values of Ψ during DR2 than DR1.This can be observed in Figure 3 where the distributions of air-filled regions at Ψ = 398 mm in DR1 is similar to Ψ = 374 mm in DR2.Likewise the distribution of air-filled regions at Ψ = 423 mm in DR1 is similar to Ψ = 398 mm in DR2.Note, this history dependence is not observed for the imbibition cycles, where the initial distribution of air within the fracture was almost identical for IMB1 and IMB2, resulting in a nearly identical filling order (Figure 3).Similar behavior was observed in MFH, but because the drainage process occurred over a narrower range of Ψ, the differences between DR1 and DR2 are less pronounced.
The history dependence of the drainage cycles in the horizontal experiments is also evident in Figure 4.In addition to similar saturations occurring at lower Ψ during secondary drainage, the initial saturation at the end of primary imbibition differs from the saturation at the onset of secondary drainage.We believe this is because there was a 12-hr gap between the primary and secondary cycles; during that time some of the trapped air dissolved leading to an increase in S A w .Note, there was negligible hysteresis in the drainage curves for FFV.This is likely due to the stabilizing effect of gravity on the advancing interface during imbibition, which led to minimal trapping of air during imbibition (Figure 3).
In addition to the smaller Ψ a,e during primary drainage for the horizontal experiment in Model MF (Experiment MFH, Figure 4), another significant difference between FFH and MFH was the distribution of the air and water phases during each sequence.Specifically, the air clusters in MFH are more compact with less roughness of the air-water interfaces.Previous scaling analyses of two-phase displacements in fractures between non-porous matrices suggest a reasonable explanation for this behavior (Glass et al., 1998(Glass et al., , 2003)).The competition between interfacial curvature in the fracture plane and across the fracture aperture have been shown to control the geometry of the air-water interfaces; smaller fracture apertures with more aperture variability lead to more tortuous interfaces than larger aperture fractures with less aperture variability.Glass et al. (2003) derived the dimensionless parameter, C/δ = 〈b〉 2 σ b λ b , where C is the dimensionless curvature number, a ratio of in-plane to out-of-plane interfacial curvature, δ is the coefficient of variation of the fracture aperture, and 〈b〉, σ b , and λ b are the mean, standard deviation, and correlation length of the aperture field.They showed that small C/δ led to tortuous air-water interfaces and as C/δ became larger entrapped regions of air became more compact.Though we cannot directly quantify C/δ for our experiments, the air-entry apertures provide an estimate of 〈b〉.Because the aperture variability in our fractures are induced by the pore-scale variations of the porous surface, both σ b and λ b likely scale with the respective pore sizes of the porous plates.Because the maximum pore sizes for MF are approximately 4 times larger than those in FF, this suggests that C/δ is an order of magnitude larger for MF than for FF.This likely explains the difference in the structure of the air-water interfaces for the corresponding experimental sequences.
It is not possible to directly measure the relative permeability of the air and water phases within the fracture during our experiments, but we can estimate these values through numerical simulations in the measured airwater distributions.For these simulations, we considered only the influence of the geometry of the air and water estimated relative permeabilities (see Text S3 in Supporting Information S1 for details).Detwiler et al. (2005) showed that the role of aperture variability on estimates of fracture relative permeability were minor relative to the influence of the distribution of the flowing phases within the fracture.We used the local cubic law to simulate flow of both air and water through the fracture for each value of Ψ (and experimentally measured phase distribution) represented in Experiments MFH, FFH, and FFV.Figures 4c and 4d show the relationship between the estimated water and air relative permeabilities, k r,w and k r,a , respectively, and the areal saturation S A w of the water phase.
The relative permeability curves (Figures 4c and 4d) are qualitatively similar to what has been measured in both porous and fractured media in other studies, suggesting the potential utility of empirical permeability-saturation relationships for modeling flow through fractured porous media.However, the potential for the development of fracture-spanning regions of either air or water can strongly influence the evolution of k r,w .This is most notable in comparing k r,w during primary and secondary drainage for FFH and MFH.The large region of air that forms at the entrance to FFH (Figure 3) caused a significant decrease in k r,w once Ψ a,e was exceeded.Conversely, in MFH, the more compact shape of the evolving air region led to a more gradual decrease in k r,w .

Concluding Remarks
We have presented the development and evaluation of a new experimental system for exploring two-phase flow processes in porous fractured media.Use of light transmission through the translucent porous fracture surface allows us to quantitatively delineate the distribution of the evolving air-water interface within the fracture.Example experiments in two different fractures demonstrated the role of changing pore pressure in the porous matrix on the distribution of air and water within the fracture, which exhibited significant hysteresis from the primary drainage through subsequent drainage cycles.
The primary advantage of this method over other approaches (e.g., 2D micromodels with a channel bounded by a porous matrix or X-ray CT in rock cores) is the ability to resolve the distribution of air and water within a fractured porous medium at: (a) spatial scales that are much larger than the scale of aperture variability and the resulting regions of entrapped phases during displacement processes and (b) temporal scales with the potential to resolve potentially rapidly evolving interfacial dynamics.In addition, the demonstration experiments presented here used a smooth glass plate as the upper fracture surface, but such experiments can be readily extended to include a rough upper fracture surface to explore the relative importance of fracture-matrix interactions and two-phase flow processes within a bounding variable aperture fracture.

Figure 1 .
Figure 1.Schematic of experimental system including: (a) an overview of the experimental setup; (b) a plan view and crosssections of the fracture test cell; and (c) a schematic showing how, during drainage and imbibition, water exits and enters the fracture through the porous matrix and air enters and exits through the fracture edges, respectively.

Figure 2 .
Figure 2. Overview of the image processing procedure, including: (a) The raw grayscale reference image, I ref ; (b) an example of a raw grayscale experimental image, I exp , corresponding to the step Ψ = 398 mm during primary drainage (DR1) process in experiment FFH; (c) the corresponding light absorbance field, A abs = ln (I exp /I ref ); (d) the relationship between the average number of segmented air clusters, N ave , and the segmenting threshold values for Experiments FFH and FFV, T FFH and T FFV , revealing distinct minima for these curves; (e) the resulting binary image (black-water, white-air) determined using the global threshold with results using the upper and lower bounds of the global threshold indicated by yellow and red lines, respectively; (f) enlarged view of region indicated by the green box in (e).

Figure 3 .
Figure 3. Phase evolution during cyclic drainage and imbibition processes of Experiments MFH, FFH, and FFV.The average value of the applied capillary head, Ψ, is used to represent each pressure step.Discrete Ψ steps are indicated by the numbers next to the color bar for each sequence and the sequences proceed from top to bottom of each color bar.Warm colors are smaller values of Ψ and cool colors are larger values of Ψ.The color bar range for each drainage cycle begins at the first indication of air entry into the fracture.

Figure 4 .
Figure 4. Relationship between applied capillary head Ψ and areal water saturation S A w for Experiment MFH (a), FFH (c) and FFV (e); and relationship between modeled relative permeability of water or air (k r,w , k r,a ) and areal water saturation S A w for Experiment MFH (b), FFH (d) and FFV (f), in which S A w , k r,w and k r,a are the average values during the last 1 min in each step.