Scale dependent dynamic capillary pressure effect for two-phase flow in porous media

Causes and effects of non-uniqueness in capillary pressure and saturation (Pc–S) relationship in porous media are of considerable concern to researchers of two-phase flow. In particular, a significant amounts of discussion have been generated regarding a parameter termed as dynamic coefficient (τ) which has been proposed for inclusion in the functional dependence of Pc–S relationship to quantify dynamic Pc and its relation with time derivative of saturation. While the dependence of the coefficient on fluid and porous media properties is less controversial, its relation to domain scale appears to be dependent on artefacts of experiments, mathematical models and the intra-domain averaging techniques. In an attempt to establish the reality of the scale dependency of the τ–S relationships, we carry out a series of well-defined laboratory experiments to determine τ–S relationships using three different sizes of cylindrical porous domains of silica sand. In this paper, we present our findings on the scale dependence of τ and its relation to high viscosity ratio (μr) silicone oil–water system, where μr is defined as the viscosity of non-wetting phase over that of the wetting phase. An order of magnitude increase in the value of τ was observed across various μr and domain scales. Also, an order of magnitude increase in τ is observed when τ at the top and the bottom sections in a domain are compared. Viscosity ratio and domain scales are found to have similar effects on the trend in τ–S relationship. We carry out a dimensional analysis of τ which shows how different variables, e.g., dimensionless τ and dimensionless domain volume (scale), may be correlated and provides a means to determine the influences of relevant variables on τ. A scaling relationship for τ was derived from the dimensionless analysis which was then validated against independent literature data. This showed that the τ–S relationships obtained from the literature and the scaling relationship match reasonably well.


