Impact of Time-Dependent Wettability Alteration on the Dynamics of Capillary Pressure

TL;DR

This paper introduces a dynamic capillary-pressure model that accounts for time-dependent wettability changes at the pore scale, significantly affecting capillarity in porous media over months to years.

Contribution

It develops a novel model incorporating wettability alteration dynamics into capillary pressure calculations, improving understanding of long-term fluid behavior in porous media.

Findings

Wettability alteration can reduce capillarity in CO2 storage.

Dynamic wettability significantly impacts capillary pressure over months to years.

The model enhances prediction accuracy for wettability-altering fluid systems.

Abstract

Wettability is a pore-scale property that has an important impact on capillarity, residual trapping, and hysteresis in porous media systems. In many applications, the wettability of the rock surface is assumed to be constant in time and uniform in space. However, many fluids are capable of altering the wettability of rock surfaces permanently and dynamically in time. Experiments have shown wettability alteration can significantly decrease capillarity in CO$_2$ storage applications. For these systems, the standard capillary-pressure model that assumes static wettability is insufficient to describe the physics. In this paper, we develop a new dynamic capillary-pressure model that takes into account changes in wettability at the pore-level by adding a dynamic term to the standard capillary pressure function. We simulate the dynamic system using a bundle-of-tubes (BoT) approach, where a mechanistic model for time-dependent contact angle change is introduced at the pore scale. The resulting capillary pressure curves are then used to quantify the dynamic component of the capillary pressure function. This study shows the importance of time-dependent wettability for determining capillary pressure over timescales of months to years. The impact of wettability has implications for experimental methodology as well as macroscale simulation of wettability-altering fluids.

Tables

Table 1. Table 1: Parameters used to simulate quasi-static fluid displacement in BoT.

parameters	values	unit	parameters	values	unit
${\sigma }_{ow}$	0.0072	N/m	no. radii	500	[-]
$R_{\min }$	10	$\mu$m	$R_{\max }$	40	$\mu$m
${\theta }_f$	80	degree	${\theta }_i$	0.0	degree
${\mu }_w$	0.0015	Pa.s	${\mu }_{\mathrm{nw}}$	0.0015	Pa.s
$r_{\mathrm{av}}$	23	$\mu$m	L	0.001	m
$\eta$	1.5	[-]

Table 2. Table 2: Estimated correlation parameter values for initial and final wetting state capillary pressure curves

Initial wetting state			Final wetting state
param.	value	unit	param.	value	unit
$c_w$	360	[Pa]	$c_w$	56	[Pa]
$a_w$	0.2778	[-]	$a_w$	0.2778	[-]
$R^2$	1	-	$R^2$	1	-

Equations

$$P_c(\cdot )-P_c^{\mathrm{st},\mathrm{i}}:=f^{\mathrm{dyn}}(\cdot ),$$

$$P_c^{\mathrm{st}}=c_w{\left(\frac{S_w-S_{wc}}{1-S_{wc}}\right)}^{-a_w},$$

$$\omega\big{(}P_{c}^{\rm st,f}-P_{c}^{\rm st,i}\big{)}=f^{\rm dyn},$$

$$P_c=(1-\omega )P_c^{\mathrm{st},\mathrm{i}}+\omega P_c^{\mathrm{st},\mathrm{f}},$$

$$\Delta P=P_l^{res}-P_r^{res}.$$

$$\Delta P>P_c^{\mathrm{mic}}(R_m,{\theta }_m),$$

$$P_c^{\mathrm{mic}}(R_m,{\theta }_m)=\frac{2{\sigma }_{ow}\cos ({\theta }_m)}{R_m},m=1,2,\cdots ,N,$$

$$q_m=\frac{R_m^2(\Delta P-P_c^{\mathrm{mic}}(R_m,{\theta }_m))}{8({\mu }_{\mathrm{nw}}{x}_m^i+{\mu }_{\mathrm{w}}(L-x_m^i))},$$

$${\theta }_m(\cdot ):={\theta }_{m,\mathrm{i}}+\phi (\cdot )\mathrm{\Delta }\mathrm{\Theta },$$

$$\phi :=\frac{{\chi }_m}{C+{\chi }_m}$$

$${\chi }_m(t,S_{n{w}_m}):=\frac{1}{T}{\int }_0^t{S}_{n{w}_m}d\tau ,$$

$$\overline{\chi }(t,S_{nw}):=\frac{1}{T}{\int }_0^t{S}_{nw}d\tau .$$

$$f(R)=\frac{\Big{[}\frac{R-R_{\rm min}}{R_{\rm av}}\Big{]}^{\eta-1}\frac{\eta}{R_{\rm av}}{\rm exp}\Big{(}-\Big{[}\frac{R-R_{\rm min}}{R_{\rm av}}\Big{]}^{\eta}\Big{)}}{1-{\rm exp}\Big{(}-\Big{[}\frac{R_{\rm max}-R_{\rm min}}{R_{\rm av}}\Big{]}^{\eta}\Big{)}}$$

$$\omega =\frac{\overline{\chi }}{{\beta }_1+\overline{\chi }},$$

$$P_{c}=\frac{\overline{\chi}}{2C+\overline{\chi}}\Big{(}P_{c}^{\rm st,f}-P_{c}^{\rm st,i}\Big{)}+P_{c}^{\rm st,i}.$$

$$\omega (S_{nw},\overline{\chi })=\frac{S_w}{\alpha (\overline{\chi })+S_w},$$

