Experimental investigation of hysteretic dynamic effect in capillary pressure-saturation relationship for two-phase flow in porous media

Well defined laboratory experiments have been carried out to investigate hysteretic dynamic effect in capillary pressure-saturation relationships for two-phase flow in homogeneous and heterogeneous (layered) porous media. Conceptually, the dependence of the capillary pressure curves on the rate of change of saturation (dSw/dt) is defined as the dynamic effect in capillary pressure relationship, which is indicated by a dynamic coefficient, τ (Pa.s) and it determines the rate at which two-phase flow equilibrium is reached, i.e., dSw/dt = 0 where Sw and t are the water saturation and time, respectively. The dependences of τ on various fluid and porous materials properties have been studied in the context of drainage; but, there is limited study for imbibition and the hysteresis of τ-Sw relationships. As such, the emphasis in this paper is on reporting τ-Sw curves for imbibition while also demonstrating the hysteresis in the τ-Sw relationships by comparing τ-Sw curves for drainage (previously reported) and imbibition (this study) in carefully designed laboratory experiments. Homogeneous sand samples composed of either fine (small particle size, lower permeability) or coarse (larger particle size, higher permeability) sand have been used for these experiments. Furthermore, a layered domain made of a find sand layer sandwiched between two coarse sand layers is used as a model of heterogeneous domain. The results of the study confirm that the τ-Sw relationships are hysteretic in nature and, as such, the speed to flow equilibrium should vary depending on whether drainage or imbibition takes place. At a particular water saturation, the magnitudes of the dynamic coefficient (τ) are found to be generally higher for imbibition, which imply that the speed to flow equilibrium at the same saturation will be slower for imbibition.