Introduction
A capillary pressure (P c ) and saturation (S) relationship is vital in the characterisation of immiscible two-phase flow phenomena in porous media (e.g., [2,[12][13][14]. This relationship is non-linear in nature and it depends on several factors, e.g., flow dynamics (steady or dynamic), fluid and porous material properties, contact angles, medium heterogeneities (e.g., permeability heterogeneity, non-uniformity in particle size), scales of observation, etc. [2,5,15,24,32,36]). In general, the forces which influence dynamic immiscible flow in porous media act in a manner so as to drive the fluids to a state of static equilibrium, i.e., when the rate of change of saturation (@S/@t is zero. During this process, the fluid/fluid interfaces in the porous medium move to re-establish a new state of equilibrium [3,11,21,34]. In practice, a state of stability in the system is utilised to obtain the equilibrium P c -S relationship for the system. In pore-scale modelling, the pressure in one of the phases is allowed to increase and a succession of equilibrium fluid configurations is computed in the pore network to obtain the P c -S relationship [3]. These are then applied to nonequilibrium conditions with fluids sometimes flowing at high flow rates. The implicit assumption in such an approach is that the disturbances to interfacial properties are rapidly dissipated [27]. At the equilibrium or quasi-static condition described above, the approach for the characterisation of the P c -S relationship is typically coupled with an extended version (i.e., an equation of two-phase flow in porous medium) of Darcy's law for single phase flow in porous medium [2,13]. In this case, the capillary pressure is defined as a function of saturation as expressed in Eq. (1). P nw À P w ¼ P c;equ ðSÞ ¼ f ðSÞ where, P nw = average pressure for non-wetting phase (ML À1 T À2 ) P w = average pressure for wetting phase (ML À1 T À2 ) P c,equ = equilibrium (steady state) capillary pressure (ML À1 T À2 ) S = wetting phase saturation (-) In the last decade, it has been argued by a number of authors that the application of Eq. (1) to non-equilibrium condition will be inappropriate as the equation lacks the requisite parameters to address the dynamic characteristics of the flow prior to the attainment of flow equilibrium when the rate of change of saturation (@S/@t) may be high. The need for a modification of the traditional relation (i.e., Eq. (1)) was the conclusion of some authors (e.g., [25,26,31]) with a view to accommodating a more complete description of the system under non-equilibrium conditions. Ever since its proposition, a phenomenal factor in the description of dynamic two-phase flow, namely a dynamic coefficient (s), has been the subjects of several publications [5,6,8,15,16,23,27,35]. The modification is mathematically expressed as in Eq. (2): where, P c,dyn = dynamic capillary pressure [ML À1 T À2 ] @S/@t = time derivative of saturation [T À1 ] While the interpretation of s in terms of its physical meaning seems to be well understood, some of the factors influencing its values are not uniquely characterised, i.e., the significance of the dependence of s on these factors seems to vary from one case to another. According to previous authors (e.g., [15,27] the dynamic effect is related to the dependence of capillary pressure-saturation relationship on the time derivative of saturation resulting from finite time needed by the fluid to neutralise the effect of the internal and external forces in order to establish flow equilibrium.
Other reports have indicated that the length of time needed to attain flow equilibrium and, hence, the dynamic effect, is larger in less permeable medium than in more permeable sample. The same is true for porous medium with higher degrees of microheterogeneity. Thus, larger values of s have been reported in the literature for fine grained or low permeability medium [10,15,23,42] and domains with micro-heterogeneities [17,33,35].
Among the other factors that affect the magnitude of s, strong impacts of fluid properties are well acknowledged [15,22,28,29].
The effect of scale/size of the domain has been linked to the s values as well. It seems that these effects are quantified more commonly in reference to numerical models of pore/core scales and application of different averaging techniques therein. However, these are rarely characterised through experimental studies. Dahle et al. [11] observe a scale dependency in the value of s in their numerical investigation of dynamic effects in P c -S relationship using a bundle-of-tubes model. Their model leads to a similar relationship as in Eq. (2) but with the inclusion of an intercept term, b.
By determining a dimensionless grouping containing s and b, respectively, they report that the dimensionless grouping shows a dependency on saturation and that s increases as the square of the averaging volume length, L. This suggests a geometrical increase in s with the length scale of the domain and that its magnitude may become arbitrarily large as we move from core to the field scale. This square of length relationship [11] was however controverted by Bottero et al. [5] in their investigation of non-equilibrium effects at different averaging windows in a porous medium domain of 21 cm height. Bottero et al. [5] determined s experimentally and obtained averages at different windows of observation, namely, 11 and 18 cm sections, within the same porous domain of 21 cm height. An order of magnitude increase in s was found at an averaging window of 11 cm as compared to the local values.
A similar order of increase was obtained in s at an averaging window of 18 cm as compared to the values at 11 cm averaging window. Thus, it can be inferred from the work of Dahle et al. [11] and Bottero et al. [5] that s generally increases with the scale of observation. The proportion of this trend to domain scale is however difficult to affirm from the current studies and there are clear inconsistencies in this regard. Das and Mirzaei [16,17] and Mirzaei and Das [36] found that the functional s-S relation follows non-linear trends at different locations within the same domain using a saturation-weighted technique for averaging domain scale s where local s increases with decrease in S. However, these authors did not find a clear relationship between the averaged and local values of s. Furthermore, insignificant differences were observed between the local and average (domain scale) values of s, which are in contrast to other findings where s is viewed as a parameter that increases with the scale of domain (e.g., [5,11]. Similarly, in the work of Camps-Roach et al. [8], an up-scaling procedure showed little or no change in s as its values remained almost identical for column scale as it was for the local measurements. Bottero et al. [5], who used averaging techniques such as those discussed by Nordbotten et al. [37][38][39] for averaging capillary pressure and relative permeability, found that the average values of s showed no dependency on water saturation within a water saturation range of 50% to 85%. This trend is contrary, within the similar saturation range, to what Das and Mirzaei [16,17] have observed on values of s and many other papers such as that by Dahle et al. [11]. Besides the above literature, there are a number of other relevant papers. For example, Bourgeat and Panfilov [7] reported a homogenisation method by which the dynamic capillary pressure effects were studied for oil-water flow in a heterogeneous porous domain. Cuesta et al. [10] reported a travelling wave solution for dynamic two-phase flow in porous media. It is evident from the literature that a significant amounts of discussion has been made to discuss the dynamic capillary pressure equation (Eq. (2)). From these discussions it can be concluded that the existing publications on the determination of s use both experiments and numerical simulations at pore or/and core scale [5,8,11,16,17]. They use various types of averaging techniques to determine averaged s at domain/larger scale and relate the local and averaged s values. While we acknowledge that the exact values of the coefficient may be different for specific two-phase system, the inconsistencies in the general trend in the results in the literature however raise many questions, e.g., are the observed behaviours of s an artefact of the averaging methods? Similarly, would the extent of the scale dependencies of s be any different if they were measured at different and separate domain scales and not averaging windows within the same domain? Considering the above arguments, one may opine that the scale dependency of a parameter, such as s in this case, should be determined using separate domains for the same conditions, e.g., material properties, boundary conditions, etc. As far as we know, there is no experimental evidence to demonstrate the scale dependency of s in this fashion.
To address these issues we aim to carry out a series of controlled laboratory experiments with separate and different sizes of porous medium with the aim of finding the trend in the magnitude of s depending on the scale. We aim to ensure that our procedure is less-dependent of mathematical averaging method for calculating the coefficient as far as possible. This approach has the advantages of assigning s to independent domain sizes, and reducing any ambiguity, which may arise from intra-domain system averaging. The experiments in this work rely on determining average quantities at the local scale (e.g., measurement sensors for water saturation, pressure). The scales of these measurements points are significantly small as compared to the size of the porous domains at which averaged quantities at the domain scale are calculated. In other words, the processes that govern the s values at the local and domain scale are likely to be different. We discuss this point further in the paper by relating the local and domain scale s values with the help of a s scaling factor.
In the context of this work, fluid property effects are important as well. Stauffer [43] emphasises the importance of fluid properties on s according to the following equation: where s is the relaxation parameter for unsaturated porous domain, a is a dimensionless parameter (defined to be 0.1), / is the porosity, K is the permeability, P d and k are the Brooks and Corey parameters related to the pore and entry pressure, and particle size distributions of the materials, respectively, l is the viscosity of the wetting phase, g is the gravity constant, q, the fluid density of the wetting phase. However, the relation does not account for the effects of the non-wetting phase viscosity and density. In the light of recent publications (see, e.g., [15,21,22,28,29], it is strongly established that the properties of the non-wetting fluid play significant influence on the displacement process and the magnitude of s. Also, the non-inclusion of the wetting phase saturation in Eq. (3)  Recently, Joekar-Niasar and Hassanizadeh [29] emphasised the dependence of s on viscosity ratio as well as an effective viscosity which depends on the values of fluid viscosities of the wetting and non-wetting phases. It was pointed out that if this ratio is less than one under drainage or greater than one in imbibition, invading front can become unstable. Conversely, the ratio that is far higher than unity in drainage or far less in imbibition results into a stable front. Viscous fingering is also said to result into less dynamic effect while stable front regime does otherwise [29]. We later show in this paper how the stable fronts arising from high viscosity ratios of the silicone oil-water system, used in our experiments for drainage process, results in higher dynamic effect as the viscosity ratio increases. The effect of viscosity was further explained by Joekar-Niasar et al. [30] and Joekar-Niasar and Hassanizadeh [29] in terms of the redistribution time needed by the fluid-fluid interface to reach equilibrium which takes longer in larger viscosity. In the same vein, higher density ratio is found to result into higher dynamic capillary pressure effect [15]. It has been observed in the literature that this promotes instability of fluid/fluid interfaces leading to the formation of fingers (e.g., [18,45]. The occurrences of instabilities and overshoots during two-phase flow in porous media have been shown theoretically by means of travelling wave analysis by Van Duijn et al. [46]. Fingering is also discussed in detail by Kissling and Rohde [40] and DiCarlo [19]). In the current work, the non-wetting fluids have significantly different viscosities but approximately the same densities. Thus, the effect of the density ratio will not be considered in the discussion. The experimental work of Goel and O'Carroll [22] reported the influence of viscosity ratio on the magnitude of s which increases with increasing viscosity of the non-wetting phase. This work utilised low viscosity non-wetting fluids, namely, 5 and 0.65 cSt silicone fluids. The authors show that the s values at various viscosity ratios collapse when normalised with effective viscosity especially for effective water saturations larger than 50%.
While the above discussion points to the importance of domain scales and fluid properties (especially, viscosity ratio) to determine the magnitude of s and its relationship to other process and system parameters (e.g., permeability), the relationship of the scale dependency of the s with viscosity ratio needs further studies. This is the main aim of this work. In particular, we aim to determine these relationships by a set of well-defined experiments with different domain scales and high-viscosity non-wetting phase fluids so as to increase the range of viscosity effects (200, 500 and 1000 cSt) studied so far. Indeed, this paper complements other published works on dynamic capillary pressure effects in relation to domain scales and fluid property. However, we expect that this work would provide a better understanding of the behaviour of the dynamic capillary pressure effect.

Porous medium, domain size and fluid property
The porous medium used in this work was prepared using sand particles (CH30) obtained commercially from Minerals Marketing Company (Cheshire, UK). The particles were compacted to prepare the porous samples following the same method of preparation and characterisation detailed in Das and Mirzaei [16]. The porous sample has an average porosity of 0.38. The intrinsic permeability of the sample is 5.66 Â 10 À11 m 2 which was measured using a constant head permeameter [2]. Silicone oil of different viscosities, namely, 200, 500 and 1000 cSt (Basildon Chemicals, Abingdon, UK) were used as the non-wetting phase while distilled water was used as the wetting phase. Silicone oil is a good non-wetting model fluid. Although their viscosities may be different their densities are almost the same. Therefore, any effect of the fluid properties on the results can be mainly attributed to viscosity difference. Properties of the fluids and porous media are shown in Table 1.

Experimental methods
Depending on the size of the domain, varying numbers of pressure transducers (PTs), which were bought from WIKA Instruments Ltd (Redhill, UK), were positioned at different levels of the three domains chosen for this work. At each level, we insert one pair of PTs where one PT is aimed at measuring the wetting phase pressure while the second PT measures the pressure of the non-wetting phase pressure. The PTs installed in our experimental cells included one pair of PTs at the centre of the 4 cm cell, two pairs at the 2 cm (top) and 6 cm (bottom) points of the 8 cm cell and, 3 pairs at the 2 cm (top), 6 cm (middle) and 10 cm (bottom) points of the 12 cm cell, respectively, where all distances are measured from the top boundary of the porous domains. membranes which were purchased from Porvair Filtration Group Ltd (Hampshire, UK) on the faces of the PTs to measure water and silicone oil pressures, respectively. They were held in place by capped fitting with a rubber seal. As the membranes are thin, they were supported by thicker Vyon filters approximately 5 mm thick (Porvair Filtration Group Ltd, Hampshire, UK). Prior to this, the membranes and filters were immersed in their respective fluid (water or silicone oil) and vacuumed for at least 24hrs to remove any entrapped air. The assembly of the PTs was also made within the respective fluid to avoid any air trap. The PTs were calibrated using a portable pressure calibrator, DPI 610 (Druck Limited, Leicester, UK). Good linearity was obtained in pressure-output voltage relationship and the linear equations of the fits obtained from Microsoft Excel (2010) were employed in the programming of the sensors. The PTs were connected to a data logger via a multiplexer (Campbell Scientific Ltd, Shepshed, UK).
Water saturation in the porous medium was determined using time domain reflectometry (TDR) probe (Campbell Scientific Ltd, Shepshed, UK). A TDR probe is positioned at the same level as with a pair of PTs. The TDR probes were earlier calibrated in a 4 cm cell and the results fitted with a polynomial function for in situ measurement of local water saturation corresponding to a level in a cell.
Transient average domain saturation is also determined using a graduated glass cylinder placed on a high precision weighing balance (±0.1 mg) which was connected to the computer and data logged in real time by the weighing balance software (A&D Company Limited, San Jose, USA). This method was adopted to determine average domain scale saturation as it represents the true saturation in the domain at any particular time from a material balance point of view. A similar approach was adopted by Bottero et al. [5] to determine the average (column-scale) fluid saturation over the whole sand column from the change of volumes of fluids in inflow and outflow burette. Having determined the porosity and the total amount of water in the porous domain with the dead volume of the system subtracted, the transient domain saturation was determined.
The sample holder was set on a metal base with wire gauze and hydrophilic filter at the base of the column. The fine wire gauze was used to prevent fine sand particles from slipping into the outflow line which might affect desaturation rate value while the hydrophilic membrane acts to prevent oil outflow. A predetermined amount of water was then poured into the cell to a certain position followed by pouring of sand through a metal sieve of appropriate size to ensure uniform sand deposition and prevent air trap. The cell was vibrated simultaneously to promote uniformity in the deposition of the sand particle. After the deposition of the sand particles was complete, the excess water was drained through the outflow valve. This gave an indication of the amount of water in the saturated domain. Hydrophobic filter was then placed over the sand to prevent water flow out through the top of the domain boundary. Finally, the column was sealed with a metal cap which provided the connection to silicone oil supply.
The silicone oil was supplied from a Marriote bottle which maintained a constant pressure on the fluid column irrespective of fluid volume. The Marriote bottle was directly connected to an air compressor (R.E.P. Air Services, Loughborough, UK) which provided the requisite pressure for injecting silicone oil in the experimental cells (also discussed in the next section).

Dynamic and equilibrium experiments on two-phase flow
As mentioned earlier, the two-phase flow experiments were conducted on 4, 8 and 12 cm cells. Pressure was supplied using  Table 1.
an air compressor which imposed a set air pressure on a column of silicone oil in a Marriote bottle (see Fig. 1). The bottle ensures that the oil invades the porous medium at a constant pressure throughout the experiment. For the dynamic drainage, a set pressure was imposed on the oil which pushed it to the top of the domain through the inlet valve. When the outlet valve remains closed, the pressure simply transmits through the fluid phases without infiltration of oil into the sand domain. This was indicated by the readings of the PTs which rose gradually till it attained approximately the same magnitude as the imposed pressure at different levels in the domain.
For the dynamic drainage processes, 10, 15 and 20 kPa pressures were imposed on top of the domains. These provide us results for two scenarios. Firstly, the s values can be calculated using a number of the pressure conditions (three values in this case) as has been done in a number of papers (e.g., [5,16,17,23]. Secondly, the chosen pressure boundary conditions provide us the possibility of obtaining s values for selected cases where the pressure gradient is the same in two domains of different size. For example, the pressure gradients are the same in 4 and 8 cm cells when pressures of 10 and 20 kPa are imposed, respectively. Similarly, the pressure gradients in 8 cm and 12 cm cells are the same when 10 and 15 kPa are imposed. These issues are discussed again in a latter section (Section 3.2).
When the outflow valve was released in our experiments, drainage of the porous domain began and it gave noticeably different P c -S profiles for different boundary conditions. Replicate runs of some the experiments were conducted to check repeatability of the experiments. The ambient laboratory temperature was at approximately 20°C.
The experiments in the cells were repeated for 200, 500 and 1000 cSt silicone oil. The equilibrium or quasi static experiments were conducted by gradual increase in air pressure imposed on the silicone oil in the Marriote bottle at less than 500 Pa in a single step. Sufficient time was allowed for the fluids to attain equilibration before the pressure was further raised. This was continued until the porous medium saturation could not be reduced anymore which we adopt as the final measured saturation. We adopt this approach so as to avoid any confusion of the final measured saturation with the typical definition of irreducible saturation in porous media. In our case, the final measured saturation is not necessarily the irreducible saturation and it simply reflects what final saturation has been obtained for the experimental conditions (e.g., pressure drop).
For each silicone oil and domain scale chosen, dynamic and equilibrium or quasi static drainage experiments were conducted. As stated earlier, the local saturation was determined using the TDR probe for the local level while an average domain saturation was determined from a graduated cylinder placed on a weighing balance whose software program was set at the same output rate (15 s) as the data-logger for sensors to ensure corresponding relation of generated data. The rate of change of saturation in this cylinder is taken to represent the whole domain desaturation rate while the capillary pressure for the whole domain was derived from saturation-weighting of the capillary pressures at different local levels in the system. Thus, at a particular output time, the domain saturation is determined from the balance measurement while the domain capillary pressure is determined from saturation-weighting of the TDR and PTs measurements at the local levels.

Calculation of dynamic coefficient (s)
The dynamic coefficients were determined for the different cases using Eq. (2). To do this, the equilibrium and dynamic drainage pressures of each phase were interpolated at selected saturations using the 'Forecast' function in Microsoft Excel (2010). The function calculates, or predicts, an unknown parameter by using known values on the basis of a linear regression. The predicted values are y-values for given x-values and the results can be plotted easily as a typical x-y graph in Excel. This readily provides the corresponding points for the parameters in the equation, i.e., P nw and P w at a particular saturation, S, for both the dynamic and equilibrium drainage experiments. The differences in the equilibrium and dynamic P c data which are required for Eq. (2) were obtained at arbitrarily chosen saturation values but covering the entire saturation range, i.e., 0-1. Similarly, the desaturation rate data were interpolated using the same function. The plots of the differences between the dynamic and equilibrium capillary pressures against the desaturation rates at the respective water saturation for different imposed conditions were fitted with a straight line using Microsoft Excel (2010). The slope of this line gives the s values. . Most papers in the literature have used linear interpolation to determine the dynamic pressure difference at a specific saturation (e.g., [5]. Also, from the saturation and time data, desaturation rate (@S/@t) was determined using the central differencing scheme. This is the most popular approach in the literature (e.g., [5,15,36]. The capillary number, lm/c, determined for the system at viscosity ratio of 200 was 1.14 Â 10 À6 at the highest superficial velocity, which was determined from the highest flow rate recorded at the outflow volume collected in bottle on the weighing balance. l is the viscosity of the oil, m is the superficial velocity and the c is the interfacial tension between the oil and water. As explained in Section 2.4, while the saturation at local levels was determined from calibrated TDR probe readings, that of the whole domain was determined from water outflow measured on a weighing balance. Taking density of water at ambient condition as 1 g cm À3 the volume collected was readily determined and the domain saturation was calculated based on difference in the initial volume of water in the porous medium and the outflow. The saturation-weighting of the local levels P c provides the domain level P c at every count of time, t n [s] (see Eq. (4)). Similar to the s determination at the local levels, Eq. (2) was used to determine s for the domain.
The numerical simulations by Hanspal and Das [23] have reported the effects of temperature on dynamic capillary pressure. At the moment, there are little experimental studies on the role of temperature on dynamic capillary pressure effects and in particular the temperature effects on the scale dependent s-S relationships. Our work in this paper does not concern the effects of temperature and therefore we measure the s only at ambient conditions. But the temperature dependency of the dynamic capillary pressure effects is an interesting topic and should be explored in detail in the future.

Averaging approach
The averaging technique to determine s for the whole domain utilised the following saturation weighted relation for capillary pressure: where P c domain j tn represents the domain representative capillary pressure at a particular saturation and the corresponding experimental time, t n [s] referring to the time at nth count of data generated, P nw [Pa] and P w [Pa] are the non-wetting and wetting phase pressures at the nth count of the generated data. S [-] refers to the wetting phase saturation and j = 1, 2, . . . m, with m being the total number of measurement levels in the domain i.e., 1 for 4 cm, 2 for 8 cm and 3 for 12 cm cells. Although we do not track the fluid-fluid interfaces within our experimental domain, in principle the higher viscosity of the non-wetting phase may result in sharper saturation front of the interfaces. To take care of this possible situation in the calculation of domain averaged P c the following approach was taken. Prior to the arrival of the oil phase at the lower measurement points in 8 and 12 cm high columns, the first term of Eq. (4) is set to zero. At these points, the calculated P c values are negative and are discarded.
In general, a sharper fluid-fluid interface is thought to have greater influence on when the interface reaches the measurement points and, hence, it may affect the calculated P c . However this was not a concern in this study. As evident from Eq. (4), the domain scale P c is calculated using the local water saturation in the porous medium determined from the TDR probes and not the saturation determined from outflows. Furthermore, the average P c is only determined when the oil phase reaches all the measurement points. In this way, we avoid any discrepancy that may arise in the equilibrium P c -S data at the domain and local levels. Such an approach has been adopted in a number of earlier studies (e.g., [5,17] and we continue the same tradition in this paper.
In addition, the following central difference scheme was used to obtain the desaturation rate for the entire domain based on the outflow volume collected in the cylinder placed on the balance: This is a well-accepted method and is frequently used in the literature [5,16].

Results and discussions
Determination of s-S relationships for two-phase flow requires prior determination of P c -S and @S/@t-S relationships for the system of interest. However, much like the s-S curves, the P c -S and @S/@t-S relationships are subject to the significance of domain scale and viscosity ratio effects. Imposition of the same boundary pressures on dissimilar domain scales lead to different impacts of pressure drop (DP) across different domain scales with decreasing magnitude of DP as the domain length increases. On the other hand, different pressures can be imposed on top of different size domains to result in the same pressure gradient across the domain. The impact of these effects on the P c -S and @S/@t-S relationships and the eventual s-S relationships are the objects of the discussion in this section. At the end, the parameters of importance identified in this paper together with relevant papers from the literature form the basis for the development of a dimensionless correlation for s using Buckingham's G-theorem.
In order to assess the repeatability of the experiments, replicate experiments were conducted which showed reproducibility under similar conditions for both the P c -S relationship as well as the desaturation rate. Fig. 2 shows a typical replicate @S/@t-S curves for the 4 cm cell and 200 viscosity ratio oil at 20 kPa imposed pressure. The repeatability indicates consistent packing of the porous medium in the cell. Viscosity ratio (l r ) was earlier defined as the ratio of the viscosity of the non-wetting phase to that of the wetting phase. In this paper, viscosity ratio and viscosity can be used interchangeably since the kinematic viscosity of water is taken to be 1 cSt. Thus, the viscosity ratio has the same magnitude as the viscosity of the non-wetting phase.
From Fig. 2 it can be seen that the desaturation rate rises fast to a high value (in the absolute sense) and declines steadily until the saturation approaches the final measured saturation. Under dynamic drainage, as the outflow from the domain continues, the desaturation rate rises to the highest value at a high saturation. This is followed by a declination as the rate decreases henceforth towards the final measured saturation. Similarly, other equilibrium and dynamic drainage experiments conducted in duplicates showed good reproducibility. All of our results show that the dynamic P c -S curves lie above the equilibrium or quasi static curves and it is uppermost for the largest imposed pressure i.e., 20 kPa. This goes to affirm the presence of dynamic capillary pressure effects and it is consistent with the findings from other authors [9,15,16,41,44].

Scale dependency of P c -S and @S/@t-S relationships for different fluid viscosity ratios
Prior to the determination of s, we discuss the P c -S and @S/@t-S relationships briefly as they form the bases for the values of s.

P c -S relationships
The trends/patterns of the P c -S relationships in this work are consistent with what have been reported by other authors [2,15] and therefore, they are not discussed in length. As an example of the P c -S curves obtained in this work Fig. 3 is included, which provides a graphical picture of the relationships for different silicone oils for the 4 cm high porous domain. In consistence with other studies, the capillary pressure rises as the water saturation reduces and it becomes almost upright as the saturation approaches the final measured saturation. Also, the dynamic drainage P c -S curves for 10, 15 and 20 kPa lie above the quasi-static drainage curve. This is the reflection of dynamic effect and is in line with the desaturation patterns. Also, worthy of note is the increase in the P c with viscosity increase. The required entry pressure is seen to similarly increase with l r . This is owing to the resistance of the porous domain to flow which increases with l r . Thus, a l r of 200 has the least resistance to enter the medium while 1000 l r has the highest. At comparable saturation, 200 l r silicone oil-water system has capillary pressure that is lesser than 500 cSt while 1000 cSt has the highest. Similar patterns are observed for the bigger domains.

@S/@t-S relationships: influences of fluid properties and porous domain characteristics
The desaturation rate (@S/@t) of the aqueous phase is one of the important factors, which determine the magnitude of s [16,17].
Therefore, a good understanding of @S/@t and its relationship with the saturation for different fluid properties (e.g., viscosity ratio) and domain scales are crucial in quantifying the scale dependency  Table 1.
of s. In this paper, the time dependent saturation data corresponding to the dynamic P c -S curves are used to determine the @S/@t-S relationships. When @S/@t = 0, the system reaches flow equilibrium for particular conditions (e.g., boundary conditions, porous medium properties). As a point corresponding to @S/@t = 0 can be easily identified in the @S/@t-S plots, these graphs can be used to distinguish between the equilibrium and non-equilibrium points in terms of saturation. In a number of previous studies @S/@t-t curves prepared from dynamic S-t curves were adopted to differentiate the equilibrium and non-equilibrium points (e.g., [17,36]. However, we adopt @S/@t-S curves in this work. We vary a number of parameters in this paper (e.g., pressure, domain size, viscosity) which affect both the saturation and time to equilibrium. However, it is the saturation values which are used later in this paper to discuss the saturation dependent s for different cases. In order to have a consistent framework for relating the @S/@t to s where we have saturation as the common parameter we chose to use @S/@t-S and not the @S/@t-t curves in this paper. We have noted with interest that Hou et al. [28] have chosen to work with @S/@t-S curves as well.
We find that @S/@t decreases as the domain size increases with a visible effect of the imposed pressure on it along the depth of the domain. This is depicted in the typical results in Fig. 4(A-F) for l r of 200. Starting from 8 cm cell in Fig. 4(A), we find that the top section of the domain has the highest desaturation rate. The bottom section has a lower desaturation rate in comparison while the whole domain scale has the least rate of desaturation. It should be recalled that the 8 cm cell is subdivided into two sections of 4 cm each. Each has the TDR probe for in situ measurement of water saturation. Similar patterns are seen in Figs. 5 and 6(A) for 500 and 1000 l r , respectively. Furthermore, Fig. 4(B) presents the @S/@t-S relationship in the 12 cm cell for 200 l r . Similar to the pattern in 8 cm cell, @S/@t at a particular saturation decreases along the depth of the domain. The top section has the highest desaturation rate followed by the middle section. The @S/@t is the least at the bottom section and remains the closest to that of the whole domain. The pattern remains similar in Figs. 5 and 6(B) for 500 and 1000 l r , respectively. It is noticeable that the @S/@t-S profiles at the middle and bottom sections are close to each other.
The fact that @S/@t changes along the depth of an experimental cell suggests that there are likely to be sharp saturation fronts and piston-like flow in the domain particularly at higher viscosity ratio. These issues were discussed recently by Hou et al. [28] using numerical simulation. Sharp saturation fronts, fingering and piston-like flow in the porous domains have been discussed also in a number of other studies as we have mentioned earlier in this paper. In our experimental studies, we do not explicitly monitor the saturation fronts in the domain, and therefore, we avoid a detailed discussion in this regard. We observe that the variations in the desaturation rates at particular saturation for different cell sizes generally reduce as the viscosity ratio increases.

Domain scale dependency of s-S relationships for different fluid viscosity ratios
As earlier stated, the forecast function in Microsoft Excel (2010) was utilised for linear interpolation in order to establish corresponding points of saturation and capillary pressure. This was done for both the quasi-static and dynamic drainage conditions as well as to determine the desaturation rate.
The estimated differences in dynamic and equilibrium capillary pressures (P c,dyn À P c,equ ) were fitted to a straight line against the desaturation rate (@S/@t) at chosen saturation values for different imposed boundary conditions. Two possible curve fitting approaches were tried. First, a direct fit of the data points which may generate an intercept on the y-axis was used. Secondly, an approach which force the fit through the origin by assuming an additional point at 100% saturation when @S/@t = 0 is chosen.
Slopes of these straight lines represent the s. The first approach is commonly used in the literature [11,16,33,36]. s results from the first approach (intercept fit) are presented in this work and the comparison of the two approaches is shown in Fig. 7.
It can be seen that the two approaches produced similar trends that only differ at some middle saturation range. In other words, the fitting method might not influence the calculated s at high   Fig. 3 for dynamic conditions. Domain scale @S/@t obtained from balance readings and local @S/@t from TDR readings. and low water saturation values. In any case, we use the first approach in this work to calculate the s values in consistent with most studies, unless otherwise stated.   Fig. 3 for dynamic conditions. Domain scale @S/@t obtained from balance readings and local @S/@t from TDR readings.
at a particular saturation. The interpretation/implication of increased or decreased s values from the physical point of view has been discussed extensively [22,27,35], and the same holds true in this case. For example, when s increases, the system requires higher capillary force (or time) to reach a state of flow equilibrium.
A graphical comparison of the results for s-S curves at different viscosities in 4 cm cell is shown in Fig. 8 Fig. 8(C) shows the relationship at 12 cm scale for 500 l r . The effect of pressure gradient is present but it seems to be lesser in impact than that for 200 l r as the domain scale values of s lie close to the values at the bottom level. Similar effect of scale had earlier been expressed by Bottero et al. [5] where they affirmed an order of magnitude increase from local to upscaled or domain scale.  [22] used a pressure cell of 10 cm diameter similar to our cell of 10.2 cm diameter but the length of their cell was 20 cm compared to ours of maximum height 12 cm. This shows that the effect of height can be very significant on the results. Furthermore, intrinsic permeability of their medium (1.53 Â 10 À11 m 2 ) is lower than that of this work and this is known to influence the magnitude of s [15]. In addition, the influence of @S/@t can still be found in the reported values. For example, Goel and O'Carroll [22] reported 9.8 Â 10 À4 and 2.9 Â 10 À3 s À1 as the highest @S/@t values for the 4.098 and 0.442 viscosity ratios, respectively, as compared to 4.1 Â 10 À3 and 3.4 Â 10 À3 s À1 at 4 and 12 cm domain scale, respectively, for 200 l r in this work. While the values are closely of the same order, the gradient in the @S/@t values which might result from the pressure gradient effect discussed above will be much higher in their work considering the length of the domain used in the work (20 cm). Thus, higher values of s might be reported than expected.
Also, Camps-Roach et al. [8] reported s in the range 3.7 Â 10 5 to 10 6 Pa s for air-water system for a cell of 10 cm diameter and 20 cm height. This can be compared with our results discussed above and the same reason can be adduced for the closeness of s despite a large difference in viscosity values. In addition, Das and Mirzaei [16] utilised similar viscosity ratio (200) as in this paper and the same 12 cm height cell, their results for s range between 10 3 to 10 6 Pa s. This is lower than the corresponding viscosity ratio and scale in this paper but the difference can be explained in the light of higher permeability medium used by Das and Mirzaei [16] which was 8.7 Â 10 À10 m 2 . By and large, our results in this paper are comparable to other studies although they are for different fluids pair and domain sizes.
In the light of the above discussion, one can discern the effects of scale and viscosity ratio on the trend in the magnitude of s. The rise in the value of s with viscosity can be attributed to decrease in desaturation rate with increased viscosity of the non-wetting phase. In other words, s value is inversely affected by the desaturation rate (@S/@t), i.e., the lower the desaturation rate, the higher the s. Thus, the value of s rises as the l r increases from 200 to change with the averaging window considered in their work. However, our observation is similar to the report by Bottero et al. [5] with respect to the domain scale. They [5] reported an order increase in up-scaled values of s for their experimental domain in comparison to the local level. Similarly, they observed that magnitude of pressure difference at local scale (0.7 cm) was close to that at the averaging window of 11 cm while the desaturation rates were more apart with 4 Â 10 À2 s À1 and 4 Â 10 À3 s À1 for the respective scales in their experiment. Thus, they concluded that the scale dependency of s is mainly because of the rate of change of saturation, whose magnitude decreases with the increase of the length scale. To explain this, we can refer to Eq. (2) where s is related to both the pressure difference and the desaturation rate. While the pressure difference can be fairly similar under different conditions the @S/@t can become drastically low under many conditions, e.g., viscosity ratios of the two-phase system, distance from injection point, porous media properties and so on. In contrast to the observation of Dahle et al. [11], s may not become arbitrarily large as the averaging length increases as observed in the above discussion. In Fig. 8, three different pressures (10, 15, 20 kPa), which may impose the same or different pressure gradients (DP) across the three domains for different domain sizes, have been used to determine the s-S relationships as typically used in the literature [17].
In order to understand the results better and to further confirm that the s-S relationship depend on the scale, the domain scale dependence of s-S curves using similar DP across the domains are presented in Figs. 9 and 10. The calculation procedures for s in this case follow the same approach as described by Goel and O'Carroll [22], Camps-Roach et al. [8], and Sakaki et al. [41]. Briefly, the approach involves the use of a single pressure head for dynamic drainage to obtain the P c,dyn in Eq. (2). The difference between the dynamic and the quasi-static P c (P c,dyn À P c,equ ) is then plotted against the @S/@t values at corresponding wetting phase saturation (S). The slope of this straight line is then defined as the dynamic coefficient. This approach is different from that employed in other experiments including the results so far in this paper where three different pressure heads are imposed for dynamic drainage and the best linear fit is drawn through the plots of ''P c,dyn À P c,equ '' against @S/@t at corresponding saturation (e.g., [17]). The results of our calculations with two similar pressure gradients are shown in Figs as the pressure gradient is increased (doubled), s increases again for a given domain size. The capillary force (or time) that a system needs to reach an equilibrium state is related not only related to the domain size or viscosity ratio or both but also the pressure gradient. As a result, we observe that as the pressure gradient across a domain is increased, s is also increased provided all other factors remain the same (Figs. 9 and 10(C)). Furthermore, the results in averaging scales and any term that affects ðP c;dyn À P c;equ Þj S =@S=@tj S should also affect the ratio. As said before the local s value depends on the processes that occur at a measurement point in the domain. Further, it is the highest for the domain scale. Thus, this ratio should be equal to or smaller than unity at all saturation. As the observation volumes for the measurements sensors (i.e., TDR probes for saturation and the pressure transducers for capillary pressures) are significantly small as compared to the size of the full domain, there may be some differences in the measured equilibrium and dynamic capillary pressures at these two observation scales. The processes which dominate the values of s in the local or domain scales are likely to be different. For the local scale s one should expect that the pore scale processes (e.g., dynamic contact angles, inertia effects in the pores, viscous effects captured in capillary pressure) are the major influence in its values. In the averaged s at the domain scale one should expect to find the s values which are governed mostly from spatially averaging of non-linear processes (e.g., pressure and saturation distribution).
To indicate the relationship between the s values at the local and domain scales, we define a s scaling factor, s r , which is simply the ratio of s at a particular measurement point to that at the domain scale, both measured for the same porous domain, fluid pair and medium permeability. We then plot s r as a function of S. This is expected to provide a better understanding of the scaling dependency of s in relation to viscosity ratio. For example, when s r is close to unity, it means that the significance of the dynamic capillary pressure effect at the local and the domain scale are similar. To obtain the s values for this section, we use all the pressures (10, 15, 20 kPa) as we have used to obtain the results in Fig. 8.
For the 4 cm high cell, there is only one measurement point at the centre of the domain, thus this scale is excluded from this discussion. For the 8 cm high cell with two pairs of measurement points, s r is higher at the bottom, than the top measurement point for 200 l r . This is illustrated in Fig. 11(A). The difference between s r at these two points appears to be nearly constant at higher water saturation but it becomes more pronounced as the saturation decreases. This implies that significance of the scale dependent s-S relationships may be different if compared at different measurement points within the same domain.
The relation of s r for 500 l r for the same domain scale is shown in Fig. 11 Fig. 11(C). The s r trend in the 12 cm high domain follows similar pattern as above. This therefore confirms the general observation we made for s r . In addition to the top and bottom measurement sections, 12 cm cell has a middle measurement point. Fig. 12(A) shows that s r remains highest at the bottom section for 200 l r . This trend is further affirmed in Fig. 12(B) and (C) for 500 and 1000 l r respectively. However, it can be seen that the ratio is the highest for 500 l r . This is connected to the fact that s values at the measurement points for 500 l r in 12 cm high cell are closer in value to the s at the domain scale than other viscosity ratios for the same 12 cm high cell. Since s is related to how far or close to equilibrium is the two-phase system [15], above behaviours can be explained   in terms of the equilibration times [27] which becomes larger as the depth of the domain from the injection point increases. Since s is higher at the bottom of the domain, it implies higher equilibration time is required at a point far from the point of fluid injection (the bottom point in this case). However the highest equilibration time will still be required at the domain scale.
3.4. Non-dimensional analysis of dynamic coefficient using the Buckingham's P theorem As discussed in a number of other publications and this work, many factors are attributed to the presence of dynamic capillary pressure effect in two-phase flow in porous media. In the light of the results in this paper the influence of domain scales is elucidated while the effects of viscosity ratio are emphasised. Overall, it can be safely stated that it is a combined effect of these effects that determines the magnitude of s and not simply the scale of the domain or the viscosity ratio or the medium permeability. Therefore, relating these variables through a functional form provides the advantage of applying this relation to directly estimate s and in particular to quantify the interplay of variables in determining the s value in systems varying in scales or other properties.
Such a functional relationship can be derived from a non-dimensional analysis of s and the relevant variables highlighted above.
This can derive from the fact that all variables in the equation should have zero dimensions. By constructing non-dimensional groups of variables expressing a relationship among different variables, it is then possible to express empirical data in the form of a function. More importantly, the lumped behaviour of a parameter (i.e., s in this case) in domains of different size can be understood.
When a dimensional analysis is applied, if the relevant dimensionless groups are the same on two systems, e.g., a smaller and a bigger porous domain, the result of the smaller system would be applicable to the larger system or vice versa. Therefore, conceptually a non-dimensional analysis of the dynamic coefficient provides an attractive way to correlate the experimental results in this work and determine further how it behaves in combination with other variables.
To carry out the non-dimensional analysis of s, the well-known Buckingham's G-theorem was employed owing to its usefulness in determining sets of dimensionless parameters from chosen variables. In principle, the theorem reduces a dimensionally homogenous relationship (i.e., every independent, additive term has the same dimension) involving n variables in m fundamental dimensions to a single relationship among n À m independent dimensionless products [20]. Thus, using the G's notation to represent the dimensionless term, the theorem can be expressed as: where G represents dimensionless term. In this work we applied a common method namely the 'method of repeating variables' to determine the G terms, where the following variables have been used based on our knowledge on the dynamic coefficient. f ðs; g; K; P d ; V; k; q r ; l r ; S; /Þ ¼ 0 Here s is the dependent variable and g; K; P d ; V; k; q r ; l r ; S; / are the independent variables. In Eq. (6), g is the gravity, K, the isotropic intrinsic permeability, P d , the entry pressure, /, the porosity, k, the pore size distribution index, V, the domain volume representing domain scales, q r and l r are the fluid density and viscosity ratios, respectively, while S is the water saturation. The variables are chosen based on the most significant variables that have been reported to be important in determining the s. We define that the effects that boundary conditions or any other parameters (e.g., pressure gradient) have on the results are directly captured by the values of s and S.
In the analysis, it is necessary to use the dimensional terms for forming the dimensionless groups. Therefore, we begin the analysis by choosing the five dimensional variables in Eq. (6), namely, s, g, K, P d and V. We can incorporate the dimensionless variables (e.g., /, S) in the derived functional form as convenient but based on the knowledge of how a particular variable affect the s.
Using the standard theorem procedures, it can be shown that two dimensionless groups of variables are needed to form a functional relationship between the above variables which are as follows: P 1 contains the dynamic coefficient, s, and it can be interpreted as the dimensionless group of s. P 2 can be viewed as the nondimensional volume of the domain (scale) which is made dimensionless with respect to the permeability of the porous domain. P 2 can be modified to include the non-dimensional parameters in Eq. (6) without affecting the dimensionality of the group. This requires some user's experience and knowledge of the relevant process. However, we have significant understanding of how the dynamic coefficient generally depends on these variables, i.e., s is directly proportional to viscosity and density ratios while it is inversely related to saturation, porosity and pore size distribution index. Using these general understanding, we modify Eq. (8) to include the non-dimensional independent parameters as follows: Expressing this in a non-linear form: where a and b are the correlation coefficients to be determined from experiments for specific case. Eq. (10) shows the dimensionless groupings which describe an interplay of variables that affects s. The influences of saturation and domain scale are not clearly discernible in Fig. 13(A) as each viscosity attains its final saturation, expressed in P 3 , differently irrespective of other terms. However, the impact of viscosity is visible as P 1 increases with increasing viscosity ratio for the same value of P 3 . It is also visible at some points that P 1 for 500 l r lies higher than for 1000 l r at the same P 3 . This may appear contrary to expectation, but its explanation lies in the fact that at the same P 3 the values of the two variables, V and S, are not the same for each of the viscosities. At this same saturation point, s value in 500 l r correspond to a certain V and S while in 1000 l r , s at this point correspond to lower V and higher S, thus giving P 1 a lower value in 1000 l r .
In order to calibrate the general scaling relationship between P 1 and P 3 (Eq. (10)) the data points for all viscosities and scales in this work are combined with P 1 and P 3 calculated from some typical literature data, namely, the data from Das and Mirzaei [16]. All of these data points are grouped together in Fig. 13(B). The best fitted curve through these points provides a scaling relationship as shown in Eq. (11).
Eq. (11) shows a power law correlation among different system variables, e.g., dimensionless s, dimensionless domain volume (scale) and the viscosity ratio for the system we considered. We appreciate that a scaling relationship such as the one shown in Eq. (11) may be derived using different forms (i.e., both power and non-power law behaviour). However, our observations on the data presented in Fig. 13 suggest the existence of an exponent and proportionality constant between the two dimensionless groups (P 1 , P 3 ) in consistent with what one would expect in a power law function. Power law functions are well understood as they occur in most studies in some form and they are useful to interpret data easily, e.g., these functions can be plotted as straight lines using logarithmic scales. Keeping these points in mind, we have used a power-law function to determine the scaling relationship for the non-dimensional form of dynamic coefficient (Eqs. (10) and (11)). To demonstrate the validity and applicability of the scaling relationship (Eq. (11)), a set of independent data from Bottero [4] were tested. Fig. 13(C) shows a plot of the dimensionless groups (P 1 and P 3 ) calculated from the data of Bottero [4]. As evident in Fig. 13(c) these data are reasonably well predicted by the scaling relationship for the same parameters used by Bottero [4]. In the figure, only one significant outlier could be found in comparison of the P 1 -P 3 curve from Bottero [4] and the predicted P 1 -P 3 relationship using Eq. (11). An attempt was also made to predict the s-S relationship for the experimental data of Bottero [4] using Eq. (11). The comparison of the original [4] and predicted s-S relationships are shown in Fig. 13(D) which shows a similar trend in the data. Again, only one significant outlier occurs at the saturation value of 0.8. These results seem to show that a scaling relationship such as the one in Eq. (11) has useful general application in the prediction of dynamic capillary pressure effect in terms of the dynamic coefficient. This paper has demonstrated that the dynamic coefficient depends on a number of parameters (e.g., domain size, viscosity ratio). The paper further shows that by constructing non-dimensional groups of quantities expressing a relationship among these different variables (i.e., P 1 and P 3 ), it is possible to summarise experimental results and determine their functional relationship. This approach has an advantage that it is possible to find out how many dimensionless groups are required to replace the variables which affect the dynamic coefficient. The empirical data obtained in this work along with the scaling relationship involving the two dimensionless groups suggest that it is possible to determine how significant the dynamic capillary pressure effect is in a particular situation in terms of the dynamic coefficient.

Conclusions
In this work, the effects of fluid viscosity ratio and scale of observation on the magnitude of dynamic capillary pressure (indicated by dynamic coefficients) have been elucidated. We have analysed two types of scale dependencies of the s-S relationships in this work, which relate to (i) the scale dependencies between small and larger domain and, (ii) the scale dependencies between the local measurement point and domain scale. It is confirmed that s increases with the scale of observation and the viscosity ratio of the fluid-fluid system. We also found that its value is inversely affected by @S/@t, which is connected to the degree of resistance to the fluid motion, namely, viscosity. In almost all cases, s is found to decrease monotonically with increase in S although we find some instances where the s-S relationship is non-monotonic.
An order increase in magnitude of s was observed as domain scale changes from 4 cm scale to region of low saturation in 8 cm domain scale. Similar order of increase with higher magnitude was recorded in 12 cm domain scale. Our results showed that the relationship between s and domain scale depends further on viscosity ratio. Also, difference in pressure gradients in different domain sizes acts to reduce the desaturation rate which in effect causes a rise in the magnitude of s.
In addition, the fluid viscosity ratio together with the experimental scale results in an order increase in the value of s as the scale increases from smaller to larger domain for the same viscosity. While a higher order increase is found with increase in viscosity ratio from 200 to 500, the change in s while l r changes from 500 to 1000 is found to be of the same order. It can be reasonably inferred that the scale dependency of the dynamic coefficient is mainly because of the rate of change of saturation, whose magnitude decreases with the increase in depth of the domain from the injection point. This can also be said of s dependency on viscosity ratios of the two-phase fluids. As the viscosity ratio increases, the desaturation rate reduces and hence s increases.
The dimensionless correlations developed for the system showed the applicability of power law for the system and can be found useful in the first hand estimation of dynamic coefficients for system of interest provided the system parameters are available. In short, the magnitude of s is similarly affected by the coupled effects of viscosity and domain scales and s may not become arbitrarily large as the averaging length increases.