$$\alpha (\overline{\chi })={\beta }_2/\overline{\chi },$$

$$P_{c}=\frac{\overline{\chi}S_{w}}{\beta_{2}+\overline{\chi}S_{w}}\Big{(}P_{c}^{\rm st,f}-P_{c}^{\rm st,i}\Big{)}+P_{c}^{\rm st,i}.$$

Full text

1 Department of Mathematics, University of Bergen, P. O. Box 7800, 5020 Bergen, Norway.

2 NORCE, Nygårdsgaten 112, 5008 Bergen, Norway.

Corresponding author: Abay M. Kassa (E-mail: abka@norceresearch.no)

1 Introduction

Wettability plays an important role in many industrial applications, in particular subsurface porous media applications such as enhanced oil recovery (EOR) and CO${}_{\textrm{2}}$ storage [1, 2, 3, 4, 5]. The wetting property of a given multiphase system in porous media is defined by the contact angle, which describes the affinity of one fluid for the surface of the medium over another. The contact angle is determined from surface chemistry and associated forces acting at the molecular scale along the fluid-fluid-solid interface [1]. In porous media applications, microscale wettability determines the strength of pore-scale capillary forces and the movement of interfaces between individual pores. At the macroscale, wettability impacts upscaled quantities and constitutive functions such as residual saturation, relative permeability, and capillary pressure, which in turn affect large-scale multiphase flow behavior.

The standard assumption is that wettability is a static property of the multiphase system. However, the composition of many fluids is capable of provoking the surfaces within pores to undergo wettability alteration (WA) via a change in contact angle. A change in contact angle can alter capillary forces at the pore scale, and thus affect residual saturations of the system [6, 7]. This effect has been exploited extensively in the petroleum industry, where optimal wetting conditions in the reservoir are obtained through a variety of means that includes chemical treatment, foams, surfactants and low-salinity water flooding (see for example [8, 9, 10, 11, 12]).

Wettability is also recognized as a critical factor in geological CO${}_{\textrm{2}}$ sequestration which exerts an important role on caprock performance [13, 14]. WA can influence the ability to prevent CO${}_{\textrm{2}}$ leakage from the reservoir by altering the capillary forces that act to trap CO${}_{\textrm{2}}$ as buoyant phase beneath a low-permeability caprock. This process is highly dependent on CO${}_{\textrm{2}}$ being a strongly non-wetting fluid, and WA may lead to conditions that allow for CO${}_{\textrm{2}}$ leakage [15, 16]. In addition, WA can affect residual saturation, results in an impact on trapping efficiency of injected CO${}_{\textrm{2}}$ [4]. Therefore, reliable quantification of wettability is key to understand the macroscale processes in subsurface applications.

Despite the fact that WA is known to impact macroscale capillarity and relative permeability behavior, few detailed measurements are available to characterize the alteration of the constitutive function themselves. Plug et al., [17] have reported brine-CO${}_{\textrm{2}}$ (gas, liquid) drainage-imbibition experiment and showed capillary instability for supercritical CO${}_{\textrm{2}}$-brine system, meaning that the capillary pressure measurements changes steadily over time. The imbibition curve also exhibited a significant deviation from the expected curve that was not explained by classical scaling arguments. The authors proposed WA as an explanation but did not explicitly measure any changes in contact angle. Additionally, recent experiments measured capillary pressure curves for a silicate sample using a fluid pairing of supercritical CO${}_{\textrm{2}}$ and brine [18]. Repeated drainage-imbibition cycles were performed over a period of 6 months, and a clear reduction in capillary pressure was recorded for each subsequent drainage cycle. The authors also attributed these deviations from the expected capillary curve to a change in wettability of the rock sample over time due to CO${}_{\textrm{2}}$ exposure, similar to aging. This hypothesis was confirmed through observations of a wetting angle increase from $0^{\circ}$ to $75^{\circ}$ after 6-months exposure. It is also reported similar $P_{\textrm{c}}$-$S$ instability and deviation in dolomite/carbonate [19], and quartz [20, 19] sands for scCO${}_{\textrm{2}}$-brine system. More literature on WA and $P_{\textrm{c}}$-$S$ measurements can be found in [14].

The above experiments reveal that capillary pressure curves are not static when significant WA occurs, despite the fact they were performed following the standard multi-step outflow procedure, i.e. where “equilibrium” is obtained after each incremental step in pressure. Therefore, the standard equilibrium capillary pressure models cannot be readily applied without additional dynamics to capture the impact of WA. Capillary pressure dynamics due to WA is fundamentally different from the well known non-equilibrium flow dynamics (for example see [21, 22]) and contact angle hysteresis, i.e. receding and advancing contact angle [23, 24] because these dynamics disappear when the system is at equilibrium.

In petroleum reservoir simulation, there have been some theoretical investigations that assume an instantaneous WA to model the constitutive relations and immiscible flow simulations [25, 26, 3, 27, 28]. The WA mechanism is incorporated by interpolation between the final and initial wetting-state constitutive relations [27, 25]. Lashgari and co-authors [26] have derived an instantaneous WA model from Gibbs and Langmuir adsorption models. The proposed WA model is coupled to $P_{\textrm{c}}$-$S$ and relative permeability curves through residual saturation. These models neglect exposure time dependency of WA and also the impact of pore-scale heterogeneity in wetting properties in the macroscale simulations. In the available literature, only one study [29] has included the effect of exposure time on WA and constitutive relations for core scale simulation. But this numerical study does not sufficiently incorporate or upscale pore-scale processes to macroscale constitutive laws.