Introduction
Mathematical description of two phase flow (immiscible fluids) in porous media requires appropriate governing equations for the conservation of fluids mass and momentum as well as other constitutive equations, e.g., capillary pressure (P c )-saturation (S w ) relationship. 1-3 Traditionally, an extended version of Darcy's law is used as the governing equation of motion for the fluid phases. The conservation of fluid mass is given by an equation for the conservation of phase saturation, i.e., the ratio of the volume of a fluid phase in a given domain to the pore volume in the domain. The relationship between P c and S w is described by various empirical models such as the  and van Genuchten 5 formulations. For these relationships, capillary pressure is generally calculated by an empirical relationship obtained under equilibrium condition between the individual fluid phase pressures as follows 1,6,7 : c w nw P P P = − (1) where, P nw and P w are the average pressures of the non-wetting and wetting phases, respectively. In equation (1), P c is defined to be a function of S w in the porous medium. Experimentally, the P c -S w curves can be obtained by taking a porous medium initially saturated with a wetting fluid (e.g., water) and then letting it to gradually drain off by increasing the capillary pressure and displacing the wetting fluid (e.g., water) by a non-wetting fluid (e.g., oil). The resulting curve representing the corresponding values of wetting phase saturation and capillary pressure at equilibrium condition (i.e., dS w /dt = 0) is known as the primary drainage P c -S w curve. In most cases, higher capillary pressure does not lead to any further displacement of the wetting phase due to phase separation and/or strong wetting phase attachment to the grain particles. The saturation at which this condition occurs is known as the irreducible wetting phase saturation, S wir . 6,8 Once the stage of irreducible wetting phase saturation is reached the saturation of the wetting phase in the sample may be increased by displacing the nonwetting phase by decreasing capillary pressure. A plot of these experimental results provides a main imbibition P c -S w curve. The drainage and imbibition P c -S w curves generally do not coincide because the pore spaces of the porous medium imbibe and drain differently. 1,6-9 Further, an imbibition experiment does not lead to full wetting phase saturation (S w = 1) at zero capillary pressure because some of the non-wetting phase will be entrapped as isolated bubbles in the pores and thereby not displaced. This saturation value is referred to as residual non-wetting phase saturation (S wr ). The primary drainage and main imbibition curves start from full saturation with wetting fluid and from irreducible wetting phase saturation, respectively. The difference in the P c -S w curves is referred to as hysteresis. The hysteresis is a result of a number interdependent pore scale mechanisms during drainage and imbibition (e.g., dynamic contact angle) and has been the subject of many studies in the literature. [9][10][11][12] Further, the hysteresis is not just a matter of the primary drainage and main imbibition curves as the P c -S w relationships may follow infinite number of scanning curves depending on when the drainage or imbibition process is reversed. However, these two curves can be considered as the limiting curves, and as such, they are useful in determining the extent of hysteresis P c -S w relationships. The present study will limit itself to these curves and will not consider primary, secondary or any other scanning P c -S w curves.
The hysteresis in equilibrium and dynamic P c -S w relationships and, in particular, the ways to measure/quantify the hysteresis are matters of great research interests in the literature. [13][14][15][16][17][18][19][20][21][22][23] There have been many discussions on the validity of Darcy's law and equation (1) for two-phase flow processes in the last two-decades. 7,20,24 It has been argued that at shorter time duration when the twophase flow may not necessarily be at equilibrium the traditional approaches as described above cannot be applied. In fact, the early experiments of Topp et al. 14 , Smiles et al. 15 , Vachaud et al. 16 and Stauffer 17 have indicated that P c -S w curves at equilibrium and dynamic conditions are different. Further, it has been shown that at short time durations the P c -S w relationships depend on the rate of change of fluid saturation (dS w /dt). The dependence of the P c -S w curves on the rate of change of saturation (dS w /dt) is known as the dynamic effect in capillary pressure relationship. 7,25 The dynamic effect has been discussed significantly in the last two decades and it has been shown to be of importance, as discussed below.
Hassanizadeh and Gray 24 presented the following relationship, which included the dynamic effect indicated by a dynamic coefficient in the conventional P c -S w relationship: where, c dyn P is the dynamic capillary pressure defined as w nw c dyn P P P − = and c stat P is the capillary pressure at equilibrium conditions (dS w /dt = 0). The coefficient τ is called the dynamic coefficient which is a measure of the rate of change in saturation and hence the speed to equilibrium condition. In the last decade a number of studies have discussed the significance of dynamic capillary pressure and τ in different circumstances. 7, 19-21, 23, 25-43 For example, hysteretic dynamic effect in the capillary pressure relationship was discussed in a theoretical study by Beliaev and Hasanizadeh 19 . Hassanizadeh et al. 20 interpreted the imbibition experiments (displacement of oil by water) of Kalaydjian 44 and reported that the average value of the dynamic coefficient in Kalaydjian's experiments was 2×10 6 kg/m.s. Hassanizadeh et al. 21 carried out a laboratory study to determine τ values for primary drainage (PD), main drainage (MD) and main imbibition (MI) for PCE and water flow in a homogeneous porous medium (permeability =10 -11 m 2 ; porosity = 0.35-0.4). In this study, τ values were reported to be in the range 5×10 4 -6×10 4 Pa.s for primary drainage, 1.3×10 5 -1.9×10 5 Pa.s for main drainage and 3.1×10 5 -7.7×10 5 Pa.s for main imbibition, which imply that the dynamic coefficient was found to be higher for imbibition than that for drainage by Hassanizadeh et al. 21 Another study which is relevant to the current paper was conducted by Sakaki et al. 23 Sakaki et al. conducted PD, MD and MI experiments and measured τ values for these three cases for a single porous domain of saturated hydraulic conductivity 0.016 cm/s. They calculated a characteristic redistribution time (τ B ) of fluids defined originally by  carried out a comparison of τ B for PD and MI which showed differences in τ B for the two cases suggesting hysteretic τ B and hence, hysteretic τ. We are not aware of any other paper, which has conducted experiments for measuring dynamic coefficient for imbibition and/or hysteretic behaviour of the dynamic coefficient. There is however some studies which attempted to simulate dynamic imbibition of porous medium using hypothetical values of τ. For example, Manthey et al. 48 used two sets of hypothetical τ values of 0, 10 5 , 10 7 Pa.s and 0, 10 7 Pa.s to simulate imbibition in homogeneous porous domains with permeabilities of 9.4×10 -10 m 2 and 1×10 -12 m 2 , respectively.
In order to have a better understanding of the hysteretic behaviour of τ-S w relationships, it seems we need more experiments, which can investigate the extent of this hysteresis. Given that both τ and dS w /dt change with S w and time, a direct comparison of τ B for PD and MI as done by Sakaki et al. (2010) may not reflect the true extent of hysteretic behavior of τ. A number of recent studies 32, 41,42 on experimental determination of τ values for drainage have indicated that τ depends on the medium properties (e.g., particle size, permeability, heterogeneity) and size of the domain. However, it is not clear how do τ for imbibition and the hysteretic behaviour of τ depend on the medium properties (e.g., particle size, permeability, heterogeneity). Note that the experiments of Hassanizadeh et al. 21 and Sakaki et al. 23 were done for homogeneous domain of a single permeability. Further, the results of Sakaki et al. 23 were based on the measurements at the middle height of a 10 cm column and it is not clear to us how do these local measurements may determine the effective τ-S w relationships at the scale of the experimental domain (core scale). This issue may become particular important if the domain is heterogeneous in nature. Another important issue that should be revisited in the context of both drainage and imbibition is the trend/shape of τ-S w curves. Many studies which reported τ values for drainage suggest that τ increases as S w decreases. 7,20,23,26,27,41,42 However some recent studies on the measurements of τ for drainage have found different trends where τ values decrease, fluctuate or remain almost constant as S w decreases, see, e.g.,  This paper aims to address the above points by using carefully controlled laboratory experiments of two-phase flow in homogeneous and layered porous media. In this work, dynamic and quasi-static capillary pressure-saturation (P c -S w ) relationships were used to obtain the dynamic coefficient (τ) as a function of water saturation (S w ) for primary drainage (PD) and main imbibition (MI). Homogeneous porous samples composed of either fine grained (low permeability) or coarse grained (high permeability) sand have been used for the experiments. Furthermore, the same coarse and fine sand particles were used to create a heterogeneous (layered) domain composed of a fine sand layer sandwiched between two coarse sand layers. The coarse and fine sand layers in the heterogeneous domain are created so that the material properties (particle size distribution, porosity, intrinsic permeability) are the same as in the homogeneous coarse and fine sand domains. The τ-S w curves are then compared for PD and MI to determine the hysteresis in the homogeneous and heterogeneous sand samples. We determine τ-S w relationships at three different heights within both the homogeneous and heterogeneous domains and, use the local τ-S w data to calculate the effective τ-S w curve for these domains. Please note that the homogeneous and heterogeneous domains used for the purpose of this paper are the same as presented earlier by Das and Mirzaei 41,42 who reported τ-S w curves for drainage for homogeneous and layered domains. Also, as discussed by Das and Mirzaei 41,42 , the layered domain in our work should be treated as a 'weakly layered' domain as the contrast in permeability between the layers is not significant.
For the purpose of this paper, additional experiments have been conducted to obtain τ-S w curves for imbibition for the same domains as used by Das and Mirzaei. 41,42 The τ-S w curves for imbibition are then compared with the curves for drainage reported earlier by Das and Mirzaei 41,42 to determine the significance of hysteresis in τ-S w curves. In effect, the results presented in this paper not only provide continuity of the results presented by us earlier but also allow easy comparison of τ-S w curves for drainage and imbibition for the same conditions.
Please also note that we are aware that the dynamic capillary pressure effect is often interpreted according to a characteristic redistribution time (τ B ) using a model by  However this paper is concerned with the dynamic effect as presented by the dynamic capillary pressure relationship (equation 2) proposed by Hassanizadeh and Gray 24 and the behaviour of the dynamic coefficient (τ). As such, we do not discuss the work related to the Barenblatt model in this paper.

