Dynamic Effects on Capillary Pressure–Saturation Relationships for Two-Phase Porous Flow: Implications of Temperature

Work carried out in the last decade or so suggests that the simulators for multiphase flow in porous media should include an additional term, namely a dynamic coefficient, as a measure of the dynamic effect associated with capillary pressure. In this work, we examine the dependence of the dynamic coefficient on temperature by carrying out quasi-static and dynamic flow simulations for an immiscible perchloroethylene–water system. Simulations have been carried out using a two-phase porous media flow simulator for a range of temperatures between 20 and 80 (cid:2) C. Simulation domains represent 3-D cylindrical setups used by the authors for laboratory-scale investigations of dynamic effects in two-phase flow. Results are presented for two different porous domains, namely the coarse and fine sands, which are then interpreted by examining the correlations between dynamic coefficient(s) and temperature, time period(s) required for attaining irreducible water saturation, and the dynamic aqueous/nonaqueous phase saturation and capillary pressure plots. The simulations presented here maintain continuity from our previous work and address the uncertainties associated with the dependency of dynamic coefficient(s) on temperature, thereby complementing the existing database for the characterization of dynamic coefficients and subsequently enabling the users to carry out computationally economical and reliable modeling studies. V V C 2012 The Authors. AIChE Journal, published by Wiley on behalf of the AIChE. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. AIChE