To our knowledge, a rigorous mathematical characterization of dynamics in $P_{\textrm{c}}$-$S$ functions introduced by exposure to a WA agent has not been previously performed. In this paper, we develop a new model for dynamic capillary pressure that is derived using direct simulation of $P_{\textrm{c}}$-$S$ curves from a pore-scale model. WA is introduced at the pore scale using a mechanistic model for contact angle change as a function of exposure time to a WA fluid. This mechanistic model is not intended as a replacement for detailed surface chemistry models, but is used merely as a convenient mathematical way to represent changing contact angle in time. We acknowledge that the field of wettability alteration at the micro- and nano-scales is an active area of research covering a broad range of complexity. We emphasize that this study does not directly address wettability alteration at pore surfaces, which is well beyond the scope of this paper.

The pore-scale flow model is represented by a cylindrical bundle-of-tubes model. This pore-scale approach is simplistic in that it does not represent interactions between pores or realistic porous geometry (i.e. converging-diverging throat diameters). Additionally, there is no natural hysteresis or residual saturation with cylindrical tubes. However, a simple pore-scale model is sufficient as a first step towards develop time-dependent correlated expressions for $P_{\textrm{c}}$-$S$. It is also important to investigate the impact of wettability alteration apart from other complicating processes that can occur in more realistic pore-network models.

2 Approach

The extended $P_{c}$-$S_{w}$ relationship introduces a dynamic component that captures the changing wettability as measured by the deviation of the dynamic capillary pressure from the equilibrium (static) capillary pressure. This relationship can be described as follows:

[TABLE]

where $P_{c}^{\rm st,i}$ represents the capillary pressure for the system given a static initial wetting state, while $f^{\rm dyn}$ represents the deviation from the static state. The static curve can be characterized by a number of well-known capillary pressure models, such as van Genuchten or Brooks-Corey. Static systems studied here are described by the Brooks-Corey model,

[TABLE]

where $c_{w}$ is the entry pressure, $1/a_{w}$ is the pore-size distribution index, whereas $S_{wc}$ is the residual water saturation. The static curve is equivalent to the equilibrium capillary pressure for systems that do not undergo wettability alteration.

The objective of this study is to characterize and quantify the dynamic term $f^{\rm dyn}$, the key term of interest in the $P_{\textrm{c}}$-$S$ model, for a system that undergoes wettability alteration. We propose an interpolation model can be applied for this system, where the dynamic term acts to shift the capillary pressure between two end wetting states in time. This type of model has been used in petroleum reservoir simulation [27, 25, 3, 26]. To obtain an interpolation model, the dynamic component in Equation (1) can be scaled by the difference between two static curves, each representing the initial and final wetting-state capillary pressure curves, to give a non-dimensional quantity we call the dynamic coefficient, which is defined as follows,

[TABLE]

where $P_{c}^{\rm st,f}$ is the final wetting state capillary pressure. This can be substituted into Equation (1) to obtain a dynamic interpolation model for capillary pressure,

[TABLE]

We note that in the model presented above, we have defined the “total” capillary pressure $P_{c}$ as simply the measured difference in phase pressures at any point in time. In a reservoir simulation, this would be capillary pressure in a given grid cell, whereas in a multi-step outflow experiment, it corresponds to the pressure drop across the experimental cell. For quasi-static displacement in a bundle of capillary tubes, $P_{c}$ is the difference in reservoir pressures, i.e., $P_{c}=P_{l}^{res}-P_{r}^{res}$, (see Figure 1). The exact nature of $\omega$ and its functional dependencies can only be determined from a full characterization of $P_{\textrm{c}}$-$S$ curves under different conditions. These curves can be derived from laboratory experiments, but this approach is costly and time-consuming. Alternatively, one may take a more theoretical approach by simulating $P_{\textrm{c}}$-$S$ curves using a pore-scale model that includes the impact of WA.

For a system that undergoes WA, a significant change in contact angle could lead to the wetting phase becoming non-wetting and vice versa. For clarity, we will continue to use $w$ subscript for the phase that was originally wetting and $nw$ for the phase that is originally non-wetting regardless of the actual state of wettability in the system.

2.1 Capillary flow in cylindrical tubes

There are various choices of pore-scale models available. The easiest to implement and analyze is the bundle of tubes (BoT) model which is a collection of capillary tubes with a distribution of radii. Herein, we describe the model and the approach for implementation of time-dependent WA within individual capillary tube.

The BoT model has been applied in many immiscible multiphase flow studies and is employed to study constitutive relations and fluid dynamics in general [30, 22, 31]. The BoT model is a popular approach due to the simplicity of implementation and the ability to study the balance of energy and forces directly at the scale of interfaces. The average behaviour of a BoT model can then be used to inform better constitutive functions at the macroscale. However, direct use of simulated capillary or relative permeabilities curves is usually not advised. This is because a collection of capillary tubes is not representative of the interacting nature of pores and pore throats. The geometry is over-simplified with respect to converging and diverging pore throats. Some of the shortcomings of BoT models can be addressed through interactive BoT models [32], or tubes of more complex geometry. These more advanced BoT models are beyond the scope of this paper.

In this paper, we consider cylindrical bundle of tubes having length, L, and is shown in Figure 1.