Materials and Methods
Detailed descriptions of the materials and methods for drainage experiments were described by Das and Mirzaei. 41,42 The experiments that were designed to obtain an understanding of the hysteresis in τ-S w relationships simply included additional experiments for dynamic and quasi-static imbibition of two-phase flow. In consistent with our previous experiments 41,42 , we performed the flow experiments ( Figure 1) in a cylindrical cell packed with the chosen porous sample. In each cycle of experiment, silicone oil was injected in the sample through a hydrophobic filter and water drained out of the cell though a hydrophilic filter. Three mini-time domain reflectometer (TDR) probes were installed at three different heights in the sample to determine water content during the experiments. Pressure transducers equipped with either hydrophilic or hydrophobic filter mounted on the wall of the cell at the same height of mini-TDR probes were used to selectively monitor water and oil pressure, respectively. Since the descriptions of the methods and materials have been given earlier by Das and Mirzaei 41,42 , we provide only brief discussions on materials and methods in this paper for its completeness.

Porous Media and Fluid Properties
The laboratory experiments were performed using Leighton Buzzard DA 14/25 sand as coarsegrained and Leighton Buzzard DA 30 as fine-grained sand (WBB Minerals, Cheshire, UK). They provide good uniformity of particles (uniformity coefficient (d 60 /d 10 ) of 1.21 and 1.32 for fine and coarse grained particles), high sphericity, high chemical purity and very low organic matter content. The nonwetting phase liquid selected for this study was silicone oil (viscosity: 200 cSt) with a negligible solubility in water, a low volatility at room temperature and a low health hazard. The aqueous phase in our experiments was distilled water. Relevant properties of the porous samples and fluids are listed in Table 1.

Sample Preparation
The flow experiments were carried out in a cylindrical cell of 10 cm internal diameter (ID) and 12 cm height filled with a porous sample (homogeneous or heterogeneous) made by sand particles. The sample was first purged of air by stirring sand in distilled water reservoir until no air bubble came out of mixture. Subsequently, the mixture was put in a vacuum unit to release any trapped air bubble for 24 hours. The cell was carefully mounted on a plate equipped with an O-ring to prevent any air inflow and a deaired and fully water saturated hydrophilic filter to facilitate water drainage and prevent oil flow at outflow. Subsequently, the column was filled with distilled water. Wet and de-aired sand poured into water whilst column was on top of a mechanical shaker to provide even and dense fully water saturated packing. Visual observation confirmed no gas bubbles in the sand pack following this procedure.
While preparing the homogeneous porous samples, either coarse grained or fine grained sand was deposited in the cell. To prepare the heterogeneous (layered) domain, a predefined mass of de-aired and wet coarse sand was put in the cell to cover 3 mm above the radius of influence of lowest mini-TDR probe. The radius of influence is the distance measured from outer probe rods, within which TDR reflection is not interfered with the reflection of other TDR probes and/or reflection of probe is not affected by the texture of sand out of radius of influence. 41 It is defined by the manufacturer to be 3 mm for the type of the mini-TDR probes used. Subsequently, a known quantity of fine sand is added to the cell to make a fine sand layer of 39 mm thickness to cover the middle mini-TDR probe. After this, another layer of coarse sand is deposited on the top of the find sand.
An X-ray test to measure packing density along the cell was also conducted before any two phase flow experiment to ensure homogeneity in packing within either the whole domain or a layer in a heterogeneous domain.
For homogeneous domain, the porosity of the sand pack was calculated from the total mass of each sand pack. For heterogonous domain, the porosities of individual coarse or fine sand layers are not measured separately and are defined to be the same as the corresponding homogeneous domains of coarse or fine sand domains. This is because the mass/volume ratios of the homogeneous coarse and fine sand samples were the same way as those of individual coarse and fine sand layers, respectively.
Before conducting two-phase flow experiment, single phase flow using water and closely following constant head permeability test was performed to measure intrinsic permeability of the sand pack. To perform this experiment, a hydrophilic filter was placed on top of the sand column and the inflow reservoir was filled with distilled water.

