Modeling of relative permeabilities including dynamic wettability transition zones
MModeling of relative permeabilities including dynamicwettability transition zones
Abay M. Kassa , Sarah E. Gasda K. Kumar F. A. Radu Department of Mathematics, University of Bergen, P. O. Box 7800, 5020 Bergen, Norway. NORCE, Nyg˚ardsgaten 112, 5008 Bergen, Norway.Corresponding author: Abay M. Kassa (E-mail: [email protected])
Abstract
Wettability is a pore-scale property that impacts the relative movement and distribution offluids in a porous medium. There are reservoir fluids that provoke the surface within pores toundergo a wettability change. This wettability change, in turn, alters the dynamics of relativepermeabilities at the Darcy scale. Thus, modeling the impact of wettability change in the relativepermeabilities is essential to understand fluids interaction in porous media. In this study, weinclude time-dependent wettability change into the relative permeability–saturation relation bymodifying the existing relative permeability function. To do so, we assume the wettability changeis represented by the sorption-based model that is exposure time and chemistry dependent. Thispore-scale model is then coupled with a triangular bundle-of-tubes model to simulate exposuretime-dependent relative permeabilities data. The simulated data is used to characterize andquantify the wettability dynamics in the relative permeability–saturation curves. This studyfurther shows the importance of accurate prediction of the relative permeability in a dynamicallyaltering porous medium.
Wettability alteration (WA) plays an important role in many industrial applications such as mi-crofluidics nanoprinting, enhanced oil recovery (EOR), and CO storage [1, 2, 3, 4, 5]. Wettabilityrefers to the tendency of one fluid over the others to spread on or adhere to a solid surface [35, 1] andis defined by the fluid-fluid contact angle (CA). This pore-scale property regulates the distributionof fluids in the pore spaces and controls the relative flow of immiscible fluids in a porous medium[38, 37, 35, 1]. This, in turn, impacts constitutive relations in the multi-phase flow systems such asresidual saturation, relative permeability, and capillary pressure at the Darcy scale [36, 39, 40, 4, 35].Investigating and upscaling the impact of WA on the constitutive relations (hereafter constitutiverelations refers to the relative permeability- and capillary pressure-saturation relations) is of greatimportance.Wettability is assumed to be static in time and uniform in space. However, wettability is a dy-namic process that depends on surface chemistry, composition of fluids, exposure time, and reservoirconditions (pressure and temperature) to name a few [41, 42, 38, 11, 10, 9, 43]. Experiments oncrude oil/brine/rock systems have shown that adsorption of active components from the crude oilis able to change the wettability of the sample porous medium from water-wet to intermediate-wetsystem [41, 44, 9]. It is also hypothesized that the oil reservoir may be more oil-wet than what is1 a r X i v : . [ phy s i c s . g e o - ph ] A ug bserved from the experiment. This is because the adsorption time of the experiment period wasmuch less than the age of the oil in the reservoir. Furthermore, CO is one of the reservoir fluidswhich contain active components that can provoke the surface within the pores to undergo a WA[45, 46, 47, 48, 49, 50, 51, 52, 53, 4].Generally, the WA process can have three phases that delineate the transition from initial to finalwetting-state conditions, e.g. initial-wet, final-wet, and dynamic-wet. The end (initial and final)wetting conditions are static in time but can be uniform and mixed in space. A mixed-wet conditioncould be created by rock mineral and exposure history differences. This is due to the fact that a poresurface exposed to the WA agent may be altered to a new wetting condition, while the unexposedsurface keeps the initial wetting state [54, 7]. This creates a mixed-wet condition even within asingle pore and was observed and explained first by Salathiel et al. [41] in the 1970s. Usually, WAis assumed to occur instantaneously and is considered as a function of the WA agent concentration.In some cases, however, the alteration process might take prolonged time in the scale of weeks andmonths [14, 18, 9, 55]. In this regard, the dynamic-wet phase can be a function of exposure time inaddition to the WA agent concentration.The WA process may result in a saturation function alteration for subsequent drainage-imbibitiondisplacements and thus cause hysteresis in constitutive relations [56, 36, 57, 58, 59]. For instance,core-flooding measurements for (supercritical or gas) CO -water system have indicated that WA-induced alteration in the residual saturation and capillary pressure curves despite the fact thatthey were measured following a standard procedure, i.e., where “pressure equilibration” is obtainedafter each increment in pressure [17, 19, 14, 20, 45, 13]. In these measurements, a steadily changein capillary pressure function over time was observed. More importantly, the capillary pressuredeviation from the initial-wet state curve could not be explained by classical scaling arguments. Theinstability and gradual change of residual saturation and capillarity through exposure time, in turn,impact the behavior of relative permeabilities.The above experiments reveal that standard constitutive models are not well suited to predictrelative permeability and capillary pressure under dynamic, long-term WA. One alternative is touse mixed-wet model, e.g. Kjosavik et al. [60] and Lomeland et al.[61], that capture the staticheterogeneity of wettability in the relative permeabilities. The main feature of these models istheir flexibility to describe hysteresis and scanning curves caused by a wettability gradient in space.Other alternatives are models designed to handle the instantaneous WA process in the relativepermeabilities. The first class of these models involves a heuristic approach that interpolates betweenthe initial and final wetting states in which the WA effect is captured as a coefficient function[25, 3, 27, 28, 62]. Interpolation models are conceptually simple, while the initial and final wettingstates are characterized by standard functions, e.g. Brooks-Corey [63] or van Genuchten [64]). Theother approach incorporates the effect of the instantaneous WA into the relative permeabilitiesthrough the residual saturation directly [26]. To date, only Al-Mutairi et al. [29] have consideredthe effect of time-dependent WA in both the relative permeabilities and capillary pressure functionsexplicitly. However, their model does not sufficiently incorporate or upscale the WA processes tocore-scale laws.Appropriate upscaling of the pore-scale time-dependent WA process connected to the capillarypressure function was the subject of our recent work [65]. There, WA dynamics were upscaledby introducing a mechanistic time-dependent CA model at the pore-level that was coupled with acylindrical bundle-of-tubes model and used to simulate capillary pressure curves for drainage andimbibition displacements. The simulated data was used to formulate and quantify a interpolation-based capillary pressure model at the Darcy scale. The new dynamic model resolves the existinginterpolation models used in the studies of reservoir simulation [28, 25, 3, 27, 62] by includingthe dynamics in time and quantifying the pore-scale WA process to the interpolation model in a2ystematic manner.One may consider employing a similar approach to [65] and an interpolation-type model tocapture the pore-scale underpinnings of WA in the relative permeability behaviors. However, time-dependent WA may impact the capillary pressure and relative permeabilities in different ways. Asobserved in [65], WA has a direct impact on the entry pressure in each pore and reflects it atthe Darcy scale. Furthermore, a small change in CA exerts a large impact on the dynamics ofthe capillary pressure function. However, the relative permeability alteration occurs when the WAaffects the pore filling/draining orders of pore-sizes. This may lead to a longer exposure time toobserve a relative permeability deviation from the initial-wet state curve. Furthermore, unlike thecapillary function, the relative permeability curves are constrained between zero and one for anychange of wettability. These features of the relative permeability may impact the modeling approachto upscale the pore-scale WA process to the relative permeability behavior.To our knowledge, a physically reliable model to characterize a prolonged exposure time-dependentWA induced dynamics in the relative permeability behaviors has not been proposed yet. This paperrevises and extends the approach discussed in [65] to develop a relative permeability model thatincludes pore-level time-dependent WA processes. Section 2 summarizes two possible approachesthat can be applied to upscale the impact of time-dependent WA in the relative permeabilities. Thefluid-fluid CA change is designed as a function of exposure time to the WA agent at the pore level tomeasure the WA process. This model is coupled with a pore-scale, triangular bundle-of-tubes, modelto simulate time-dependent WA induced relative permeability curves. These curves are presented inSection 3 and are used to evaluate the modeling approaches hypothesized in Section 2. A time-dependent WA may introduce a dynamic term in the relative permeability–saturation ( k rα - S )relationship. This dynamic term can be measured by its deviation from the static initial wetting-state as: k rα ( · ) − k irα ( S α ) := f d α ( · ) , (1)or, the dynamic term can be correlated with the parameters of the standard models k rα ( · ) = k irα ( S α , a ( · ) , b ( · ) , . . . ) , (2)where a ( · ) and b ( · ) are fitting parameters that change along exposure time, whereas f d α representsthe WA induced dynamic component, and the subscript α ∈ { w, n } represents the wetting andnon-wetting phases.In this study, we explore both approaches in Eqs. (1) and (2) to quantify and characterize f d α inthe relative permeabilities for a system that undergoes a WA. From Equation (1), we propose aninterpolation model following our previous work [65], where the dynamic component is designed tointerpolate between two end wetting-state curves. To obtain an interpolation model, the dynamiccomponent in Eq. (1) can be scaled by the difference between the initial and final wetting-state rel-ative permeability curves. The resulting quantity is non-dimensional and referred to as the dynamiccoefficient , ω α (cid:0) k frα − k irα (cid:1) = f dα , (3)where the superscript i and f represents relative permeabilities at the initial and final wetting-statesrespectively. This can be substituted into Eq. (1) to obtain dynamic relative permeability models k rα = (1 − ω α ) k irα + ω α k frα , (4)3here ω α , the dynamic coefficient, is responsible for capturing the wettability dynamics at themacroscale. Similar model to Eq. (4) were employed to include the impact of instantaneous WA intothe relative permeability curves [25, 3, 27, 28, 62].The approach in Eq. (2) relies on a systematic inclusion of the dynamic term f d α into the rel-ative permeability function through the model parameters. This can be done by formulating theparameters a ( · ) and b ( · ) as a function of exposure time and WA agent in similar fashion as ω α . Thisapproach is motivated by the fact that the parameters in the standard relative permeability modelsare adjusted to different values when wettability changes from one state to the other.Both the initial and final wetting-state curves can be characterized fully by the well-knownrelative permeability models such as van Genuchten[64] or Brooks-Corey [63], Purcell [66], the LETmodel [61] or a model proposed by Kjosavik et al. [60]. For the sake of brevity, we focus on theBrooks-Corey (BC) and LET models in this study. The BC relative permeabilities can be derivedby integrating the capillary pressure over the capillary tubes [34, 60]. After the integration of theBC capillary pressure, one can obtain relative permeabilities k irw = S a w w , and k irn = (1 − S a n w )(1 − S w ) m n , (5)for the water-wet system and k frw = (1 − S a w w )(1 − S n ) m w , and k frn = S a n n , (6)for a hydrophobic system (see [60]), here a w , a n , m w , and m n are data fitting parameters in which thesubscripts w and n indicate the wetting and non-wetting surfaces. Particularly the m ’s are known tobe tortuosity exponents. In 2005, Lomeland et al. [61] have proposed relative permeability modelswith three parameters L, E, and T for two-phase flow system. Their correlation models read as: k irw = S L w w S L w w + E w (1 − S w ) T w , and k irn = (1 − S w ) L n (1 − S w ) L n + E n S T n w , (7)where L α , E α , and T α are data fitting empirical parameters. A detailed description and explanationof the parameters can be found in [61]. The important characteristics of the LET model (7) is itsflexibility to predict the relative permeability curves for any type of wettability conditions. We notethat the parameters, a and b , in Eq. (2) are associated with data fitting parameters in Eqs. (5)-(6)and (7).The WA induced dynamics in k rα - S relation should be characterized in a unified manner inorder to evaluate the behaviors of ω α ( · ), a ( · ), and b ( · ) along exposure time to the WA agent. TheWA induced k rα - S data can be measured from laboratory experiments. However, this approachis expensive in terms of time. Thus, we follow a theoretical approach in which we simulate time-dependent k rα - S data from a pore-scale model. We employed a triangular bundle-of-tubes to represent the pore-scale model. Note that one can alsouse a simpler pore-scale model, i.e., a cylindrical bundle-of-tubes model, to characterize the impactof WA on k rα - S relations. However, polygonal pores are advanced in the way that they can representphysical processes such as the establishment of mixed wettability within a single pore. This maylead to different fluid distributions within a pore, and establishment of non-wetting fluid layers inthe corners of the pore space and drainage through layers [54, 67, 72].A bundle-of-tubes model is a collection of capillary tubes with a distribution of radii as depictedin Fig 1. The tubes in Fig 1 are connected with the wetting (right with pressure P res.r ) and non-4igure 1: The fluid displacement scenario in a bundle of tubes that is connected with wetting andnon-wetting phase reservoirs. The right column shows fluid distribution during primary drainage.For complete (drainage-imbibition cycles) fluid configurations, see Fig 2.wetting (left with pressure P res.l ) phase reservoirs. Once the fluid movement is initiated in thetubes, the fluid configurations in each tube can have the form as in Fig 1. Here, A b,m , A c,m , and dx represent the bulk area covered by the non-wetting fluid, the corner area still covered by the wettingphase, and change of fluids along the tube respectively, whereas α is half angle and θ m is fluid-fluidCA. Detailed calculation of these areas and the angle β m is discussed below.Let the boundary pressures difference be defined as∆ P = P resl − P resr , (8)and the tubes in the bundle are filled with the wetting phase initially. To displace the wetting phasefluid in the m th tube, the pressure drop has to exceed the local entry pressure [22]∆ P > P c,m . (9)If condition (9) is satisfied, the non-wetting fluid starts to displace the wetting phase, and thevolumetric flow rate in the m th can be approximated by the Lucas-Washburn flow model [68], q m = G m ( R m , θ m )[∆ P − P c,m ]8[ µ nw x intm + µ w ( L − x intm )] , (10)where, µ nw and µ w are non-wetting and wetting fluid viscosities, respectively, the superscript int stands for fluid-fluid interface, q m = dx intm /dt is the interface velocity, R m is the m th tube inscribedradius, and θ m is the fluid-fluid contact angle at pore m . Here, θ m is a general contact anglerepresentation and can be specifically defined as θ r,m if the interface is receding, θ a,m if the interfaceis advancing, and θ h,m if the interface hinges in the corner of the pores. The G m in Eq. (10) representsthe conductance of the fluids and we follow the work of [67] to pre-compute the conductance. Theinterface is assumed to be trapped when it reaches the outlet of the tube, thus q m = 0 when theinterface reaches the boundaries.The wettability change and/or gradient in the polygonal pores may create distinct fluid configu-rations, for example see Fig 2. These configurations can be encountered during drainage (e.g. A, D,and E) and imbibition (e.g. B, C, and F) displacements in each pore in the bundle. The bold surfacein configurations B through F is to show that part of the surface is exposed to a WA agent at somepoint and experiences a wettability change. Configuration B occurs when the wettability is alteredup to a CA value satisfying θ < π − α , where α = π is the corner half-angle and configuration Cotherwise. Configuration F may occur when the non-wetting fluid invades configuration E. This isthe case when the receding CA is sufficiently smaller than the previous advancing contact angle.5igure 2: Fluid configurations for primary drainage, imbibition, and secondary drainage (water inlight color and CO in a dark color). We adopted the configurations from [72]. The bold lines alongthe sides indicate altered wettability.Below, relevant aspects of the polygonal pores for the relative permeabilities measurement willbe stressed, which includes entry pressure and the areas covered by the two fluids. The capillarypressure at the pore-level is given by P c,m = p n − p w = σr m , (11)where p n and p w are the non-wetting and wetting phase pressures respectively, r m is radius of arcmeniscus (AMs) that separates the bulk from the corner fluid, and σ is fluid-fluid interfacial tension.The radius of curvature r m is determined from the minimization of Helmholtz free energy of thesystem [69]. For isothermal, constant total volume, constant chemical potentials, and incompressiblesystem, the minimization of the change in Helmholtz free energy can be simplified to [69, 43, 70] P c,m dV n = σ ( dA nw + cos( θ m ) dA ns ) , (12)where θ m is fluid-fluid surface angle, dV n is the change of non-wetting fluid volume, dA nw , and dA ns represent the change in area of fluid-fluid and fluid-solid interfaces respectively. In [72, 71, 74] theentry pressure curvature, r m , and the fluid layer existence criterion were calculated from Eqs. (11)and (12) for each of the configurations in Fig 2. Here, we follow and implement the work of [72] tocalculate the radius of curvature r m and the associated entry pressure.Once we obtained r m for the m th tube, fluid volumes (the bulk and corner fluids) can be calculatedin terms of r m and θ m . To calculate these quantities, we numbered the fluid-fluid interfaces in orderfrom the apex if there exist more than one interface in the corner, and we apply the indicatornotation I k = (cid:26) , if interface k separates bulk nonwetting and corner water − , if interface k separates bulk water and corner nonwetting . (13)6he bulk cross-sectional area A kb,m in each tube in the bundle is defined as, A kb,m = R m tan α − r m b km sin( β km + α ) + 3 r m β km , if I k = 1 , R m tan α − r m b km sin( β km − α ) − r m β km , if I k = − , (14)where b km = r sin( β km )sin( α ) , and β km = (cid:40) π − α − θ km if I k = 1 π + α − θ km if I k = − . (15)Though we used a general notation for CA θ km above, θ km may be replaced by, θ r,m and θ a,m whenthe interface recedes and advances respectively, and θ kh,m if the interface, separating the bulk andcorner fluids, is hinging. The hinging contact angle (if exists) changes with the entry pressure P c,m according to: θ kh,m = arccos( P c,m b km sin( α ) σ ) − α if I k = 1arccos( P c,m b km sin( α ) σ ) + α if I k = − . (16)The fluid area that occupies the corner regions can be estimated by A c,m ( θ m ) = 3 r m (cid:16) θ m + α − π θ m ) (cid:16) cos( θ m )tan( α ) − sin( θ m ) (cid:17)(cid:17) , (17)where θ m is an argument to determine the appropriate area. The corner surface covered by waterin configuration C is obtained by A c,m ( θ h,m ), whereas the non-wetting fluid layer in configurationE is calculated from A c,m ( π − θ a,m ) − A c,m ( θ h,m ). The same approach can be applied if there existmany layers in the corner of the pore.As clearly seen above the areas, A b,m and A c,m are dependent on wettability, i.e., fluid-fluid CA.For example, fluid configuration C and D occurs only if the condition θ m ≤ π − α is satisfied duringthe drainage displacement. Otherwise, the non-wetting phase occupies the cross sectional area ofthe tube, and the entry pressure calculation is reduced to the well known Young-Laplace equation.On the other hand, the non-wetting fluid layer occurs in the corner when the condition θ m > π + α is satisfied. This implies that a dynamic change of CA can also determine the fluid distribution ina single pore. Above, we observed that the entry pressure in Eq. (12), fluid distributions in Eqs. (14)-(17), and aconductance in Eq. (10) are wettability evolution dependent. In this section, we introduce a WAmechanism at the pore level to examine its impact on the saturation distribution. We recall thatthere are many factors that provoke the surface within the pore to undergo a wettability change.Here, we consider the effect of exposure time and fluid-history on the CA change with the assumptionthat: • The CA of the pore surface is altered through exposure time to the WA agent and the alterationis permanent. That means the wettability is not restored to the original wetting condition whenthe WA agent is displaced by the other fluid unless the displacing fluid has a composition thatgradually restores the original wetting condition.7
The WA becomes quasi-static in time if the WA agent is removed from the pore before the finalwetting-state is reached. If the agent is reintroduced at some later point, alteration continuesuntil the final state.According to the assumptions above, the bulk surface area A b,m is supposed to be altered dynamicallyin time, whereas the corner surface area A c,m may keep the initial condition. Thus, we introduce afunctional form to describe a WA mechanism for any arbitrary tube m as θ m ( · ) := (cid:40) θ im for A c,m ,θ im + ϕ ( · )∆Θ for A b,m , (18)where ∆Θ = θ fm − θ im , θ fm , θ im are the final and initial contact angles respectively. The WA model(18) is designed to evolve from an arbitrary initial wetting state to the final wetting condition sothat ϕ is used to interpolate between end wetting conditions and has a value between zero and one.Theoretical investigation and detailed laboratory measurement on time-dependent WA is verylimited. Furthermore, WA is a complex process, where surface free energies, surface mineralogy, fluidcomposition and exposure time interact. However, adsorption of the WA agent onto the surface areais a natural process in CA change. Such adsorption type wettability evolution is observed for CAmeasurements [48, 50, 49, 75, 77, 76]. These all give an insight to model ϕ in Eq. (18) according tothe adsorption of the WA agent and can be given as follows ϕ := χ m C + χ m , (19)where C is a non-dimensional parameter that controls the speed and extent of alteration from initial-wet to final-wet system. The derivation of Eq. (19) can be found in our previous work [65]. Thevariable χ m is a measure of exposure time and is defined as χ m := 1 T (cid:90) t A b,m x intm V p,m dτ, (20)where T is a pre-specified characteristic time, x intm is the fluid-fluid interface position along thetube length L , A b,m is the pore surface area that covered by the non-wetting fluid, and V p,m is thepore volume. Without losses of generality, the characteristic time T is set to be the time for onecomplete drainage displacement under static initial wetting condition, which can be pre-computedfrom Eq. (10).For a given interface position x intm , one can determine the required time to reach the specifiedinterface position from Equation (10). The obtained time is used in Eq. (20) to calculate the exposurehistory of the m th pore to the WA agent. According to Eq. (20), each individual pore surface area A b,m will experience CA change based on the exposure time to the altering fluid, which may bedifferent for different pores depending on the local saturation history. This process would givesrise to a non-uniform wetting condition across the bundle until all pores have reached the finalwetting-state. Moreover, the bulk surface area A b,m is the only surface that undergoes WA, and thecorner surface area that covered by water is not subject to WA. This results in mixed-wet condition.However, this paper do not consider the wettability gradient across the length of the tubes becausethe time to drain is assumed to be fast compared to the exposure time for WA to occur. The wettability dynamics described in Eq. (18) are coupled into a triangular bundle-of-tubes modelto simulate relative permeability curves according to Algorithm 1. Here, the relative permeability8or phase α is calculated as k rα = Q α µ α L K A T ∆ P α , (21)where Q α is the volumetric total flow rate of phase α , K is absolute permeability of the bundle, A T is the cross-sectional area of the bundle. We have repeated algorithm 1 for a few numbers of Algorithm 1
A single drainage-imbibition cycle. Fluid and rock properties are given according toTable 2
Drainage displacement
Set P max c % the maximum capillary pressure while ∆ P < P max c do Calculate P c from Eq. (12) if ∆ P > P c then Calculate S nw and χ = T (cid:82) t S nw dτ Calculate K rα Update θ m from Eqs. (18) and (20) end if Increase ∆ P end while Imbibition displacement
Set minimum entry pressure P min c while ∆ P > P min c do Calculate P c from Eq. (12) if ∆ P < P c then Calculate S nw and χ Calculate K rα Update θ m from Eqs. (18) and (20) end if Increase ∆ P end while drainage-imbibition cycles provided that the WA process is completed within these displacements.The flow rate (or the pressure drop ∆ P ) is controlled in an arbitrary manner in order to gain θ f foreach tube within a few numbers of drainage-imbibition cycles. When we reduce the increment ofeach ∆ P , the flow rate becomes very slow. This imposes a prolonged exposure to the WA agent andthus χ m grows and eventually results in a large change in CA. For example, each ∆ P increment isreduced by three order of magnitude for the last cycle compared to the first cycles to complete thealteration process.The relative permeabilities–saturation “data points” are obtained in each drainage-imbibitioncycle, and the obtained data is presented in the following Section. The generated k rα - S curves areused to quantify the dynamics in the relative permeabilities which caused by time-dependent WA.The goal is to develop a correlation model that involves only a few parameters. Finally, the relationbetween these parameters and changes in the pore-scale WA model parameter C is studied andexamined. 9 Simulation results
The two-phase flow simulation tool at the pore-scale (Algorithm 1) is implemented in MATLAB.The pore-scale model consists of parallel triangular tubes that connect the non-wetting and wettingreservoirs. Each tube in the bundle is assigned a different radius R, with the radii drawn from atruncated two-parameter Weibull distribution [67] R = ( R max − R min ) (cid:110) − δ ln (cid:104) x (1 − exp( − /δ )) + exp( − /δ ) (cid:105)(cid:111) /γ + R min , (22)where R max and R min are the pore radii of the largest and smallest pore sizes respectively, and δ and γ are dimensionless parameters. The rock parameters and fluid properties are listed in Table1. These parameters are coupled to the bundle-of-tubes model to simulate fluid conductance andparameters values unit parameters values unit σ R min µ m R max µ m θ fm
180 degree θ im µ w µ nw δ γ In section 2, we point-out that end wettingstate relative permeabilities are the foundation to charac-terize the dynamic relative permeability curves. Thus, it is natural to examine the end-state k rα - S relations before quantifying the dynamic relative permeabilities using the same pore-size distribu-tion and fluid properties in Table 1. Thus, we simulated static k rα - S data by fixing the wettabilityat pre-specified initial θ im and final θ fm values in each tube, see Table 1. The simulated curves areplotted along the saturation path in Fig 3.We correlated both the BC (in Eqs. (5)-(6)) and LET (in Eq. (7)) models with the simulated(initial and final wetting-state) static k rα - S curves, and the result is compared in Fig 3. The fittedparameters for the correlation models are found in Table 2. Both the BC correlations and LETmodels give an excellent match to the simulated phase relative permeability curves under staticconditions.As a complement, we have done experiments on the sensitivity of model parameters (not shownhere) for different pore-size distributions by setting the wetting condition to be static. In thisexperiment, we found that the parameters L n and T w in the LET model can be unity for any pore-size distribution. Furthermore, from the end wetting-state correlation results, we observe that E α is the only parameter that depends on wettability change, whereas the BC model parameters aresensitive for wettability change, see Table 2. More importantly, we have noticed that the parameters L w and T n have the same value and thus, we can represent them with a single parameter (say λ ).Furthermore, the parameter E n is the inverse of E w . In this regard, we can reduce the LET model10odel Parameters Initial value ( θ i ) Final value ( θ f ) a w m w a n m n L w E w T w L n E n T n S w K r w Initial state dataBC modelFinal state dataLET model (a) S w K r n Initial state dataBC modelFinal state dataLET model (b)
Figure 3: The BC and LET correlation models compared to the initial and final wetting-state relativepermeability curves: (a) wetting phase and (b) non-wetting phase relative permeabilities.to a two-parameter (but in each phase) model and we call it reduced LET model which read as k rw = E n S λw (cid:0) E n S λw + 1 − S w (cid:1) − , and k rn = (1 − S w ) (cid:0) − S w + E n S λw (cid:1) − . (23)The model in Eq. (23) is a two parameter model with only one parameter that varies along thewettability change. This makes the reduced LET model more reliable than its counter part, i.e., theBC model. Note that the parameter λ can also vary with wettability. In this case, a better matchwith the simulated relative permeability data can be obtained. We demonstrate five drainage-imbibition cycles and simulate dynamic k rα - S relations in each cycle(see Fig 4), while the tubes in the bundle are altered through time (see Fig 5) following the CAmodel (18). These data are generated with a pore-scale parameter C = 10 × − . Note that theCA change may be halted (temporarily) in the pores if the displacement is to configuration D afterimbibition. In this case, the drainage curve may follow the previous imbibition path. However, the11mbibition curve may show a deviation from the previous drainage curve if wettability is alteredsufficiently. S w K r w Initial stateSimulated dataFinal state 0.10.20.30.40.50.60.70.80.9Time (a) S w K r n Initial stateSimulated dataFinal state 0.10.20.30.40.50.60.70.80.9Time (b)
Figure 4: Simulated dynamic relative permeability curves ((a) wetting phase and (b) non-wettingphase) with respect to wetting phase saturation. The static k rα - S curves for the initial and finalwetting states are shown as a reference. The color code shows the k rα - S dynamics within a year ofexposure time.According to the CA distribution in Fig 5 and k rα - S curves in Fig 4, the wettability that rangesfrom strongly to weakly water-wet does not affect the pore filling and draining orders of the pore-sizes. As a consequence, the k rα - S relations for drainage/imbibition displacements follow the samesaturation path of the initial-wet condition. This implies that a WA induced k rα - S hysteresis maynot occur in a straight bundle-of-tubes model for a wettability range of 0 ◦ ≤ θ m ≤ ◦ . Similar k rα - S relations are reported for a pore-network model with a wettability range of θ m = 0 ◦ to 45 ◦ [36]. Thisis in contrast to the capillary pressure–saturation relation, where a small change of CA impacts the P c - S path significantly [65]. However, imbibition/drainage displacement may not necessarily occurin monotonically increasing/decreasing order of pore-sizes when the wettability of the pores (some)are altered to intermediate/weakly hydrophobic. This results in a relative permeability–saturationpath deviation as observed in Fig 4.We observe that the wetting and non-wetting phase relative permeabilities steadily increase anddecrease (see Fig 4) respectively in each subsequent drainage-imbibition displacements when thewettability evolves from hydrophilic to hydrophobic conditions. This occurs because the relativelylarger pores start allowing the water (originally wetting phase) to imbibe through before the smallerpores when the medium is altered to intermediate/hydrophobic system. On the contrary, the orig-inally non-wetting fluid prefers the smaller pores to flow in during drainage. As a consequence,a mobility reduction and improvement respectively for the non-wetting and wetting phase relativepermeabilities are observed. However, any additional drainage-imbibition cycle would follow alongthe static curve for the final wetting state once the final CA ( θ f ) is reached.Analogous to the capillary pressure–saturation relation [65], time-dependent WA introduces dy-namic hysteresis in the k rα - S relations for a bundle-of-tubes model as shown in Fig 4. The k rα - S hysteresis imposes a non-unique relation between the relative permeabilities and saturation. Thisis one of the challenging features of WA during the quantification of the dynamics in relative per-meabilities. To eliminate the hysteresis observed in the k rα - S relation, we projected the simulatedrelative permeability data onto the temporal domain χ in Fig 7. The temporal domain, χ is the12 Radius [ m] C on t a c t ang l e [ deg r ee s ] Figure 5: Dynamic CA evolution as a function of exposure time to the WA agent per each tube.This CA distribution was recorded at the end of each drainage-imbibition cycle and the color codeshows the CA dynamics during a year of exposure time.measure of the exposure history in averaged sense which defined as: χ = 1 T (cid:90) t S nw dτ. (24)Unlike k rα - S , k rα - χ is uniquely related but non-monotonically i.e., it raises to one and decreases tozero in time along each drainage-imbibition cycle.The WA process also affects the corner fluid distribution in each drainage/imbibition displace-ment. Fig 7 shows the evolution of (average) corner wetting and non-wetting saturations. Thewetting phase saturation decreases through time while the layer saturation grows. This is becausewettability is altered from water-wet to hydrophobic and thus, the originally non-wetting fluid prefersto be in the corners. However, when (all) pores become more hydrophobic, the corner water increasesand bulges-out during imbibition. This process may result in a layer collapse and has shown in thefifth imbibition, where the layer saturation increases during larger pores were imbibed and decreasewhen smaller pores imbibed over exposure time. Here, we have checked the establishment of cornerfluids during the calculation of k rα - S relations in each drainage-imbibition displacements. The interpolation approach, similar to Eq. (4), was applied successfully to capture time-dependentWA mechanisms in the capillary pressure curves [65]. Thus, it is important to test the potential ofthe interpolation-based model (in Eq. (4)) to predict time-dependent dynamics in k rα - S relation.The dynamic coefficient ω α is calculated according to Eq. (2) and plotted in Fig 8. In Fig 8b, weobserve that the dynamic coefficient ω α is related to the exposure time χ non-monotonically. Thus,it is challenging to propose a functional relationship between ω α and χ directly. On the other hand,the dynamic coefficient ω α in Fig 8a is increasing with respect to the exposure time to the WA agent.13 .2 0.4 0.60.10.20.30.40.50.60.70.80.9 k r w [ - ] (a) k r o [ - ] (b) Figure 6: The relative permeability data: (a) wetting phase relative permeability and (b) non-wetting phase relative permeability as a function of χ for the non-uniform WA case. The color ofeach data point indicates the time elapsed in years. S w (a) S w (b) Figure 7: Corner fluid saturations: (a) wetting phases saturation and (b) non-wetting phase satu-ration The color of each data point indicates the time elapsed in years.However, the ω α - S curves in Fig 8a have piece-wise functional forms i.e., zero and Langmuir-typefunction of phase saturation, that are altered with exposure time. This imposes another challengeto find a smooth model to predict the curves along with the saturation. But, one can design apiece-wise function to correlate the ω α − S data. From Fig 8a, we observe that the starting pointof the Langmuir functional form is exposure time-dependent along the saturation path. Thus, the ω α − S can be represented as, ω α ( S w , χ ) := (cid:40) S w − S ∗ w ( χ ) S w − S ∗ w ( χ )+ β ( χ ) , for S w > S ∗ w ( χ ) , , otherwise (25)14 .2 0.4 0.6 0.8 S w (a) (b) Figure 8: The scaled dynamic deviation from the initial wetting-state relative permeability curvesas a function of wetting phase saturation (a) and exposure time χ (b).Parameter Vale Parameter Value a a b -2.522 b -2.086 c c S ∗ w and β are time-dependent. The variable S ∗ w is used to transform the starting point ofthe Langmuir part ω α along the saturation to zero, and β is used to determine the curvature of thecurve. From Fig 8a, we observe that the parameters S ∗ w and β are decreasing functions of exposuretime, χ .We matched the designed model in Eq. (25) with the dynamic coefficient data, ω α − S w , todescribe the relations S ∗ w − χ and β − χ . The obtained relations have the form, S ∗ w = a χ b + c , and β = a χ b + c , (26)where a i , b i and c i for i = 1 , k rα = (cid:40) S w − S ∗ w ( χ ) S w − S ∗ w ( χ )+ β ( χ ) (cid:0) k frα − k irα (cid:1) + k irα , for S w > S ∗ w ( χ ) k irα , otherwise . (27)The dynamic model (27) is compared with the simulated relative permeability in Fig 9. Accordingto Fig 9, the designed model predicts beyond the initial wetting state at joint-point of the initialwetting-state model and the designed interpolation model. This shows that the designed model (27)is badly correlated with the simulated k rα - S data. Furthermore, the non-smoothness behavior of ω α − S w results in a piece-wise phase relative permeability model with many dynamic parametersto be calibrated. 15 .2 0.4 0.6 0.8 S w K r w Simulated dataInterpolation model 0.10.20.30.40.50.60.70.80.9Time (a) S w K r n Simulated dataInterpolation model 0.10.20.30.40.50.60.70.80.9Time (b)
Figure 9: Comparison of the interpolation model (27) with the simulated (wetting (a) and non-wetting (b)) relative permeabilities.The second approach discussed in Section 2 relies on establishing a relation between the modelparameters and χ to upscale the effect of pore-scale wettability evolution in relative permeabilities.This is supported by the result reported in Table 2 for end wetting-state curves in which the modelparameters are dependent on the wetting condition of the porous domain in addition to the pore-sizedistribution. As a consequence, we can formulate dynamic correlation models from standard k rα - S relations, i.e., BC models in Eqs. (5)-(6) or reduced LET models in Eq. (23). So far, we noticed thatonly E α is sensitive for wettability change in the reduced LET model, whereas both parameters aresensitive in the BC model (see Table 2). This implies that the BC model is more expensive thanthe reduced LET model to upscale the effect of WA on relative permeabilities. Thus, we choose thereduced LET model to represent the dynamic relative permeability curves. Here, we rearrange thereduced LET model as k rw = L n ( χ, E n ) S λw L n ( χ, E w ) S λw + 1 − S w , and k rn = 1 − S w − S w + L n ( χ, E n ) S λw , (28)where λ and E n are determined from the initial wetting-state correlation. Here, L n is designed tochange with CA. Thus, we coupled the reduced LET model (28) with curve fitting tool in MATLABand matched with the k rα - S data to study the functional dependencies between the L n and χ . Theobtained relation is linear and has the form L n ( χ, E n ) = a n χ + E n , (29)where E n is as given in Table 2 and a n is dynamic fitting parameter for wetting and non-wettingphase relative permeabilities and determines the slope of the relative permeabilities along exposuretime. For this particular simulation this parameter is estimated to be a n = 3 .
57 for all dynamicdrainage-imbibition cycles reported above. Now the dynamic term L n in Eq. (29) can be substitutedinto the reduced LET model (28) to give the dynamic relative permeabilities: k rw = ( a n χ + E n ) S λw − S w + ( a n χ + E n ) S λw , and k rn = 1 − S w − S w + ( a n χ + E n ) S λw , (30)where the parameters E n and λ are pre-determined from the initial wetting-state correlation. Theproposed dynamic relative permeability models in Eq. (30) and the simulated k rα - S data are com-16ared in Fig 10. From Fig 10, we observe that the proposed dynamic model (30) correlated well S w K r w Simulated dataModefied LET model 0.10.20.30.40.50.60.70.80.9Time (a) S w K r n Simulated dataModefied LET model 0.10.20.30.40.50.60.70.80.9Time (b)
Figure 10: Comparison of the modified LET model and the simulated phase relative permeability:(a) wetting phase and (b) non-wetting phase.with the simulated relative permeability curves. The proposed relative permeability model is single-valued regardless of the number of drainage-imbibition cycles. However, this single-valued model isnot well predictive around the junction points, particularly for low wetting saturations. However,the obtained correlation result shown in Fig 10 is acceptable and a significant improvement on theinterpolation model result shown in Fig 9. Furthermore, the modified LET model prediction can beimproved by letting the parameter λ to vary along the exposure time. However, this may (at least)double the number of parameters in the model that need to be calibrated.If we compare the two dynamic models (i.e., the piece-wise interpolation model in Eq. (27) andthe reduced/modified LET model in Eq. (30)), the reduced LET model is more efficient and simplerto implement in the Darcy flow than the piece-wise interpolation model. Because, the reduced LETmodel is smooth along time and saturation. We also note that the number of parameters in themodified relative permeability model in Eq. (30) is reduced by half from the original LET model.These all make the reduced dynamic model more reliable than using a model consisting of multipleparameters that change in each cycle (or hysteresis models). Thus the analysis below will concernonly on the modified LET model. In this section, we will investigate the response of the upscaled model parameter a n to the changein the CA model parameter C in Eq. (18). The parameter C controls the extent of WA at thepore-level, whereas a n determines the WA induced dynamics in the relative permeabilities. We havesimulated different drainage-imbibition k rα - S curves by varying C to draw a relation between a n and C . To do so, we determine the parameter a n from each k rα - S data that was simulated byconsidering different values of C . Then, we correlate the estimated parameter values of a n with thechosen values of C , which can be read as a n = P C + P , (31)where P and P are fitting parameters. The correlation result is plotted in Fig 11, where theparameters are estimated to be P = − P = 6. Therefore, we can predict the upscaled17ynamics of the relative permeabilities directly from the pore-scale WA process. -5 a n [ - ] DataModel
Figure 11: The relation between pore-scale wettability parameter C and the correlation parameter a n in Eq. (30).The relation in Eq. (31) can be substituted into the dynamic relative permeability model (30) tocomplete the upscaling process. The resulting dynamic relative permeability models are saturation,exposure time, and pore-scale WA parameter dependent. The pore-scale parameter has to be esti-mated from the calibration of CA model (18) with experimental data. Validating the underlying CAchange model is beyond the scope of this paper. Rather, we consider the CA model as a reasonablebasis to perform and analyze the upscaling process. The saturation path that used to generate the relative permeability data in Fig 12 can be consideredas one of the arbitrary paths in the domain S w × χ . This saturation path is dependent on thehistory of the exposure time to the WA agent (degree of WA change) and the reversal-point ofsaturation. This implies that the relative permeability-saturation dynamics may behave differentlyif one chooses a different saturation path that entails a prolonged exposure time for a fixed saturationprofile and/or flow reversal at intermediate saturation.Here, we simulate as many as possible k rα - S curves that involve different reversal-points andexposure history within the S w × χ . Note that we use the pore-scale model parameter C = 10 × − .The simulated data is plotted in Fig 12a and 12b for phase relative permeabilities. These arbitrarycurves are used to test the potential of the modified LET model in Eq. (30). To do so, we apply thecalibrated dynamic relative permeability model (30) to generate the k rα - S - χ surface. The absolutedifference between the simulated data and the surface generated by the calibrated model is depictedin Fig 12c. According to the results in Fig 12, we can justify that a single saturation history pathis sufficient to calibrate a dynamic model that can be applied to any saturation-time path. In contrast to the dynamic capillary pressure model presented in [65], the interpolation-based ap-proach is poorly correlated with the simulated relative permeability curves. However, further in-18 a) (b)(c)
Figure 12: Top: Simulated relative permeabilities obtained by taking multiple paths in the S w × χ space, for the wetting (left) and non-wetting (right) phases. Bottom: The difference between thedynamic model with the simulated datavestigation may be needed to re-evaluate the potential of capturing the WA process in the relativepermeabilities based on the interpolation approach. We also examined the BC model to upscale theWA induced dynamics in the relative permeabilities though we do not report it here fully. For thesake of brevity, we have highlighted only the response of the BC model to the wettability change inSection 3.3 for end wetting conditions. In general, other models that involve more than two param-eters (sensitive to CA change) could be calibrated with reasonable accuracy. But, these parametersneed to be adjusted in each drainage-imbibition displacements. This complicates the modeling pro-cess and the resulting model may involve many parameters that may impose an extra challenge toanalyze the WA impact on the flow dynamics at the Darcy scale.In this study, we have developed a dynamic relative permeability model by modifying an existingstatic-wettability model i.e., the LET model. The original LET model consists of three parametersin each phase that need to be adjusted differently for different wetting conditions. First, we didnumerical experiments to study the dependency of these parameters on the pore-size distribution,while the wettability was kept constant. We found that T w for wetting and L n for non-wetting phaserelative permeability models are constant for any pore-size distribution. From this, we reduced theLET model to a model (the reduced model in Eq. (23)) that involves two parameters in each phase.19e further investigate the sensitivity of the reduced LET model parameters to the WA dynamics,and we found that only E α is dependent on CA change dynamics. We then draw a clear relationbetween E α and the exposure time χ which allows us to come up with a single-valued dynamic relativepermeability model that represents any arbitrary drainage-imbibition cycle. We note that the modelinvolves two types of parameters. The first one is pore-size distribution dependent parameters E n , L w , and T n which are determined a priori from the initial wetting-state correlation. Knowing thesevalues, the parameter a n is the only parameter that controls the WA induced dynamics in the relativepermeabilities for both phases.The proposed model is at the Darcy scale, that allows for a change in relative permeability as afunction of averaged variables such as saturation ( S w ) and exposure time to a WA agent ( χ ). Thisimplies that the relative permeability in a grid block that is exposed to the WA agent may changeover time even for constant saturation profile. However, if the grid block is not exposed to the WAagent, χ is zero for this particular grid block. In this case, the dynamic model predicts the initialwetting-state curve. However, the developed model may continue further after the end wetting-statecurve has attained which is in contrast to the interpolation-type model, where prolonged exposuredoes not contribute to the dynamics once the final wetting curve is met, see [65]. Thus it is importantto propose a strategy to ensure that the relative permeability dynamics do not to cross the end-statecurve. Above, the WA induced dynamics in the relative permeabilities is represented by L n ( χ ) = a n χ + E n (32)where E n is known from the initial wetting-state correlation. The final wetting state is attainedwhen a n χ + E n = E fn is satisfied. From this, we can estimate the exposure time needed to reachthe final wetting state (say χ max ) such that the relative permeability is represented by the endwetting-state curve. After knowing this we can set the dynamic variable as χ = (cid:40) T (cid:82) to S nw dτ, if χ < χ max χ max , if χ ≥ χ max (33)This controls the unnecessary dynamics once the final wetting-state curve is predicted. Nevertheless,the dynamic term in the model pushes the relative permeability towards the higher and lower endof the curve for the wetting and non-wetting phases respectively.Previous studies [25, 3, 27, 28] represent the impact of instantaneous WA on relative permeabili-ties by an interpolation model which matched directly to core-scale data in a heuristic manner. Thisstudy revealed that the interpolation model is not the best approach to upscale the pore-scale WAprocess. Rather, we have shown the potential of a modified LET model to capture the underlyingWA process at the pore-scale represented by a triangular bundle-of-tubes. The proposed modelsare smooth and simple to use for practical applications. Most importantly, the models are designedto eliminate the relative permeabilities hysteresis induced by CA change during drainage and im-bibition displacements. Similar to the developments in [65], we have quantified the link betweenthe pore-scale model parameter C and the core-scale parameter a n . We estimated the core-scale(dynamic) parameter a n by varying the pore-scale parameter C . According to the simulation re-sults, we have shown that a very simple scaling can relate a pore-scale process with the core-scale.This result implies that knowing the mechanism that determines the CA change at the pore-levelcan be used to predict the macroscale dynamics without performing pore-scale simulations. This isan important and valuable generalization for making use of experimental data to inform core-scalerelative permeability-saturation relations. 20 Conclusion
In this paper, we developed a dynamic relative permeability model that includes the pore-scaleunderpinnings of WA in the relative permeability–saturation relationships at the Darcy scale. Wefound that the developed model (i.e., the modified LET model in Eq (28)) is simple to use and canpredict WA induced changes in the relative permeabilities. The modified LET model shows a goodagreement with the simulated relative permeability data. Furthermore, this model is independent ofthe saturation-time paths generated by any drainage-imbibition cycles. More importantly, the WAdynamics in the relative permeabilities is controlled by a single-valued parameter that has a clearrelationship with the time-dependent CA change model parameter.
Acknowledgement
Funding for this study was through the CHI project (n. 255510) granted through the CLIMITprogram of the Research Council of Norway.
References [1]
Bonn, D. & Eggers, J. & Hindekeu, J. & Meunier, J. & Rolley, E.
Rev. Mod. Phys. , 739-805.[2] Iglauer, S. & Rahman, T., & Sarmadivaleh, M. & Al-Hinai, A. & Fernø, M. A. &Lebedev, M.
SPE J. ,1916–1929.[3] auYu, L. & Kleppe, H. & Kaarstad, T. & Skjæveland, S. M. Netw. Heterog. Media , 149–183.[4] Iglauer, S. & Pentland, C. H. & Busch, A. wettability of seal and reservoirrocks and the implications for carbon geo-sequestration. Water Resour. Res. , 729–774.[5] Blunt, M. J. , 197–207.[6] Ahmed, A. & Patzek, T. W.
Water Resour. Res. , 1-12.[7] Blunt, M. J.
SPE J. , 494–510.[8] Morrow, N. R. & Lim, H. T. & Ward, J. S.
Effect of crude-oil-induced wettabilitychanges on oil recovery . SPE J. , 89-103[9]
Buckley, J. S. & Liu, Y. & Monsterleet, S.
SPE J. , 54-61.[10] Jadhunandan, P. P. & Morrow, N. R. 1995 Effect of Wettability on Waterflood Recovery forCrude-Oil/Brine/Rock Systems. SPE Reservoir Engineering , 40–46.2111] Haagh, M. E. J. & Siretanu, I. & Duits, M. H. G. & Mugele, F.
Langmuir , 3349–3357.[12] Singh, R. & Mohanty, K.
SPE J. , 1126-1139.[13] Kim, Y. & Wan, J. & Kneafsey, T. J. & Tokunaga, T. K. and Brine: Pore-scale studies in micromodels. Environ.Sci. Technol. , 4228–4235.[14] Tokunaga, T. K. & Wan, J. capacity. Rev. Mineral. Geochem. , 481-503.[15] Chiquet, P. & Broseta, D. & Thibeau, S.
Geofluids , 112–122.[16]
Chalbaud, C. & Robin, M. & Lombard, J. & Martin, F. & Egermann, P. & Bertin,H. storage. Adv Water Resour , 98–109.[17] Plug, W. J. & Bruining, J. -water system undervarious pressure conditions. Application to CO sequestration. Adv Water Resour Wang, S. & Tokunaga, T. K. and brine in limestone/dolomite sands: Implications for geologic carbon sequestration incarbonate reservoirs. Environ. Sci. Technol. , 7208–7217.[19] Wang, S. & Tokunaga, T. K. & Wan, J. & Dong, W. & Kim, Y. trapping. Water Resour. Res. , 6671–6690.[20]
Tokunaga, T. K. & Wan, J. & Jung, J. & Kim, T. W. & Kim, Y. & Dong, W. and brine in sand: High-pressure P c ( S w ) controller/meter measurements and capillary scaling predictions. Water Resour. Res. , ,4566–4579.[21] Hassanizadeh, S. & Celia, M. & Dahle, H.
Vadose Zone J , 38–57.[22] Dahle, H. K. & Celia, M. A. & Hasanizadeh, S. M.
Transport porousmed , 5–22.[23] Krumpfer, J. W. & McCarthy, T. J.
Faraday Discuss. , 103–111.[24]
Eral, H. B. & ’t Mannetje, D. J. C. M. & Oh, J. M.
Colloid Polym Sci , 247–260.[25]
Delshad, M. & Najafabadi, N. F. & Anderson, G. A. & Pope, G. A. & Sepehrnoori,K.
SPEJ. , 361-370. 2226] Lashgari, H. R. & Xu, Y. & Sepehrnoori, K.
SPE , 1–17.[27]
Andersen, P. Ø. & Evje, S. & Kleppe, H. & Skjæveland, S. M.
SPE J. , 1261–1275.[28] Adibhatia, B. & Sun, X. & Mohanty, K.
SPE , 1-15.[29]
Al-Mutairi, S. M. & Abu-Khamsin, S. A. & Hossain, M. E. Flooding Process.
SPE , 1-12.[30]
Bartley, J. T. & Ruth, D. W.
Transport porous med , 161–187.[31] Helland, J. O. & Skjæveland, S. M.
SPE , , 171–180.[32] Dong, M. & Dullien, F. A. L. & Dai, L. & Li, D.
Trends Anal. Chem. , 1–18.[33] Skjæveland, S. M. & Siqveland, L. M. & Kjosavik, A. & Thomas, W. L. H. &Virnovsky, G. A.
SPE , 60–67.[34]
Xu, W. S. & Luo, P. Y. & Sun, L. & Lin, N.
Plos One , 1–9.[35] Falode, O. & Manuel, E.
Journal of Petroleum Engineering
Ahmed, A. & Patzek, T. W.
Water Resources Research Bobek, J. E. & Mattax, C.C. & Denekas, M.O.
SPE
Anderson, W.
Journal of Petroleum Technology Pentland, C. H. & El-Maghraby, R. & Iglauer, S. & Blunt, M.J. .2011Measurements of the capillary trapping of super-critical carbon dioxide in Berea sand-stone.
Geophys. Res. Lett. Iglauer, S. & Paluszny, A. & Pentland, C. H. & Blunt, M. J. imaged with Xray microtomography. Geophys. Res. Lett. Salathiel, R. A.
J.Petrol. Technol. Treiber, L. E. & Archer, D. L. & Owens, W. W.
SPE J. Morrow, N. R.
Ind.Eng. Chem. Anderson, W. G.
J. Pet. Tech. Wang, S. & Edwards, I. M. & Clarens, A. F. -brine-mineral interface: Implications for geologic carbon sequestration. Environ Sci Technol Bikkina, P. K.
Int. J. Greenhouse Gas Control Yang, D. D. & Gu, Y. G. & Tontiwachwuthikul, P. at high pressures and elevatedtemperatures. Energy Fuels Dickson, J. L. & Gupta, G. & Horozov, T. S. & Binks, B. P. & Johnston, K. P. /water/glass interface. Langmuir Iglauer, S. & Mathew, M. & Bresme, F. interfacial tensions and brine-CO -quartz contact angles and their effects on structural andresidual trapping mechanisms in carbon geosequestration. J. Colloid Interface Sci. bvol386 405–414.[50]
Jung, J. W. & Wan, J. and ionic strength effects on wettability ofsilica surfaces: Equilibrium contact angle measurements. Energy Fuels Espinoza, D. N. & Santamarina, J. C. geological storage. Water Resour Res. Farokhpoor, R. B. & Bjørkvik, J. A. & Lindeberg, E. & Torsæter, O.
Int. J. Greenhouse Gas Control Saraji, S. & Goual, L. & Piri, M. & Plancher, H. /Water/Quartz Systems: Simultaneous Measurement of Contact Angle and Interfacial Ten-sion at Reservoir Conditions. Langmuir
Kovscek, A. R. & Wong, H. & Radke, C. J.
AIChE J. , 1072–1085.[55] Powers, S. E. & Anckner, W. H. & Seacord, T. F.
J. Environ. Eng.
Vives, M. & Chang, Y. & Mohanty, K.
SPE J. . 260–267.[57] Delshad, M. & Lenhard, R.J. & Oostrom, M. & Pope, G.A.
SPE Reser-voir Evaluation and Engineering . 328–334. 2458] Spiteri, E. J. & Juanes, R. & Blunt, M. J. & Orr, F. M. Landry, C. J. & Karpyn, Z. T. & Ayala, O.
WaterResour. Res. Kjosavik, A. & Ringen, J. K. & Skjæveland, S. M.
SPE J Lomeland, F. & Ebeltoft, E. & Hammervold, T. W.
Society of Core Analysts .[62]
Sedaghat, M. H. & Azizmohammadi,S.
Comput Geosci
Brooks, R. H. & Corey, A. T. van Genuchten, M. T. Soil Sci Soc Am J Kassa, A. M. & Gasda, S. E. & Kumar, K. Radu, F. A.
Adv Water Resour .[66]
Li, K. & Horne, R. N.
Water Resour. Res. Hui, M. & Blunt, M. J.
J. Phys. Chem. B , 3833–3845.[68]
Washburn, E.
Phys. Rev Helland, J. O. & Skjæveland, S. M.
WATER RESOURCESRESEARCH . Bradford, S. A. & Leij, F. J.
J. Contam. Hydrol. van Dijke, M. I. J. & Sorbie, K. S. Journal of Colloid and Interface Science , Helland, J. O. & Skjæveland, S. M.
SPE J.
Ma, S. & Mason,G. & Morrow, N. R.
Colloids and Surfaces
Jafari, M. & Jung, J. -water con-ditions: Implication on geological carbon dioxide sequestration. Geochem. Geophys. Geosyst. Morton III, S. A. & Keffer, D. J. & Counce, R. M. & DePaoli, D. W. & Hu M. Z.C.
Journal of Colloid and Interface Science
Davis, A.N. & Morton III, S. A. & Counce, R.M. & DePaoli, D.W. & Hu, M.Z.-C.
Colloids Surf. A221