These tubes are designed to connect wetting (right) and non-wetting (left) reservoirs with pressures labeled as $P_{\rm r}^{\rm res.}$ and $P_{l}^{\rm res.}$ respectively. Let the reservoir pressures difference be defined as:

[TABLE]

Equation (5) indicates that the pressure drop across the bundle is equal to the measured capillary pressure $P_{c}$ defined previously in Equation (1).

We consider the tubes to be filled with water and are water-wet initially. To displace the water phase, the pressure drop has to exceed the entry capillary pressure [22]:

[TABLE]

where $P_{c}^{\rm mic}(R_{m},\theta_{m})$ is given by the Young’s equation

[TABLE]

where $R_{m}$ is the $m$-th tube radius, $N$ stands for number of tubes, $\sigma$ is fluid-fluid interfacial tension (IFT) and $\theta_{m}$ measures the wetting property of the rock surface (it may be static or dynamic in time per the $m$-th tube). As long as condition (6) is satisfied, the fluid movement across the length of a given tube can be approximated by the Washburn flow model,

[TABLE]

where $q_{m}$ is the interface velocity, and $\mu_{\rm nw}$ and $\mu_{\rm w}$ are non-wetting and wetting fluid viscosities, respectively. The above expression for interface velocity is equivalent to the rate of change of the interface location $x^{i}_{m}$ in tube $m$, which can be defined as ${dx_{m}^{i}}/{dt}$. From equation (8), one can then determine the required time to reach a specified interface position, and subsequently the time to drain each individual tube in a multi-step outflow experiment.

2.2 Pore-scale time-dependent wettability model

In this study, the WA agent is defined as either the non-wetting fluid itself or some reactive component therein that is transferred to the wetting fluid across a fluid interface. We consider an alteration process that continues until the ultimate wetting state is reached. The alteration is permanent, but can also be halted at some intermediate wettability state if the WA agent is removed from the system before the final wetting state is reached. If the agent is reintroduced at some later point, alteration continues until the final state.

To this end, we introduce a general functional form of pore-scale WA mechanism,

[TABLE]

where ${\rm\Delta\Theta}=\theta_{m,\rm f}-\theta_{m,\rm i}$, $\theta_{m,\rm f}$ and $\theta_{m,\rm i}$ are the ultimate and initial contact angles respectively. Note that $\varphi$ can be any function form that is properly designed to capture the WA processes over a period of time. As discussed in the introduction, the process of contact angle change at the pore-scale involves surface chemistry acting at micro and nano-scale. The mechanistic model in Equation (9) is not intended to capture the full complexity at the pore-scale, but is simply a convenient mathematical form that gives smooth change in contact angle over time between an initial and final contact angle. Detailed laboratory work would needed to further justify the use of this model, which is beyond the scope of this paper.

Though it is well known that WA impacts the reservoir system, theoretical investigation and laboratory measurement on time-dependent WA is very limited. However, coupling the instantaneous WA processes to the Langmuir adsorption model have been successfully implemented in EOR simulation technology [27, 25, 3, 26]. In this paper, the parameter $\varphi(\cdot)$ in equation (9) is assumed to evolve according to Langmuir-type adsorption as a test case. Thus, $\varphi$ can be defined as follows:

[TABLE]

where $C$ is non-dimensional parameter that controls the speed and extent of alteration from water-wet to intermediate-wet or hydrophobic system. The non-linear variable, $\chi_{m}$ is a measure of the local exposure time of each individual tube $m$, defined as

[TABLE]

where $T$ is a pre-specified characteristic time. In this paper, the total time to complete the possible drainage-imbibition cycles from the initial to final wetting state is used as a characteristic time.

For convenience, we may also define $\overline{\chi}$ as a measure of exposure time across the entire bundle (i.e. REV),

[TABLE]

Given the above model in Equation (10), we may consider two WA approaches, namely the uniform and non-uniform alteration models.

In the case of uniform WA, pores are assumed to be altered identically in time due to transfer of the WA agent by dissolution from neighboring pores. In this case, the local exposure time is equivalent to the average exposure time, such that $\chi_{m}=\overline{\chi}$ in Equation (10). Alternatively, we may restrict WA to only those pores that have been flooded by the non-wetting fluid. In this case, we must keep track of the local exposure time of each tube within the bundle to evaluate Equation (10). This process gives rise to a non-uniform or mixed alteration of the pores until such time that all pores have reached the final wetting state.

2.3 Simulation approach

The uniform and non-uniform approaches for contact angle evolution above are coupled into the BoT model following the algorithm pictured in Figure 2.

For a drainage cycle, the pressure drop is increased incrementally. Upon each increase, the respective tubes are drained according to the entry pressure inequality. The Washburn model, Equation (8), is used to calculate the elapsed time to drain each tube. Afterwards the contact angle (and thus $P_{c}^{\rm mic}$) is updated according to Equation (9). For the uniform model, the contact angle change is applied uniformly, while in the non-uniform model, the change is applied according to local exposure time. The average saturation is recorded and the cycle is repeated until all tubes have drained.

During imbibition, the process is reversed. For the non-uniform model, the contact angle continues to change locally until the non-wetting phase is displaced by the wetting phase. For the uniform model, the alteration occurs continuously throughout imbibition according to the average exposure time $\overline{\chi}$.

Repeated drainage and imbibition cycles are run according to this algorithm until the final wetting state is reached. At the completion of the numerical experiment, we obtain a set of $P_{\textrm{c}}$-$S$ “data points” that can be plotted in the usual way.