Quasi-static and dynamic two-phase flow experiments
In the designed cell (Figure 1), three mini-TDR probes are installed at three different heights in the domain to measure water content during the experiment. Pressure transducers equipped with either hydrophilic or hydrophobic filters are mounted on the wall of the cell at the same height as mini-TDR probes which are used to monitor water and oil pressure, respectively. 41,42 Locally measured water saturation and, water and oil pressures are used to calculate local saturation and capillary pressure, respectively, which are then used to construct either dynamic or quasi-static P c -S w curves. These curves are subsequently used to calculate dynamic coefficient.
The quasi-static two-phase flow experiments during drainage were conducted as follows. The outflow valve of the cell was opened to provide very low flow rate during experiment and it was levelled with the top of the sand to overcome hydrostatic head pressure gradient to minimize gravity effect 41 . The initial pressure of silicone oil was zero within sand column. Then, the pressure of silicone oil on top boundary was gradually increased by increasing air pressure at the top of an oil reservoir, which was the Marriotte bottle. This resulted in infiltration of silicone oil through hydrophobic filter at the top of cell. Injected oil also displaced water out of sand. The experiment was carried out until steady-state flow condition was reached; that was, a constant rate water flows from the outflow valve. Water saturation was calculated considering initial water content and outflow water volume and compared with the mini-TDR readings. Silicone oil and water reservoir and pressure transducers pressure facilitated P c measurement for comparison. Calculated S w and P c from the reservoirs and recorded sensors provide one point of the equilibrium P c -S w curve. The experiment was continued until a new steady-state was reached for higher fluids pressures. A second point for P c -S w curve was thus obtained. This procedure was repeated until displacing oil front reached bottom hydrophilic filter; that was, measured water content at lower TDR probe reached irreducible water content. At this point, the sample was assumed to have reached its irreducible wetting phase saturation (S wir ).
The procedures adopted to determine the dynamic P c -S w curves for drainage are as follows. Silicone oil pressure at the top of domain was increased to a high pressure at once and maintained constant by constant air pressure on Marriotte bottle provided through a manual pressure regulator. Table 2 shows the applied boundary condition for different dynamic flow experiments. The drainage displacement was continued until the saturation at the lower mini-TDR probe reached its S iw . Mini-TDR measured water content and pressure transducer recorded oil and water pressure provide data to calculate three local dynamic P c -S w curves for each boundary pressure. Dynamic experiment was repeated for four different pressures.
Following the drainage experiment, imbibition experiments are performed in the same domain, which is at its irreducible water saturation. For this, the direction of fluid displacement is reversed. Water is injected from bottom of the domain through the hydrophilic filter and oil is ejected through the hydrophobic filter at the top of the cell. Imbibition experiment ends when water front reached hydrophobic filter at the top boundary. In this condition, the domain has reached its residual saturation (S wr ). The same approach was adopted for obtaining both the quasi-static and dynamic P c -S w curves.

Calculation of dynamic coefficient
Equation (2) shows that if P c,dyn -P c,equ and d w S /dt are known at a S w value, τ can be determined. From the experimental data for P c,dyn , P c,equ , S w and d w S /dt at different S w , the local τ-S w relationships are determined at different heights in the porous sample using equation (2). Furthermore, the average τ and S w values are calculated for the whole domain by using the following equations (Das and Mirzaei, 2012a,b): Where w S is the volume weighted water saturation, which is also defined as the average/effective water saturation for the whole domain.
is the measured saturation at measurement volume V i at corresponding measurement height i, and i=1,2,3,…,n is the number of measurements heights in which the τw S curves are calculated, i.e., the number of τw S curves, in which n=3.