Introduction
Flow of immiscible fluids in porous media is of special importance in soil science as well as in chemical, environmental, construction, and petroleum industries. In particular, the immiscible flow of nonaqueous phase liquids (NAPLs) commonly occurs in various subsurface conditions, groundwater remediation operations, and in oil recovery processes. In such systems, capillary pressure plays a crucial role in determining the motion of fluids within the porous media. [1][2][3][4][5][6][7][8] In addition, temperature variation may affect such flows where the hydraulic properties of the porous media such as hydraulic conductivity and water retention are temperature dependent. [9][10][11][12] Conventional approaches for modeling of immiscible two-phase flow in porous media involve the use of an extended version of Darcy's law 3,13,14 for multiphase flows in conjunction with the use of constitutive relationships between capillary pressure, saturation, and relative permeability (P c -S-K r ) based on quasi-static conditions. 15,16 These P c -S-K r relationships are highly nonlinear in nature and depend on flow hydrodynamics (dynamic/static) conditions, capillary and viscous forces, contact angles, grain size distribution, surface tension, boundary conditions, fluid properties, and length scales of observation. The conventional relations of the capillary pressure are functions of wetting-phase saturation (P c -S relationships) and differences between the average pressures for oil/nonwetting and the water/wetting, 15,16 which can be mathematically expressed as to the water or wetting-phase saturation. Steady-state P c -S-K r relationships have been previously studied by a large number of authors, some of whom, including [17][18][19][20][21][22] argue that they are insufficient to characterize the behavior of two-phase flows in porous media fully, as they fail to capture the dynamic flow. Hassanizadeh and Gray 17,18 proposed a generalized P c -S relationship with the inclusion of dynamic capillary pressure as shown below P c;dyn À P c;equ 8 : 9 ; s j ¼ Às @S=@t ð Þ s j (2) where s is a dynamic coefficient, qS/qt is the time derivative of saturation, P c,dyn is the dynamic capillary pressure, P c,equ is the steady-state capillary pressure, and the measurements for all variables is taken at the same saturation value. The dynamic coefficient can be determined from the slope of Eq. 2, which exhibits a linear relationship in the form of a straight line. As reported in the literature, Eq. 2 has been subsequently used by many workers. Hassanizadeh et al., 22 Das et al., [23][24][25] and Mirzaei and Das 26 have reviewed these works in significant detail. Das et al. 23 determined the effects of fluid properties on the dynamic coefficient without considering the effects of temperature. This work extends the work by Das et al. 23 as discussed below.
In the context of investigating multiphase flows within porous media and the quantification of dynamic coefficient and its dependence on various physical parameters, a thorough review of literature reveals that a majority of the studies conducted in this area are applicable for two-dimensional domain, whereas studies in three-dimensional cases being reasonably limited. In addition, none of the studies reported above has been attempted to quantify the effects of temperature on the dynamic coefficients directly. However, some authors report the effects of temperatures on the two-phase porous flow behavior. 5,[27][28][29][30][31][32][33][34][35][36][37][38][39][40] Sinnokrot et al. 27 determined the effects of temperature on capillary pressure curves through measurement of drainage and imbibition capillary pressures for three consolidated sandstones and one limestone core. Lo and Mungan 41 proposed that oil recovery operations are more efficient at higher temperature. They show that the oil relative permeability increases as the residual oil saturation is decreased. They carried out experiments at room temperature and 300 F for Berea sandstone cores, concluding that the change in relative permeability curves is caused by viscosity reductions. Nutt 28 showed that the displacement of oil by other fluids is affected by capillary bundle size (pore size volume), pore size distribution, interfacial force, interfacial tension, and viscosity ratio, besides showing that at higher temperature faster oil recovery rates are observed. Davis 29 experimentally investigated the effects of temperatures on the P c -S relationships for water-oil and water-air modes by using a mixture of hydrocarbon oil in two silica sands with different grain sizes, concluding that the temperature changes have a major effect on the residual wetting and nonwetting phase saturation while having a negligible effect on capillary pressure. Davis 29 argued that in water-oil systems, the residual wetting saturation increases as the temperature increases, whereas the residual nonwetting saturation decreases as the temperature increases. Davis 29 also reported that an increase in temperature also increases the nonwetting to wetting relative permeability ratio, thereby enhancing the fluids movement by reducing the displacement pressures, which is a major advantage in oil recovery and water remediation techniques.
Grant and Salehzadeh 31 carried out calculations of several models that incorporate the temperature effects on wetting coefficients and capillary pressure functions. They also con-sidered temperature effects on liquid-gas phase interfacial tension and liquid-solid interfacial tension (surface tension). They found that linear relations of temperature and liquid-gas interfacial tension fit well with the reference interfacial tension data from Haar et al. 42 They noted that the capillary pressure function sensitivity on temperature is mostly due to capillarity. Gatmiri and Delage 32 later numerically studied coupled thermodynamic flow behavior of saturated porous media using the finite element method. A mathematical model was developed to deal with thermal variations in saturated porous media. A new concept called thermal void ratio state surface was introduced to include the thermal effects and the stress state level influence on volume changes. Variation of water permeability, water and solid unit weight due to thermal effects, and pore pressure changes were included. They concluded that the variation of permeability tensor in their equation incorporates the effects of temperature changes on dynamic viscosity, thus affecting the capillary pressure-saturation relationships. She and Sleep 33 studied the temperature dependence of capillary pressure-saturation relationships for water-perchloroethylene (PCE) and water-air systems in silica sand with temperatures ranging from 20 to 80 C. They reported that for water-PCE system the irreducible water saturation increases and the residual nonwetting saturation decreases as the temperature increases. This result was quite similar to the results of Sinnokrot et al. 27 and Davis 29 even though they used different fluid pair (water-CHEVRON 15 white oil) and different sand domain (consolidated natural sandstone). In terms of the capillary pressure, She and Sleep 33 concluded that it decreases when the temperature increases. However, the authors showed that there are other effects besides interfacial tension and contact angles that can have a significant role in the temperature dependence of capillary pressure-saturation relationships. Grant and Salehzadeh 31 modified Laplace's equation to determine the capillary pressure functions. Then, by relating Parker and Lenhard 43 hysteretic P c -S w relationships with equations from Grant and Salehzadeh, 31 temperature dependence of the P c -S w relationships was established.
Narasimhan and Lage 34 investigated temperature dependency on viscosity and hence on the global pressure drop of the flowing fluid in porous media. They showed that with an increase in temperature the viscosity and the global pressure drop across the domain decrease. Muralidhar and Sheorey 35 investigated the displacement of oil and the saturation pattern, which lead to viscous fingering under isothermal and nonisothermal water flow conditions. Investigations were carried out for both homogenous and heterogeneous domains with high-and low-permeability values with water injection temperatures of 50 and 100 C. They conclude that there are three forms that can occur in the presence of both viscous and capillary forces; stable displacement, viscous fingering, and capillary fingering and the fact that higher temperature affects the water saturation in the reservoir and reduces the capillary pressure. They show that the oil-water flow has a strong dependency on the permeability distribution.
Grant 36 proposed a modified mathematical model that takes into account the effects of temperature on capillary pressure-saturation relationships in homogeneous and heterogeneous porous media by using one-parameter and two-parameter models. Grant's 31,35,36 models are modified forms of the Van Genuchten 44 equation and one-parameter equation model presented by Grant and Salehzadeh. 31 The model takes into account the interfacial tension as a function of temperature. It was determined that the two-parameter model is better than the one-parameter model for a system with a larger temperature range. Grant verified his extended model with temperatures ranging from 273.16 to 448 K. Hanyga and Jianfei 37 investigated the dependency of thermal effects on immiscible two-phase flow and implicitly stated that the capillary pressure-saturation relationships with thermodynamic consideration in nonisothermal systems are significantly different with the same relationships in isothermal systems. Schembre et al. 40 carried out experiments that predicted the relative permeability of heavy-oil in situ aqueous phase saturation phase profiles during high-temperature imbibition. In their study, they compared experimental data with the simulated results and used two different nonwetting phases. They carried out experiments at 120 and 180 C. Their results showed that the oil saturation and the water relative permeability decrease when the temperatures increase, resulting in a decrease in capillary pressure.
As it can be inferred from the studies reported above, although there are a number of studies that suggest the importance of temperature effects on P c -S relationships, very little information concerning the temperature dependency of the dynamic effects on capillary pressure relationships can be found. The major concern of this work is therefore to address this issue by carrying out numerical simulations of a two-phase flow system (PCE-water) in homogeneous porous domain. The study is important as we need to understand how quickly or slowly a two-phase flow system reaches equilibrium, to determine the range of validity of dynamic capillary pressure theory 17,18 and to obtain values of the dynamic coefficient (s) at different temperatures. These effects are determined in our simulations through the use of constitutive relations for fluid and material properties (interfacial tension, residual saturation, viscosity, density, and relative permeability) as a function of temperature. A large number of simulations have been conducted to demonstrate how the combined effects of dynamic flow (dependent on different pressure boundary conditions) and temperature variations can affect the capillary pressure relationships. Results presented in this article correspond to temperatures ranging between 20 and 80 C for two different materials namely coarse and fine sands.