2.4 Correlation approach

Once the capillary pressure and associated dynamics have been simulated, the resulting curves used to quantify the dynamic coefficient in the interpolation function in Equation (4). The goal is develop a correlated model that involves only a single parameter, and this parameter should have a clear relationship to the pore-scale WA process such that changes in the pore-scale WA model parameter $C$ can be directly accounted for in the dynamic coefficient.

We note that for the interpolation model, the dynamics are contained in the coefficient $\omega$ that interpolates between two static capillary pressure curves. Therefore before correlating this model, we determine the functional form of the static wetting states using a simple power-law model.

3 Results

In this section, we present the results of simulated capillary pressure experiments using pore-scale models represented by a bundle of tubes following the approach described above. First, the static capillary pressure is simulated for both the initial and final wetting states. Then, simulations are performed on the model subject to WA. Two cases are performed: uniform and non-uniform alteration. The correlated $P_{\textrm{c}}$-$S$ models for the static cases are presented, followed by the correlated models for the dynamic capillary pressure curves.

3.1 Bundle of tubes model set-up

The pore scale is described by a BoT model (see Section 2.1). Each tube in the BoT is assigned a different radius R, with the radii drawn from a truncated two-parameter Weibull distribution [31]

[TABLE]

where $R_{\rm max}$, $R_{\rm min}$, and $R_{\rm av}$ are the pore radii of the largest, smallest, and average pore sizes, respectively, and $\eta$ is a dimensionless parameter. The rock parameters and fluid properties are listed in Table 1.

These parameters are coupled to the BoT model to simulate capillary pressure curves.

3.2 Static capillary pressure for end wetting states

A starting point for the dynamic capillary-pressure models presented in Section 2 (Equations (1) and (4)) is characterizing the capillary pressure curves for the end wetting states. Given the same tube geometry and fluid pairing described above, the capillary pressure functions are simulated under static conditions for both the initial and final wetting states using the BoT model.

Only a single drainage experiment is needed in the static case to fully characterize the capillary pressure curve. This is due to the simplicity of the BoT model and lack of residual saturation in cylindrical tubes. Therefore, hysteresis is not possible if the contact angles in the tubes (and other parameters) are held constant.

The simulated static curves can be correlated with the Brooks-Corey model (2). The resulting correlations can be found in Figure 3, while fitted parameters for the Brooks-Corey model can be found in Table 2.

Figure 3 compares the Brooks-Corey formula (2) and the capillary pressure curves associated with static contact angles (initial and final wetting states).

The Brooks-Corey correlation gives an excellent match to the simulated $P_{\textrm{c}}$-$S$ data under static conditions. We observe that the pore-size distribution index $a_{w}$ for the initial and final wetting states are the same for both end wetting states, which is expected since the same distribution of tube radii is used in both cases. While, on the other hand, the coefficient $c_{w}$ decreases by a factor of 0.85 from the initial to final wetting state, corresponding to a decrease in macroscale capillary entry pressure. Leverett-J scaling theory [34] predicts that entry pressure scales by $\cos\theta$, which agrees nicely with the reduction in $\cos\theta$ by a factor of 0.83 for a contact angle change from 0 to 80 degrees.

We reiterate that for the static case where no wettability alteration occurs, the Brooks-Corey model describes both drainage and imbibition for the bundle of tubes.

3.3 Dynamic capillary pressure for uniform wettability alteration

Here, we present results for the first dynamic case (Figure 4) in which wettability is altered uniformly throughout the bundle by considering the pore-scale parameter $C=0.005$. Two dynamic drainage-imbibition cycles are simulated using the quasi-static BoT approach over a period of approximately 7 months. During this time, the capillary pressure steadily decreases from the initial to the final wetting state with each subsequent cycle. Having reached the final wetting state, any additional drainage-imbibition cycle would follow along the static curve for the final wetting state. We observe that wettability-induced dynamics also introduces an apparent hysteresis in the $P_{\textrm{c}}$-$S$ data. We recall that a simple BoT model cannot exhibit hysteresis under static wettability conditions.

The dynamics in capillarity are coupled to a continuous change in contact angle between the two wetting states, which is shown in Figure 5 for the same two cycles. We observe that $\theta$ undergoes a rapid increase by the end of the first drainage, followed by a slower evolution towards the final wetting state. Since there is no spatial variability in wettability at the pore scale in the uniform case, smaller tubes are altered at the same rate as larger tubes. This is evident in the associated $P_{\textrm{c}}$-$S$ data where we observe substantial alteration of the first drainage path at low wetting saturation. The alteration from the initial wetting condition is fast, leading to non-trivial and almost non-monotonic dynamics in capillarity at early time.

The saturation history is an important factor, defining the length of exposure in individual tubes. Figure 6a shows that saturation history is driven by the time required to displace fluid through each tube in sequence, i.e. from largest to smallest during drainage and vice versa during imbibition. From integration of this curve, we obtain the exposure time $\overline{\chi}$ (Figure 6b) which controls the contact angle change shown previously in Figure 5. We recall that for the uniform model, contact angle change at the pore-scale (Equation (10)) is driven by the average saturation not the local saturation.

Figure 6a also indicates the time taken for each drainage-imbibition cycle, with the first cycle takes approximately 1 month to complete. Since the displacement process is speeding up with increasing contact angle, up to 12 cycles would be required in this manner before the end wetting state is reached. In the interest of presentation, we have controlled the pressure drop during the final imbibition experiment so that exposure time was lengthened in order to reach the final wetting state by the end of the second cycle.