Transient saturation (S w -t) profiles
In this section, we present the S w -t curves obtained in this study. Although the trend of these curves are well understood, it is important to discuss them in the present context as they help to understand when the flow system reaches equilibrium (∂S w /∂t =0). These curves are also needed to calculate ∂S w /∂t and the dynamic coefficient using Equation (1). Figure 2 shows some typical S w -t curves at three measurement heights during primary drainage and main imbibition in coarse sand domain. The drainage S w -t curves in the domain have been discussed by Das and Mirzaei (2012) and as such, these curves are only included in this paper to demonstrate the differences between the S w -t curves for drainage and imbibition. As expected, the water saturation in the domain decreases during the drainage and increases during the imbibition cycle. Figures 2(a) and 2(c) use the time elapsed from the start of the experiment while Figure 2(b) uses the time from the moment the first change of saturation is recorded. Figure 2(a) demonstrates that the water saturation in the domain (i.e., top, middle and bottom of the domain) at the end of an imbibition cycle for a pressure of 8 kPa is almost the same everywhere. However, the superimposed S w -t curves in Figure 2(b) for the same boundary condition show that although the saturation is almost the same throughout the domain, saturation profiles at different measurement heights may vary. Consequently, contrary to the drainage displacement, the S w -t curves do not completely overlap during imbibition. Figure 2(c) displays the S wt curves for imbibition displacement in the same domain for different boundary conditions of 8, 9, 10 and 11 kPa. Figures 3(a) shows typical drainage and imbibition S w -t curves in fine sand domain for 8 kPa boundary pressure. At first glance, the S w -t curves at different measurement heights seem to follow a consistent trend. However, superimposed S w -t curves in Figure 3(b) show that although the curves follow a consistent trend and end up at almost similar saturation value, towards the end of an experiment, the shape of the S w -t curve at upper the measurement point varies from the other curves slightly. A close comparison between Figures 2(b) and 3(b) shows that contrary to the overall behaviour of all S w -t curves in coarse sand domain, the curves in fine sand domain overlie for a large range of saturation. The S w -t curves at different measurement heights in fine sand for 9, 10 and 11 kPa boundary pressure also follow similar trends. Again, for a large range of saturation, the curves overlap. Figure 3(c) shows the S w -t curves at three different measurement heights for various boundary conditions. The figure clearly shows that saturation distribution along the homogeneous fine sand domain in imbibition is fairly uniform, which is consistent with the observations made for homogeneous coarse sand. The fact that the S w -t curves for both drainage and imbibition at different heights within the coarse or fine sand domains in this work are not significantly different is not completely surprising given that the domain size is small and the distances between two measurement points are small.
As in the homogeneous domains, imbibition experiments were carried out in the layered domain ( Figure 4) for oil pressures of 8, 9, 10 and 11 kPa. Figure 4(a) displays the drainage and imbibition transient saturation profiles in the domain for 8 kPa boundary pressure in real time scale. However, the superimposed curves in Figure 4(b) make the comparison of S w -t curves easier. The Figure shows that the shapes of S w -t curves at different measurement heights are similar at first, which follow a consistent increasing trend. As the experiments progress and water saturations increase, the S w -t curves at different measurement heights show slightly different behaviour. This is obviously due to the variations in the properties of the sands in different layers, which affect the fluid flow behaviour. The same trends in S w -t curves for 9, 10 and 11 kPa boundary pressure are found. S w -t curves for imbibition for all pressures are shown in Figure 4(c).