Mathematical Model
The complete description of the two-phase flow behavior representing the PCE-drainage process can be characterized by a system of hydrodynamic model: constitutive equations for capillary pressure as a function of entry pressures and saturation conditions, etc., and additional models for the physical properties of the fluids as a function of temperature. The governing and the constitutive model equations used in this work have been listed in Tables 1 and 2, respectively.
As evident in the tables, the two-phase flow in porous media is modeled using an extended version of Darcy's law, 45,46 which represents the momentum balance of the phases under consideration described by Eq. 3, where q c is the Darcy's flux, K rc is the relative permeability, l c is the viscosity of related phase, !P c is the driving pressure gradient, and !z is the upward unit vector. It is defined here that the flow is Newtonian and incompressible. The mass balance of the two phase (wetting and nonwetting) is given by the continuity equation (4), where / is the porosity, q c is the fluid density, S c is the saturation for corresponding phase, q c is the fluid flux, c : w,nw refer to the wetting and nonwetting phases, respectively. ; þ r q c q c 8 : Brooks-Corey function S ew ¼ S w À S wi 1 À S wi 8 > : 9 > ; for 0 S wi 1 (8) Relative permeability of wetting phase Relative permeability of nonwetting phase Liquid-gas interfacial tension Liquid-liquid interfacial tension Water viscosity Oil density Material properties such as porosity, tortuosity, and compressibility parameters are considered to be independent of temperature variations, whereas the change in material particle density can be safely neglected [47][48][49][50] in the temperature range considered. The relation between the wetting-phase saturation and the nonwetting-phase saturation is given by Eq. 5, where S w and S nw correspond to the average saturation of the wetting and nonwetting phases, respectively. The Brooks-Corey functions represented by Eqs. 6-8 are used to determine the capillary pressure-saturation relationships, 51 where S ew is the effective saturation of the wetting phase, S rw is the irreducible saturation of the wetting phase, P d is the displacement or entry pressure, and k is the Brooks and Corey 51 pore size distribution index. The Burdine porosity distribution along with the Brooks-Corey 51 formulation is used to obtain the relative permeability of multiple phases 51 using Eqs. 9 and 10, where K rw and K rnw are the relative permeability for the wetting and nonwetting phases, respectively. The simulator used in this work uses three different porosities namely effective, diffusive, and total. Effective porosity refers to the interconnected pore spaces associated with the gross fluid flow. Diffusive porosity refers to all the interconnected pore spaces, whereas the total porosity refers to both the isolated and connected pore spaces. Isolated pore spaces are assumed to be filled with liquid water. All saturations are defined with respect to the diffusive porosity. The liquid contained within the pore space, represented by the difference between the diffusive and effective porosities, equals the residual water content. This residual moisture content is a function of capillary pressure and consequentially not strictly an immobile fluid.
The density and viscosity of the fluids representing different phases, water residual saturation, and interfacial surface tension (liquid-liquid and liquid-gas) are considered to be temperature variant. Change in water residual saturation is described by Eq. 11, where a r and b r are the constants for water residual saturation function and can be referred to in the work of She and Sleep. 33 In this work, Eq. 11 in conjunction with the use of Brooks-Corey 51 relations, i.e. Eqs. 6-8, is used to account for the effects of temperature on residual saturation. The interfacial surface tension for a liquid-gas system can be mathematically modeled by Eq. 12 given by Grant,36 where c lg is the water-air interfacial tension and a lg 0 , a lg 1 , and c lg are constants for liquid and gas interfacial tension. The interfacial surface tension between the two immiscible liquid phases can be modeled using Eq. 13 given by She and Sleep 33 and Grant and Bachman, 52 where c lo is the water-oil interfacial tension, a and b are constants for water-oil interfacial tension. Water viscosity as a function of temperature is given by Eq. 14 described in the works of Yaws 53 and Oostrom et al. 54 Water density as a function of temperature is modeled using steam tables [55][56][57] represented in the form of Eq. 15, where a 1 …, a 12 , A 11 …, A 22 , are constants for the water density function, q l is water density, P w r is the reduced water pressure, T w r is reduced water temperature, and v w c is the critical specific water volume. The oil viscosity as a function of temperature variation is given by Eqs. 16 described in the works of Reid et al. 58 Time derivative of saturation @S @t  Quasi-static and dynamic simulations using computational fluid dynamics-based two-phase flow simulator STOMP have been carried out within coarse and fine sand domains after which the averaging procedures have been used to evaluate the capillary pressure curves and the corresponding dynamic coefficients for a range of temperatures between 20 and 80 C.
Simulator used, grid refinement, and convergence The two-phase flow behavior has been simulated in this work by using the ''water-oil'' mode of the simulator ''Subsurface Transport Over Multiple Phases'' (STOMP) (http:// stomp.pnl.gov/, White and Oostrom, 59 Nichols et al. 60 ), which is written in FORTRAN. The water-oil (W-O) mode of the simulator is a fully implicit, integrated finite difference code, which solves the model equations for the water and oil components simultaneously, was used to simulate aqueous and NAPL flow through the porous media in many previous studies. The code has been previously used and validated for a range of engineering problems. [61][62][63][64] The governing equations of the two-phase porous flow model are discretized using the finite volume technique, whereas the temporal discretization procedure follows the implicit-theta time-stepping approach. Simulation results were obtained on a computational grid, which was refined by increasing the total number of cells by 25%, and no major differences were observed. Convergence was achieved by altering the maximum Newton-Raphson iteration between 16 and 50 and by applying a tolerance value of 10 À6 . Results indicate that grid refinement has negligible effect on the final capillary pressure-saturation curves. A maximum time step of 600 s and time acceleration factor of 1.125 have been chosen to yield stable solutions.
Quasi-static and dynamic simulation procedure 1. Quasi-static simulations are carried out by initializing the PCE pressure to zero value. Subsequently, the PCE pressure at the top of domain is gradually increased until the capillary pressure (P c ) equates to the porous domain entry pressure (P d ) when PCE starts to infiltrate the porous domain. Simulations are carried out until saturation at each grid point does not change with time or the variation is smaller than the tolerance limit (10 À6 ). The simulation is continued by increasing the PCE pressure at the top for the attained P c -S curves. Simulation is stopped when the residual saturation is attained at P c equal to 11,000 Pa.
2. Dynamic simulations of the two-phase flow are carried out by imposing a high PCE pressure at the top of the domain for every simulation until the residual saturation is attained (refer Das et al. 23,65 and Mirzaei and Das 26 for details).
3. To carry out the postprocessing of the data, an inhouse FORTRAN program that performs volume averaging to determine saturation weighted average of capillary pressure is used. For determining the dynamic coefficient knowledge of saturation time, derivative data and the capillary pressure data under dynamic and steady-state conditions are processed from the simulation output. Averaged data for (P c,dyn À P c,equ ) and (qS/qt) are then fitted on to the straight line by carrying out linear regression analysis. If the R 2 values are in the vicinity of 1 (R 2 ! 0.95), we define that the data fit well to a straight line. The final capillary pressuresaturation curves are then generated using Microsoft Excel. As the set of output data is huge, Tecplot Macro files that can relatively reduce the times for postprocessing are generated.

Averaging procedure to compute the dynamic coefficient
The dynamic coefficient (s) can be determined by using Eq. 2, which in turn requires the estimation of average capillary pressure for the domain. The averaging procedures that have been used in this work for the calculation of (s) can be found in detail in the works of Das et al. [23][24][25] and Mirzaei and Das. 26 Basically, the average capillary pressure and water saturation are calculated after each time step for both steady-state and dynamic capillary pressure curves. The equations used for averaging of the field variables for generating the capillary pressure curves have been listed in Table 3. The averaging calculations have been implemented using a numerical integration module of Tecplot 360 and an in-house FORTRAN program.
The average values represent the saturation-weighted average of fluid pressures in individual nodes within the domain. Average capillary pressure is obtained using Eqs. 18 in conjunction with the use of saturation-volume relation for the wetting and nonwetting phases represented by Eqs. 19, where t n (s) is an arbitrary nth time step, \P nw [|t n and \P w [|t n (N m À2 ) are the volume-averaged nonwetting-and  wetting-phase pressures, P nwj and P wj (N m À2 ) are the nonwetting-and wetting-phase pressures, V wj and V nwj are the volume of the wetting and nonwetting phases on the jth node of volume V j , and S nwj and S wj (-) are the saturation corresponding to nonwetting-or wetting-phase pressure at time t n (s) in an arbitrary jth node, where j ¼ 1, 2, 3,..., m, with m being the total number of nodes in the computational domain. Average water saturation is computed by averaging the saturation values obtained on the node through use of Eq. 20. Finally, the time derivative of saturation (qS/qt) is evaluated using the average water saturation values at different time levels in conjunction with the use of Eq. 21, which represents a central-differencing scheme, where S w |t nþ1 is the average wetting-phase saturation at time step (n þ 1).

Description of Simulation Test Cases
Simulations carried out in this work imitate the pressure cells used for investigating the drainage process through generation of P c -S relationships. The details of the 3-D domains, material properties, fluid physical parameters, and initial and boundary conditions used in this work for simulating various test cases at different temperatures are described in briefly below.

Computational domain
A 3-D cylindrical column shown in Figure 1 similar to a previous setup 26 has been used in this work.
The diameter of the column is 10 cm, whereas the vertical length in z-direction is 12 cm. The computational grid consists of (4 Â 4 Â 24) cells with 416 nodes. The nodal spacing(s) for the domain are listed in Table 4. As heterogeneity effects have not been included, the domains have been considered to be fully isotropic.

Material and physical parameters
The material properties used in this work are defined not to vary with temperature and correspond to two different sands (coarse and fine), which are listed in Table 5. It should be noted that the same parameters 23-26 that provide a viable way to benchmark our simulation results at a reference temperature of 20 C before gaining confidence in simulations conducted at additional temperatures.
The various temperature-dependent fluid physical parameters based on the constitutive equations (11)(12)(13)(14)(15)(16)(17) are listed in Table 6. In addition, the constants required for evaluating the fluid physical parameters in Eqs. 11-17 are presented in Table 7. Simulations have been conducted for a range of temperature values between 20 and 80 C, where the effects of temperature variation have also been included on the initial and bottom boundary pressure values. Table 8 lists the simulations that have been carried out at different temperatures for both coarse and fine sands. At any given temperature, four dynamic simulations and one steadystate simulation have been carried out with a view to generate numerically stable and reliable capillary pressure curves.

Initial and boundary conditions
Capillary pressure-saturation relationships and the dynamic coefficient can be determined by simulating the flow through and the pressure cell experiments for drainage process. 26,63,66,67 The sand samples are assumed to be fully saturated with water and placed inside a reservoir filled with a nonwetting phase, PCE in this case which is a common DNAPL contaminant. At the top face of the porous domain, we impose different PCE pressures for different dynamic cases listed in Table 9.
Constant pressure boundary conditions for PCE flow and no-flow conditions for water at the top of the domain are applied. On the sides of the domain, no-flow condition for both phases has been imposed. At the bottom boundary, a constant pressure for water and no-flow condition for PCE is applied. The outflow is assumed to be in level with the top of the domain. Simulations carried out are for the vertical displacement of the phases (direction aligned with gravitational force) and therefore the effect of gravity is included.  However, no attempts have been made in this work to specifically investigate the effects of gravity on the dynamic coefficients.

Results and Discussion
Results obtained from our simulations of the drainage of water by PCE in coarse and fine sands within 3-D domains at different temperatures have been presented in the form of P c -S curves for four dynamic and one steady-state flow. Thereafter, we discuss the effects of temperature on dynamic coefficients for both the coarse and fine sands. 3-D contours of saturation have also been presented for both the domains to show the effects of temperature on the times required for the flow to reach equilibrium. The results presented in this study are therefore expected to reduce the uncertainty associated with the quantification of flow relaxation parameters under the effects of varying temperature.
Dynamic and steady-state capillary pressure curves Figures 2a, 3a, 4a, and 5a illustrate the dynamic and steady-state P c -S curves for coarse sand at 20, 40, 60, and 80 C, respectively, whereas Figures 2b, 3b, 4b, and 5b illustrate the dynamic and steady-state P c -S curves for fine sands at 20, 40, 60, and 80 C, respectively. Observing P c -S curve in Figure 2a, we see that the initial capillary pressure is 224 Pa. At high saturation between 1 and 0.35, the capillary pressures in dynamic cases (1)(2)(3)(4) gradually rise as the saturation values increases.
However, there is a sudden pressure drop in the steadystate case from 270 to 249 Pa between the saturation values of 0.449 and 0.416. This occurs at the very start of the drainage process between 6 s and 2 h. At saturation value of 0.325 a pressure jump is observed for all the dynamic cases. The jumps correspond to simulation times of 60 s for dynamic case 4 (PCE head 135 cm) and 1200 s for dynamic case 1 (PCE head 50 cm). After this jump, the capillary pressure further increases as the saturation decreases. As discussed by Mirzaei and Das 26 and Das et al., 65 these jumps are known as Haines 68 jumps.
In Figure 2b, the initial capillary pressure is 800 Pa, which is higher than for the coarse sand simulation at the same temperature.
No pressure jumps are observed as in the case of fine sand. In the dynamic and the steady-state cases, stable increase in P c is observed from 1200 s to 2 h. At lower saturation values, from 0.4 to residual saturation, P c increases rapidly, whereas the rate of decrease in saturation values is much slower. This occurs only after 9-10 h. Overall, the increment of P c in each dynamic case is not as drastic as in the case of coarse sand.
Similar trends are observed for the P c -S curves plotted in Figures 3a, 4a, 5a and 3b, 4b, 5b for coarse and find sands at 40, 60, and 80 C, respectively. In Figure 3a for the coarse sand case at 40 C, the initial capillary pressure is 230 Pa for the dynamic and steady-state cases. However, unlike in the 20 C case in Figure 2a, the steady-state P c does not show a rapid drop. The dynamic cases exhibit a pressure jump at  saturation value of 0.39. It can be noticed that the capillary pressure at 40 C is relatively higher than the pressure at 20 C for the same values of saturation. Figure 3b represents the P c -S curves for fine sand at 40 C where the initial capillary pressure is 823 Pa. The pressure for each case starts increasing rapidly from saturation of 0.356 until near the residual saturation. This occurs between 9 and 20 h after the PCE pressure is injected. Again in this case, no pressure jumps are observed. The trend at 40 C is similar to the one observed at 20 C.
In Figure 4a for coarse sand at 60 C, the initial P c is 236 Pa, and a similar pattern for the P c -S relations is observed at 20 and 40 C. The pressure jump for the dynamic case at a saturation value of 0.47 and a constant increase of P c in steady-state case is observed. The capillary pressure value at 60 C is higher than the values at 20 and 40 C for the same saturation.
The initial capillary pressure is 844 Pa in Figure 4b for the fine sand P c -S curves at 60 C. The pressure for each case at 60 C is higher than the corresponding P c values for the fine sand cases at 20 and 40 C for the same saturation values. The P c values for each case increase rapidly at the saturation value of 0.447 from the start of the drainage process, which occurs between 10 and 30 h. Figure 5a represents the P c -S curves for the coarse sand at 80 C. The pressure jump in dynamic cases occurs at saturation value of 0.5. The initial capillary pressure is 246 Pa, and the pressure values in 80 C are higher than the values computed in 20-60 C coarse sand cases at the same saturation values. After the pressure jump, the difference in pressure values between different dynamic cases becomes more significant. Finally, at a physical time 1800 s, the saturation decreases very rapidly.
In Figure 5b for fine sand at 80 C, it is observed that capillary pressure for each case increases significantly after saturation value of 0.449 until it reaches the residual saturation. Moreover, the saturation decreases for all dynamic cases much slower after 0.449. This occurs after around 10-30 h from the start. The qP c values are higher than qP c values at 20, 40, and 60 C.

Effects of temperature on dynamic coefficient
In Table 10, the computed dynamic coefficient values for coarse and fine sands at 20 C have been listed. The relation between dynamic and steady-state capillary pressure differences and the time derivative of determines the dynamic coefficient with the inclusion of an intercept term.
For the coarse sand at 20 C, the values of P c (dynamic)-P c (steady-state) range from 1.86 to 198 Pa corresponding to the saturation values of 0.28 and 0.191, respectively. The slope increases for lower value of saturation. Saturation time derivative ranges from À3.19 Â 10 À5 to 1.28 Â 10 À5 s À1 , whereas the dynamic coefficient ranges from 4.2 Â 10 6 to 9.1 Â 10 6 Pa s for the chosen saturations. The dynamic coefficient increases whilst qS/qt decrease as the saturation decreases until the residual saturation. R 2 values lie between 0.978 and 0.998, which means that at 20 C for the chosen saturation values the relationship of (qS/qt) to dP c is quite linear. For the fine sand case, the values of qP c range from 0.872 to 193 Pa corresponding to saturation values of 0.31 and 0.0714, respectively. Saturation time derivative ranges between À1.27 Â 10 À6 and À7.5 Â 10 À8 s À1 , whereas the dynamic coefficient ranges from 6.4 Â 10 6 to 3.5 Â 10 9 Pa s for the same saturation value. The dynamic coefficient increases whilst qS/qt decreases as the saturation decreases to near the residual saturation. R 2 values vary between 0.62 and 0.98, signifying a nonlinear relationship between qP c and qS/qt.
In general, the calculated dynamic coefficients in the 20 C case are higher for fine sands in comparison to the coarse sand case at same temperature. Similar observations can be made by looking at Tables 11-13 for 40, 60, and 80 C, respectively.
Elevated temperatures result in a higher capillary pressure and a smaller saturation time derivative. Dynamic coefficient at higher temperature is higher than its value at lower temperature for the same water saturation. These trends are more significant after a critical saturation (low saturations near the residual saturation) at each temperature. The R 2 values show that the data in some cases may not be linear. This is more significant in the case of fine sand domain where the displacement process is very slow and the effects of difference in injection pressure are not as significant as in the coarse sand case. The values of the dynamic coefficient in fine sand for different temperature ranging from 6.4 Â 10 6 to 1.1 Â 10 11 Pa s; however, in the coarse sand the dynamic coefficient ranges from 4.2 Â 10 6 to 5.33 Â 10 8 Pa s. Figures 6a, b represent the dynamic coefficient plots as a function of temperature for coarse and fine sand cases, respectively, which corresponds to the data listed in Tables  10-13.
It is noticed that the dynamic coefficients for fine sands are relatively higher than for coarse sands at the same temperature and saturation values. Reviewing Eq. 2 and the results for the dynamic coefficients and the time derivative saturation, we observe that as the dynamic coefficient increases the time required for attaining the residual saturation decreases for each dynamic and steady-state case. This trend is consistent with our observations at different temperatures, where the times required to attain the flow equilibrium have been listed in Tables 14-17.   In terms of different sand types, the time derivative saturation is mostly affected by the relative permeability and the entry pressure. In fine sand, a low value of relative permeability (5 Â 10 À12 m 2 ) causes the infiltrated fluid (PCE) to displace the existing fluid (water) very slowly. In addition, the infiltrate fluid needs higher pressure or larger time at the same pressure to overcome the entry pressure and start infiltrating the pore space.
As expected, our results demonstrate that the relationships between capillary pressure and saturation are nonlinear for water-PCE system in homogeneous porous media, specifically at high saturation values or when (qS/qt) is very high. This nonlinearity is more significant in dynamic case simulations, when the injection PCE pressure is much higher than the entry pressure. In addition, this nonlinearity is also observed near residual saturation, when the capillary pressure increment is high whilst (qS/qt) is very small. In this study, we only report saturation after an hour when the saturation values have already decreased to half. As it can be observed from Tables 10-13, as we find that some R 2 values are not very close to unity, the existing dynamic capillary pressure-saturation relationships might need further modifications. The results obtained at 20 C (for coarse sand) are similar to those by Mirzaei and Das, 26 confirming the applicability of dynamic P c -S for certain ranges of saturation values. Although the nonlinearity decreases when (qS/qt) or the saturation decreases, results show that nonlinearity occurs again near the residual saturation even though with a lesser degree. In terms of the proposed relations for dynamic coefficient, we observe that there is a possibility that the intercept terms may increase with a decrease in saturation or in (qS/qt). However, these relations may not be linear and need further investigation along with the relations at saturation range close to saturated condition.
We also observe that an increase in temperature gives a slight increase in capillary pressure at the same water saturation. This is quite different with some other studies that have been suggested that the capillary pressure decreases with elevated temperature. 31,36 The reason for this difference may be due to the fact that we impose a lower pressure for water at the bottom of the domain for elevated temperatures. This may lead to higher capillary pressures at earlier time steps. The computed saturation time derivatives are smaller at higher temperature; hence, the times needed to attain the residual saturation get longer. If we consider Eqs. 6-10, this trend may occur because we applied an increment in the water residual saturation that may lead to slow reduction of saturation. However, as a fluid flow in porous media involves many factors such as viscous forces, buoyancy, fluid-solid surface tension, and coupled effects from many physical parameters (material and fluid properties), we cannot justify that only one physical phenomenon dominates. The decrease in fluid viscosities and fluid/fluid interfacial surface tension will generally result in higher mobility of the fluids. In the present case, the decrease in fluid densities and through imposing different water residual saturations may result in this effect. Several factors that may be important to characterize the flow movement such as contact angles and wettability are ignored by the STOMP simulator for water-oil mode. However, as we impose a higher water residual saturation for higher temperature, the time taken to reach the water residual saturation decreases.

Effects of material properties on dynamic coefficient
Material properties in a porous domain are correlated with one another. Consequently, the material properties such as porosity, intrinsic relative permeability, density, and pore size distribution tend to have a lumped effect on the dynamic coefficient. The boundary pressures for same temperature range are same for both fine and coarse sands. The fine sand has lower relative permeability and higher entry pressure than in the case of the coarse sand. However, we impose the same water residual saturation and similar porosity for both fine and coarse sand at a given temperature.
Higher capillary pressure is observed in the fine sand case for the same water saturation. The dynamic coefficient values for the same water saturation are higher in fine sand than in coarse sand, whereas the saturation time derivative is high in fine sand. These results are similar to the experimental results obtained by Mirzaei and Das. 26 The water saturation needs more time to reach the new saturation in fine sand. This slower movement of fluids and slower displacement process is observed in fine sand.  High capillary pressure and lower Darcy's velocities result as the infiltrate fluid (PCE) needs more time to reach the entry pressure. Once PCE attains the entry pressure it begins to infiltrate. A much lower permeability in fine sand also acts as a resistance to the fluids movement whilst the capillary pressure is very high. This results in high capillary pressures and longer times for PCE to displace water in the fine sand test cases, which are tabulated in Tables 14-17.
In addition, we have observed pressure jumps in coarse sand cases (Figures 2a, 3a, 4a, and 5a) at some particular saturation value when the simulation period is generally less than 30 s. These are characterized as Haines 65,68 jumps, which occurs because of the presence of viscous and/or gravity forces along with the low entry pore space pressure and higher intrinsic permeability. This combination results in a jump when the PCE pressure outcomes the entry pressure of pore space and starts to displace the entrapped water inside each pore space. This explains why in the case of fine sand where the intrinsic relative permeability is much lower and the entry pore space pressure is much higher, the Haines 65,68 jumps do not occur.

3-D contours of aqueous phase saturation
The 3-D saturation plots presented in Figures 7 and 8 correspond to the changes occurring in the aqueous-phase saturation for coarse and fine sands, respectively, at 20 and 80 C.
These plots correspond to the Tables 14 and 17, which illustrate the times required to attain the residual saturation at different temperature values for coarse and fine sands under dynamic conditions.
As it can be observed in the 3-D contour plots of saturation (Figures 7 and 8), the effects of temperature mostly occur at the earlier time steps, i.e., at 30 s, 60 s until 1 h, when sudden displacements are observed. Thereafter, the displacement gets slower, which is marked by a slower decrease in water saturation. At 20 C and at 3 and 60 s, PCE displaces a lot of water at the bottom of the porous domain. When the temperature is increased, water amount that is displaced in the first 30 and 60 s decreases. At 80 C, water amount that is displaced in the first 30 and 60 s is reduced. Gravitational effects and viscous forces along with low water pressure may play an important role slowing the displacement as the temperature increases. The water displacement is not distributed evenly at each node   in the cylindrical domain. Initially, the water content in the bottom part is displaced, and the displacement proceeds upward with water in the sides of domain being displaced earlier than the water in the middle part of the domain.

Conclusions
Simulations have been carried out in this work to quantify the effects of temperature on the dynamic capillary pressure and subsequently dynamic effects on two-phase porous flow in cylindrical domain. This is important to understand how quickly and slowly a two-phase flow system reaches equilib-rium, to determine the range of validity of dynamic capillary pressure theory 17,18 and to obtain values of a dynamic coefficient (s) at different temperatures. Series of simulations have been carried out within the range of 20-80 C for two different materials (coarse and fine sands) under dynamic conditions to determine the values of the dynamic coefficient. The dynamic coefficients values for coarse sand ranged between 4.2 Â 10 6 and 5.33 Â 10 8 Pa s, whereas for the fine sand case the range was between 6.4 Â 10 6 and 1.1 Â 10 11 Pa s. We observe that with an increase in temperature the dynamic coefficients increase, and at any given saturation value large time scales are required for nonaqueous phase flow (NAPL) displacements.  The dynamic coefficients are also found to be nonlinear functions of saturation and temperature. This is contrary to the results by Stauffer, 69 which suggest that the dynamic coefficient is independent of saturation. However, our results are consistent with recent studies. 22,26,70,71 The temperature effects as reported in this work have never been reported before, as far as the authors know.