We have observed above that capillary pressure is no longer a unique function of saturation as it is in the static case for this simple geometry. However, if $P_{\textrm{c}}$-$S$ is transformed to the temporal domain by plotting against $\overline{\chi}$, we obtain a unique function with respect to exposure time. However, Figure 7a shows that the capillary pressure relates with the temporal domain $\overline{\chi}$ non-monotonically i.e., the capillary pressure increases and decreases with each drainage-imbibition cycle. Similarly, the associated contact angle change is plotted along with exposure time (Figure 7b), which is exactly the function given in Equation (10) at the pore level. The value of $\theta$ is smoothly increasing in time according to the pore-scale model.

3.3.1 Dynamic capillary pressure model development

Following the approach discussed in Section 2, we applied Equation (3) to determine the dynamic coefficient $\omega$. As shown in Figure 8a, we observe that $\omega$ is a non-unique function of wetting-phase saturation but with values that are continuously increasing as the dynamic capillary pressure approaches the final wetting state. This non-uniqueness behavior makes challenging to propose a functional form for $\omega$-$S_{w}$ relation.

On the other hand, Figure 8b shows that $\omega$ is smoothly increasing, uniquely related to $\overline{\chi}$, and mimics the functional form of the pore-scale model in Equation (9) (see Figure 7b). Given this insight, we propose an adsorption-type model to fit the interpolation coefficient $\omega$ to the average exposure time:

[TABLE]

where $\beta_{1}$ is a fitting parameter obtained from the best fit to simulated data in Figure 8b. For this particular case, the calibrated parameter is estimated to be $\beta_{1}=0.01$.

The proposed model for $\omega$ in Equation (14) is incorporated into the dynamic capillary pressure model proposed in Equation (4) and compared against the simulated dynamic capillary pressure curves. The resulting fit, shown in Figure 9 indicates that a one-parameter model for the dynamic coefficient agrees well with the simulated data. The correlation coefficient for this comparison is $R^{2}=0.9921$.

We hypothesize that the value of $\beta_{1}$ in this model is dependent on the dynamics of wettability alteration at the pore-scale. Therefore, we investigate how the parameter $\beta_{1}$ evolves with changes in the pore-scale parameter $C$ in Equation (9). We repeat the capillary pressure simulations for different values of the pore-scale parameter $C$ and determine the correlated value of $\beta_{1}$ in each case. Figure 10 shows that the interpolation model parameter is directly proportional to the pore-scale model parameter, with a proportionality constant of 2. Thus, the relationship $\beta_{1}=2C$ can be used to predict the macroscale parameter directly from knowledge of the pore-scale process.

The final dynamic capillary pressure model can now be obtained by incorporating the above model for $\omega$ with the established relationship for $\beta_{1}$ into the general form of the interpolation model in Equation (4):

[TABLE]

In its final form, the dynamic capillary pressure model in Equation (15) is solely a function of the pore-scale parameter $C$. This parameter must be determined from laboratory experiments for a given sample exposed to a WA agent.

3.3.2 Dynamic model applicability

We note that the saturation history in Figure 6a used to generate the $P_{\textrm{c}}$-$S$ curves in Figure 4 can be thought of as one arbitrary path of an infinite number of possible paths. If a different path had been chosen, such as a flow reversal at intermediate saturation or a prolonged exposure time at a given saturation, it would result in entirely different capillary pressure dynamics.

We generate many different $P_{\textrm{c}}$-$S$ curves by taking numerous different paths in the saturation-time domain. The resulting simulated data forms a surface with respect to saturation and exposure time as shown in Figure 11a. We then apply the calibrated dynamic model from Equation (15) to the same saturation-time paths used to generate the $P_{c}$-$S$-$\overline{\chi}$ surface. The difference between the calibrated model and the simulated data is shown in Figure 11b. A good comparison of the dynamic model to simulated data demonstrates that model calibration to a single saturation-time path is robust enough to be applied to any possible path.

3.4 Dynamic capillary pressure for non-uniform wettability alteration

We present the second case where WA evolves throughout the bundle non-uniformly at the pore level with the pore-scale parameter $C=1\times 10^{-5}$. In this case, we performed four WA induced dynamic drainage-imbibition displacements and are shown in Figure 12.

The first drainage curve follows exactly along the static curve for the initial wetting state, which is expected since the entry pressure to each individual tube is not altered before the tube has drained. Only after the tubes are drained by the WA agent can the contact angle change in each tube. Therefore, it is only the imbibition data that show the corresponding reduction in capillarity. At the end of imbibition, the system becomes quasi-static in the sense of alteration, which means the wettability of a given tube stops being further altered until a new drainage cycle has begun. This implies that each new drainage curve must follow along the previous imbibition curve and thus, imbibition and drainage paths overlap in Figure 12. This effect is in contrast to the uniform WA case, where each drainage and imbibition curves follow different paths.

The capillary pressure instability in Figure 12 is attributed to the local change of wettability in each tube. The larger tubes drain first and are imbibed last, resulting in a long local exposure time compared to the smaller tubes that drain last and imbibe first. Therefore, the initial wetting state persists in the smaller tubes. Figure 13 shows the resulting heterogeneous contact angle change in time across the bundle.