Dependence of the rate change of saturation (dS w /dt) on time (t)
The time derivative of saturation (dS w /dt) is the gradient at a point on a S w -t curve. As explained before (Das and Mirzaei, 2012a,b), one needs the values of dS w /dt at different saturation and time to determine a dynamic coefficient values. The importance of the dS w /dt-t curves and their implications are also discussed earlier (e.g., Hassanizadeh et al., 2002;Das and Mirzaei, 2012a,b) and are not discussed in detail in this paper. In general, the rate of change of saturation (dS w /dt) changes with saturation and time and this determines when a two-phase flow equilibrium is reached (dS w /dt = 0). This section reiterates this point and reports some typical dS w /dt-t curves for main imbibition and compare them with the curves for primary drainage with a view to demonstrate the hysteresis in these curves which would be related to the hysteresis in the τ-S w relationship (discussed later). The values of dS w /dt are approximated as shown below, S + are the average wetting phase saturation at t n-1 , t n and t n+1 calculated using equation 3 in different heights. It should be noted that the dS w /dt values are also related to the mobility ratios, which depend on the relative permeabilities of the fluid phases and the fluid viscosity ratios . In other words, if the mobility ratios change then dS w /dt and, hence, the dynamic coefficient may also change. In the present paper, no attempt is made to address this point as it requires experiments with different fluid pairs.
Figures 5(a) and 5(b) show some typical imbibition and drainage dS w /dt-t curves in the coarse sand domain for 8 kPa boundary pressure. As expected, the drainage dS w /dt have negative values (Das and Mirzaei, 2012a) while the values of dS w /dt during imbibition are positive. This is logical as S w decreases with time in drainage and increases during imbibition. As the flow reaches equilibrium dS w /dt values for both drainage and imbibition reach almost zero. These points are further illustrated in Figure 5(b), which shows the superimposed imbibition and drainage dS w /dt-t curves. The magnitudes of dS w /dt vary for drainage and imbibition ( Figure 5(b)) at different heights within the domain (Figure 5(c)). The differences in the dS w /dt values at the same time and saturation give rise to different dynamic effects, which is discussed later. Figure 5(c) shows the dS w /dt-t curves for different pressures of 8, 9, 10 and 11 kPa. Within the pressure range tried in this work, dS w /dt values do not seem to be significantly different. The observations related to the dS w /dt-t curves for coarse sand domain seem to be valid for fine sand too ( Figure 6) and the layered domain (Figure 7). Dynamic and quasi static capillary pressure (P c ) -saturation (S w ) curves Figures 8(a) and 8(b) present the hysteresis in the dynamic drainage and imbibition P c -S w curves at 8 and 11 kPa pressures. Similar patterns are also observed for oil pressures of 9 and 10 kPa. Figure  8(c) shows the hysteresis in quasi-static P c -S w curves for drainage and imbibition. Figures 8(a) and 8(b) show the curves at the upper, middle and lower measurement heights of coarse sand domain. As demonstrated in the figures, the drainage P c -S w curves lie above the imbibition P c -S w curves, demonstrating hysteresis in the drainage and imbibition curves. In general, the hysteresis in the P c -S w curves can arise from the fact that capillary pressure depends on the adhesive forces between the fluids and the minerals which make up the pore walls, the pressure difference between oil and water phase on the pore spaces and, hence the curvature of meniscus between oil and water in the pore space, etc. Although the porous medium and fluids are the same in drainage and imbibition, it is difficult to pinpoint exactly which of these factors is the least or most significant in causing the hysteresis in this study. However, it is clear that the difference of water and oil pressure and, consequently the saturation distribution in filling the pore spaces during drainage and imbibition result in different capillary pressure behaviour resulting in different P c -S w curves during drainage and imbibition. Similar interpretation is made from Figure 8(c), which demonstrates that the quasi static P c -S w curves are also hysteretic in nature.
The imbibition experiments following the drainage experiments in coarse sand domain has been conducted for fine sand (Figure 9), which also show hysteresis in P c -S w curves. Figure 10 shows some typical graphs to compare the dynamic and quasi static drainage and imbibition P c -S w curves at the three measurement heights in layered media. As in the case for coarse sand domain, hysteresis in imbibition and drainage P c -S w curves in heterogeneous domain is present in these figures.
The hystereses in P c -S w curves in Figures 8 -10 mean that the average saturation in the domain is not the same for drainage and imbibition at the same capillary pressure, which is due to different pore filling and emptying rates. There have been a number of pore-scale studies in the literature, which may explain the hysteresis in terms of various pore scale physics (e.g., Tsakiroglou and Payatakes 1998;Tsakiroglou and Avraam, 2002;Liu, et al., 2011Liu, et al., , 2012. However in the context of this study, we measure averaged P c -S w curves at the core scale which are then used to determine the hysteretic dynamic coefficient (τ) -saturation (S w ) relationships. Difference in average S w along with different dS w /dt at the same P c or/and time level give rise to hysteretic dynamic effect in capillary pressure relationships as discuss in the next section.

Hysteretic dynamic coefficient (τ) -saturation (S w ) relationship
In order to understand the hysteretic behaviour of the dynamic coefficient, the τ-S w data for homogeneous and layered domains for primary drainage and main imbibition are presented in Figures  11 and 12. We also present a table (Table 2) which summarises the important data ranges and statistics of our experimental results for the dynamic coefficient (τ). The table would be invaluable in following the discussions in this section. Figure 11(a) displays the drainage and imbibition τ-S w curves at three different heights in the coarse sand domain. The effective τ-S w curves are also presented in the figure, which are determined using equations (3) and (4) as mentioned earlier. The hystereses in the local and the effective curve τ-S w curves are found, which imply that the energy required to achieve flow equilibrium (dS w /dt = 0) at the same water saturation is different for drainage and imbibition. As discussed earlier, the dS w /dt values vary slightly for drainage and imbibition (Figures 5-7) at the same saturation and time which cause the hysteresis in the τ-S w curves. Analysing the pore scale processes is beyond the scope of this paper. However, it can be safely stated that these differences in dS w /dt values and, hence, the τ-S w curves, are due to the difference in the pore-scale flow processes during drainage and imbibition, e.g., difference in pore filling and emptying rates, and mechanisms (e.g., percolation-type piston drive in drainage, snap-off and cluster growth in imbibition), wettability, etc. What we observe in the τ-S w curves is the lumped affect of all possible pore scale processes. In consistent with the conclusion of Hassanizadeh et al. (2005), the imbibition τ values in Figure 11(a) are found to be bigger than the drainage τ values at the same S w , meaning that the speed to flow equilibrium during imbibition is slower than the drainage at that S w . The imbibition τ-S w curves for coarse sand show that the dynamic coefficient varies non-linearly with the water saturation similar to the drainage τ-S w curves (Das and Mirzaei, 2012a,b). But, as the water saturation increases from irreducible water saturation, the imbibition τ-S w curves break away from that for drainage and they remain almost constant. Similar trends are observed for fine sand domain. Table 2 shows that our imbibition τ values for both coarse and fine grained sands are consistent with the results of Hassanizadeh et al (2005) who reported τ values in the range 3.1×10 5 -7.7×10 5 Pa.s for main imbibition. In the imbibition experiments of Sakaki et al (2010), the dynamic coefficient at the middle of their domain was found to be relatively constant in the S w range of 0.3 -0.8 with the τ values ranging between 1.0×10 6 and 6.0×10 6 Pa.s. Our local and effective τ values for imbibition are smaller than those obtained by Sakaki et al. (2010). However, the slopes of our imbibition τ-S w curves for both coarse and fine sands are small in the range 0.5<S w <0.8 which imply that we also observe almost horizontal τ-S w curves in consistent with the results of Sakaki et al. (2010). This practically means that the energy required for the two-phase flow system to reach equilibrium is almost the same in the middle saturation points during imbibition.
In Figure 11(b), we demonstrate the local and effective imbibition τ-S w curves for the coarse and fine sand domains. In effect, the figure shows if the effective τ-S w curves for homogeneous porous media at the core scale are any different from the local τ-S w relationships while also demonstrating any effect of the particle size and permeability on the imbibition τ-S w curves. The results show that there are not significant differences among the local and effective τ values at different saturations for most practical purposes. This observation is somewhat expected as the domain size is not very big (12 cm long) and it is homogeneous in its properties (e.g., particle size and shape). Further, the sample volumes over which each set of pressure transducers and TDRs measures the τ-S w curves are small which are also equal. Similar qualitative conclusion has been made by others (Camps-Roach et al., 2010;Das and Mirzaei, 2012a) for drainage in homogeneous domains who have found that there are not significant differences among τ-S w curves at different levels of a core scale domain. Whether or not one obtains significantly different τ values at different points in a larger homogeneous domain or for different boundary conditions is not certain and it should be explored separately. Bottero et al (2011) have calculated effective τ for primary drainage in a bigger homogeneous domain of 21 cm; however, it is not clear to us if τ values at different points within the domain are different and, if so how significantly. Figure 11(b) suggests that our τ values obtained for imbibition are higher in fine sand than the values at the corresponding heights in coarse sand at a given saturation. As such, the average τ values in Table 2 are also higher for fine sand than those for coarse sand. These values are smaller than those of Kalaydjian (1992) but similar in magnitude to those obtained by Hassanizadeh et al (2005) and Sakaki et al (2010).
To make a further analysis of the drainage and imbibition τ-S w curves in fine and coarse grained porous domains, plots of the effective and local τ-S w relationships measured at different domain heights are shown in Figures 11(c) and (d). Figure 11(c) shows the drainage and imbibition τ-S w curves at the middle measurement heights of both coarse and fine sands. On other hand, Figure  11(d) shows the effective drainage and imbibition τ-S w relationship in the two sand domains. The figure show that both the drainage and imbibition τ-S w curves follow non-linear functional dependence and there is clear hysteresis where the imbibition τ values are generally larger as compared to the τ values for drainage at the same saturation. It suggests that for the imposed boundary conditions the two flow system would require more energy to reach equilibrium during imbibition, as such, one should expect hysteresis in the τ-S w curves.
As discussed earlier, the τ values depend on a number of other factors, e.g., permeability of the domain, fluid properties, domain size, temperature, etc. Therefore, the results of this study cannot be compared directly with values from other study, which may use different fluid, porous medium or domain size. Nevertheless, comparisons of our results with those obtained by others reveal interesting facts. Hassanizadeh et al. (2005) used a domain of permeability 10 -11 m 2 while our fine grained domain has a permeability of 3.1×10 -10 m 2 . Hassanizadeh et al. (2005) used a fluid pair of PCE and water while we used silicone oil and water as fluids with properties as shown in Table 1. Although there are significant differences in media and fluids properties, it seems that our results are similar to those of Hassanizadeh et al. (2005). For example, both studies obtain an imbibition τ value of approximately 3×10 5 Pa.s at 70% water saturation for the cases mentioned above. This implies that the speed to flow equilibrium during imbibition in these two systems at the conditions above would be almost the same despite difference in the properties of the flow system. Figure 11 that the τ values during imbibition are higher for fine grained sand (i.e., the domain with smaller particles and less permeability) than coarse grained sand is consistent with previous studies which measured τ for drainage Camps-Roach et al. 2010;Das and Mirzaei, 2012a). The uniformity coefficients, which indicate the uniformity of the particle shape, were 1.6 and 2.3 for the particles that Camps-Roach et al. (2010) have used. The particles that we have used have uniformity coefficients of 1.21 and 1.32 for fine and coarse particles, respectively, and these have better uniformity in shape than the particles of Camps-Roach (2010) Das and Mirzaei (2012a,b) and the current paper imply that at a particular saturation the τ value is expected to be higher for fine sand domain (i.e., domain with smaller particles and less permeability) for both drainage and imbibition for a range of particles shape and sizes.

The observation in
Although most results in this paper are consistent with other results in the literature, the shapes of the τ-S w curves for drainage and imbibition (and hence their interpretation) are not necessarily consistent with all studies. Sakaki et al. (2010) report that the τ values from their main imbibition experiments are slightly smaller than those in the primary drainage. This is opposite to the observation of this paper (Figure 11) as well as of Hassanizadeh et al. (2005).
The results for both PD and MI in our study suggest that the τ values are inversely proportional to the water saturations. Although the values of the dynamic coefficient may vary in different studies and similar trends in τ-S w curves are found by most Dahle et al., 2005;Mirzaei and Das, 2007;Das et al., 2007;Fučík et al., 2010;Sakaki et al., 2010;Das and Mirzaei, 2012a,b), some papers have been published recently where different trends of τ-S w curves have been observed (Camps-Roach et al., 2010;Sakaki et al., 2010). During drainage, Camps-Roach et al. (2010) find that the τ values decrease as water saturation decreases. Sakaki et al. (2010) report that their τ values increase as water saturation decreases during drainage. However, the imbibition τ values of  fluctuate, decreasing and then increasing with decrease in water saturation. Bottero et al (2011) reported τ values for PD for a saturation range of 0.55 -0.85 and found that τ was almost constant in that saturation range. It is not exactly clear to us at this moment why there are distinctly different trends of τ-S w curves. We are likely to revisit this question in another paper in the future. At this moment, it is suffice to say that the τ-S w curves of this paper are consistent with the majority of results found in the literature. Figure 12(a) displays the imbibition and drainage τ-S w curves at different measurement heights in the heterogeneous (layered) domain. As shown in the figure, imbibition τ-S w curves lie higher than the drainage τ-S w curves for the domain, which is consistent with the results obtained for homogeneous domain (Figure 11). Contrary to the drainage τ values, which increase as saturation decreases, the imbibition τ values decrease sharply initially as saturation increases from irreducible water saturation but as the imbibition proceed and saturations increase further, the dynamic coefficients decline slowly and they become almost constant. In Figure 12(a), the effective τ-S w curves for both drainage and imbibition are also shown. For both drainage and imbibition, the effective curves fall between the curves for coarse sand (upper and lower measurement point) and fine sand (middle measurement point) for the heterogeneous domain used in this work. The results suggest that the effective τ-S w curves are dominated by the τ values of the coarse sand as it occupies the maximum volume of the domain. The results also show that the curves are hysteretic in nature (similar to the local the τ-S w curves).
In Figures 12(b) -(d) we attempt to compare the τ values during imbibition at different measurement levels of the homogeneous and layered domains. As the properties of the coarse and fine sand in these experiments as well as the sample size, boundary condition, fluid properties are the same; this comparison may reveal interesting information. In Figures 12(b) and 12(d), we present the imbibition τ-S w curves at the upper and lower measurement heights of the homogenous coarse sand and layered domain, which are occupied by coarse sand. As demonstrated in these figures, although the imbibition τ-S w curves for the same sand type but different samples show similar trends, the imbibition τ-S w curve of coarse sand in heterogeneous domain lie slightly higher than the imbibition τ-S w curve of coarse sand in homogeneous domain. This means that the fluid displacement in the upper and lower layers (coarse sand) in heterogeneous domain needs more energy to reach equilibrium, which is logical considering the presence of fine sand layer at the middle of the heterogeneous domain. On the contrary, the imbibition τ-S w curves at the middle measurement point of the homogeneous fine sand and heterogeneous domain show slightly different behaviour (Figure 12c). When saturation is low at the beginning of the experiment, imbibition τ values in the fine sand layer at the middle of heterogeneous domain seem to be higher. However, as saturation increases, the two curves almost overlie and at some saturation values imbibition τ for fine sand layer in heterogeneous domain is less than the imbibition τ data at corresponding height in the homogeneous fine sand domain. Again, this is due to higher permeability of upper and lower coarse sand layers in the heterogeneous domain, which makes faster displacement than that in homogeneous fine sand domain filled with lower permeability sand. In general, however, we believe, there is not a significant difference between the τ values at the different heights in the homogeneous and heterogeneous domain in this study due to small size of the domain and weak permeability contrast between coarse and fine sand layers.