We observe that the impact of fluid history and exposure time on the wetting property of the bundle and the associated capillarity is distinct from the uniform WA case. Here, the dynamic capillary pressure continues to rise to the initial wetting state capillary pressure curve for low water saturation. This is due to the reason that the smaller tubes were exposed to the WA agent for shorter time than the larger tubes. This situation persists until the bundle approaches a uniform wetting condition corresponding to the final wetting state.

As we noticed for the uniform WA case, we also observe $P_{\textrm{c}}$-$S$ hysteresis induced by non-uniform WA in Figure 12. Therefore, we can transform the $P_{\textrm{c}}$-$S$ curve to the temporal domain $\overline{\chi}$ in the same way as before. The temporal domain $\overline{\chi}$ is shown in Figure 14a, and is obtained by integrating the averaged non-wetting fluid history in Figure 14b over exposure period of 33 months. For the sake of clarity and the same reason as the uniform WA case, we controlled the pressure drop after the third cycle so that the tubes experience prolonged exposure time in order to complete the alteration processes within few number of cycles.

The averaged variable $\overline{\chi}$ calculated from the saturation history is then used to transform the capillary pressure curve into temporal space. In Figure 15, we observe a unique relation between the capillary pressure and $\overline{\chi}$ as hypothesised above. However, unlike the $P_{\textrm{c}}$-$S$ relation, the $P_{\textrm{c}}$-$\overline{\chi}$ relation is not monotonic which rises and falls along exposure time.

3.4.1 Dynamic capillary pressure model development

Given the capillary pressure results above, the dynamic coefficient $\omega$ is calculated according to Equation (3), and plotted as a function of $S_{nw}$ and $\overline{\chi}$ in Figure 16a and 16b respectively. Unlike the uniform WA case, seen in Figure 8b, the coefficient $\omega$ is not monotonically increasing in $\overline{\chi}$, but instead it rises and falls. This makes it challenging to propose a functional relation between the dynamic coefficient $\omega$ and $\overline{\chi}$ directly as we did for the uniform WA case.

On the other hand, $\omega$ in Figure 16a has a well-behaved curvature for each drainage-imbibition cycle along the saturation history. Further, the dynamic coefficient as a function of saturation is increasing with increased exposure time. Given these insights, we proposed a model for the dynamic coefficient,

[TABLE]

where $S_{w}=1-S_{nw}$, and $\alpha$ controls the curvature of the $\omega$-$S_{nw}$ curve for each drainage-imbibition cycle. Since $\omega$ is increasing function of exposure time, $\alpha$ should decrease along the averaged variable $\overline{\chi}$.

The function form of $\omega$ in (16) is then matched with the $\omega$-$S_{nw}$ data to analyze the dynamics of $\alpha$ along $\overline{\chi}$. The obtained $\alpha$-$\overline{\chi}$ relation is decreasing as hypothesized and in particular has a form that can be governed as,

[TABLE]

where $\beta_{2}$ is non-dimensional fitting parameter. For this particular simulation the parameter $\beta_{2}$ is estimated to be 0.004 for the four of drainage-imbibition cycles.

The suggested models (16) and (17) are substituted into the interpolation model (4) to give:

[TABLE]

The calibrated dynamic capillary pressure model is compared with the simulated time-dependent capillary pressure data and depicted in Figure 17. We observe that the proposed dynamic model agrees well with simulated dynamic capillary pressure curves. Thus, we have obtained a single valued single parameter model that describes the evolution of dynamic capillarity over multiple drainage-imbibition cycles rather than using a model consisting of multiple parameters that change with each cycle (or hysteresis models).

3.4.2 Dynamic model applicability

As we point out in the uniform WA case, the saturation paths above are arbitrary paths out of infinitely many $P_{\textrm{c}}$-$S$-$\overline{\chi}$ paths. Keeping this in mind, we generate a data for dynamic capillary pressure by considering possibly many paths in $S_{w}\times\overline{\chi}$ domain, resulting in a capillary pressure surface shown in Figure 18a.

Following the paths used to generate the simulated $P_{\textrm{c}}$-$S$-$\overline{\chi}$ surface, we applied the calibrated model (18) to predict the simulated capillary pressure data. Figure 18b shows the actual difference between the simulated data and the model predicted data. The comparison demonstrates that a single parameter value model calibration to a single arbitrary path in $S_{w}\times\overline{\chi}$ domain can be applied to any possible path in the domain. In general, this simple model is efficient to predict the the impact of non-uniform WA at the pore scale on the capillary pressure at the macroscale for any given exposure time and saturation path.

Given the successful WA induced dynamic capillary pressure model, we perform the same exercise as before where we vary the pore-scale parameter $C$, from Equation (9), and observe its impact on the macroscale parameter $\beta_{2}$ in the interpolation model above. The relation between those parameters is obtained and depicted in Figure 19. A power law model, $\beta_{2}=b_{1}C^{b_{2}}$ is observed between these parameters and the model parameters are estimated to be $b_{1}=3.3\times 10^{6}$ and $b_{2}=1.8$. The model for $\beta_{2}$ can be substituted into model (18) to characterize the dynamics of the capillary pressure directly from the pore scale phenomenon.

To summarize, we have obtained a reliable model to describe the impact of the pore-scale wettability alteration on capillary pressure dynamics.

4 Discussion and Conclusions