Conclusions
Most previous studies, which attempted to relate the dynamic effects in capillary pressure relationships to process parameters (e.g., grain size, permeability, temperature, domain size, saturation, fluid properties, boundary conditions) find that τ values depend on these parameters. As such it is expected that the τ-S w curves for both drainage and imbibition cycles, and hence the hysteresis in the τ-S w curves, will be affected by these parameters. However, it seemed there were not enough experiments to provide a detailed understanding of these effects. This paper has attempted to address this point by carrying out a number of laboratory experiments. Dynamic and quasi-static capillary pressure-saturation (P c -S w ) curves for primary drainage and main imbibition were determined in this work which were then used to obtain the dynamic coefficient (τ) as function of water saturation (S w ). Both homogeneous and heterogeneous (layered) domains were used in this work. The results in this paper clearly show that the energy required for, or the speed to, flow equilibrium at the same water saturation varies depending on whether a drainage or imbibition takes place, which confirms that the dynamic effect in capillary pressure-saturation relationships for two-phase flow is hysteretic in nature. We believe further work is needed to relate the properties of the fluids and other important factors (e.g., boundary condition and domain size) to the hysteresis in the τ-S w relationship. This is beyond the scope of this paper. Nevertheless, it is envisaged that the purpose of this paper, i.e., to investigate further the hysteresis in τ-S w relationships in homogeneous and heterogeneous domains through laboratory experiments, has been achieved. Tsakiroglou, C.D. and Avraam, D.G. (2002). Fabrication of a New Class of Porous Media Models for Visualization Studies of Multiphase Flow Process, J. of Material Sci., 37, 353. Tsakiroglou, C.D. and A.C. Payatakes (1998). Mercury intrusion and retraction in model porous media, Adv. Colloid Interface Sci., 75, 215-253. Tsakiroglou C.D., M.A. Theodoropoulou, and V. Karoutsos, "Nonequilibrium Capillary Pressure and Relative Permeability Curves of Porous Media." AIChE J. 49(10): 2472-2486(2003. Topp, GC, Klute, A and Peters, DB (1967). Soil Sci. Amer. Proc. 31, 312-314. Vachaud, G, Vauclin, M, Wakil, M. (1972). Soil Sci. Amer. Proc., 36, 531-532. van Genuchten, M.Th. 1980 Adamson and Gast (1997) a: Water-air system b: Silicone oil -water system 1