We have proposed a new interpolation based model for dynamic capillary pressure caused by time-dependent wettability alteration. The model has been correlated to pore-scale simulations of capillary pressure-saturation curves where a time-dependent model for contact angle change has been applied directly at the pore scale. The resulting correlation shows that a simple interpolation-based model is highly effective at capturing WA induced dynamics in macroscale capillary pressure by a single value parameter. In the interest of completeness, we explored other types of models to capture capillary pressure dynamics, including the mixed-wet model of Skjæveland et al [33]. For brevity, we do not report the results of that separate study herein. We found that although other models could be calibrated with reasonable accuracy, they all involved more than one calibration parameter (up to four) that need to be adjusted in each drainage-imbibition cycle. Therefore, the single-parameter single-valued interpolation model presented in this study is the preferred model due to its reliability.

The proposed interpolation model is a macroscale model that allows for change in capillary pressure as a function of macroscale variables–saturation and exposure time to a WA agent. We recall that exposure time is simply the integration of saturation history over time. The model consists of three main components–two capillary pressure functions at the initial and final wetting state and a dynamic interpolation coefficient that moves from one state to the other. The initial and final capillary pressure functions can be determined a priori from static experiments using inert fluids. In this study, the initial and final states are represented by classical Brooks-Corey functions. The dynamic coefficient is thus the only variable correlated to dynamic capillary pressure simulations. In this study, we have shown that the coefficient can be easily correlated to saturation and exposure time via a single parameter.

The resulting correlation is dependent on the conditions for wettability alteration. Two conditions are explored, the first is dissolution of the agent into the wetting phase that alters the contact angle uniformly across the REV. The second is contact angle change only through direct contact with an altering agent as a separate non-wetting phase, resulting in non-uniform contact angle change throughout the REV. The differences in the two wettability alteration mechanisms changes the complexity of the resulting capillary dynamics. In the uniform case, the dynamic coefficient can be correlated to exposure time through a sorption-type model, which seems to be a natural result given the contact angle change at the pore scale is also based on a sorption model. In the non-uniform case, the dynamic coefficient has a similar sorption form with increased exposure time, but now with the product of saturation and exposure time as the dynamic variable. The additional complexity is needed to draw the capillary pressure curve back to the initial wetting at low saturation (a region of the $P_{\textrm{c}}$-$S$ curve dominated by smaller pores where the contact angle takes longer to change).

An important result of this study is quantifying the link between the pore- and macroscale. We showed that by varying the parameter that alters the speed and extent of contact angle change in each individual pore, we could predict the resulting impact on dynamic capillary pressure. In fact, in both the uniform and non-uniform cases, there is a very simple scaling from the pore-scale and macroscale parameters. In the uniform case, the two parameters are directly proportional, while in the non-uniform case, the macroscale parameter scales with the pore-scale parameter via a power law. The implication of this result is that by knowing the mechanism that controls contact angle at the pore-scale, which can be obtained by a relatively simple batch experiment, we can quantify a priori the macroscale dynamics without having to perform pore-scale simulations. This is an important generalization and valuable for making use of experimental data to inform macroscale constitutive functions.

The success of the interpolation model in this study is a natural development from previous studies that incorporates the interpolation model directly into reservoir simulation of wettability alteration. In those studies, the model was matched directly to macroscale data in a heuristic manner. The contribution of this study is to quantify the pore-scale underpinnings to the interpolation model through a direct and systematic manner, thus confirming the validity of this type of model for use in macroscale simulation. The insights gained in this study into capillary pressure dynamics due to wettability alteration. Further study is needed to calibrate the macroscale dynamic parameter to actual data.

There are additional aspects that were not addressed in this study. First, is whether the effect of time-dependent wettability alteration can be quantified in a similar manner for relative permeability. For brevity, we did not explore this additional constitutive relationship in this paper. This can be explored in future work. An additional quantity of interest is residual saturation. Quantifying dynamics and stability of residual wetting and non-wetting saturation is an active area of research. The scope of this work involved a bundle-of-tubes pore-scale model, which due to its simple nature cannot model trapped phases. Therefore, a more realistic pore-network model would need to be implemented, which is the subject of ongoing work.

The results of this study have important implications for future laboratory investigations. Previous experimental studies have shown that wettability alteration can cause a significant decrease in capillarity that occurs at the same time scale as the experiment. Our work further highlights this time-dependent effect and shows that proper quantification of the dynamics requires capillary pressure data from both inert and reactive fluid pairings. In addition, there should be separate batch measurements on mineral surfaces to characterize the mechanistic model for contact angle change at the pore scale, such as the one proposed herein. Also important is whether the alteration mechanism occurs by dissolution of a WA agent or by direct contact with a non-wetting phase. Such data are important for determining the macroscale parameters needed in the interpolation-based model.

At the reservoir scale, there is a potential for dynamics in capillary pressure to have a significant impact on fluid flow that should be accounted for during simulations. For instance, for CO${}_{\textrm{2}}$ storage could there be a change in wetting that could degrade the integrity of the sealing caprock and allow for seepage of CO${}_{\textrm{2}}$ out of the storage reservoir, and thus further investigation is warranted. This will require implementation of dynamic capillary pressure functions into reservoir simulation, with possible wettability reversal that could prove challenging for numerical consistency and robustness. Modified simulation methods can then be used in the future to quantify the impact of capillary dynamics on reservoir processes for both hydrocarbon recovery and CO${}_{\textrm{2}}$ storage.

Acknowledgement

The research is supported by CHI project funded by the CLIMIT program of the Research Council of Norway.