Modeling of relative permeabilities including dynamic wettability transition zones

Wettability is a pore-scale property that impacts the relative movement and distribution of fluids in a porous medium. There are reservoir fluids that provoke the surface within pores to undergo a wettability change. This wettability change, in turn, alters the dynamics of relative permeabilities at the Darcy scale. Thus, modeling the impact of wettability change in relative permeabilities is essential to understand fluids interaction in porous media. In this study, we include time-dependent wettability change into the relative permeability--saturation relation by modifying the existing relative permeability function. To do so, we assume the wettability change is represented by the sorption-based model that is exposure time and chemistry dependent. This pore-scale model is then coupled with a triangular bundle-of-tubes model to simulate exposure time-dependent relative permeabilities data. The simulated data is used to characterize and quantify the wettability dynamics in the relative permeability--saturation curves. This study further shows the importance of accurate prediction of the relative permeability in a dynamically altering porous medium.


Introduction
Wettability alteration (WA) plays an important role in many industrial applications such as microfluidics nanoprinting, enhanced oil recovery (EOR), and CO 2 storage [1,2,3,4,5]. Wettability refers to the tendency of one fluid over the others to spread on or adhere to a solid surface [35,1] and is defined by the fluid-fluid contact angle (CA). This pore-scale property regulates the distribution of 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 as residual 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 constitutive relations refers to the relative permeability-and capillary pressure-saturation relations) is of great importance.
Wettability is assumed to be static in time and uniform in space. However, wettability is a dynamic process that depends on surface chemistry, composition of fluids, exposure time, and reservoir conditions (pressure and temperature) to name a few [41,42,38,11,10,9,43]. Experiments on crude oil/brine/rock systems have shown that adsorption of active components from the crude oil is able to change the wettability of the sample porous medium from water-wet to intermediate-wet system [41,44,9]. It is also hypothesized that the oil reservoir may be more oil-wet than what is observed from the experiment. This is because the adsorption time of the experiment period was much less than the age of the oil in the reservoir. Furthermore, CO 2 is one of the reservoir fluids which 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 final wetting-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 condition could be created by rock mineral and exposure history differences. This is due to the fact that a pore surface exposed to the WA agent may be altered to a new wetting condition, while the unexposed surface keeps the initial wetting state [54,7]. This creates a mixed-wet condition even within a single pore and was observed and explained first by Salathiel et al. [41] in the 1970s. Usually, WA is 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 and months [14,18,9,55]. In this regard, the dynamic-wet phase can be a function of exposure time in addition to the WA agent concentration.
The WA process may result in a saturation function alteration for subsequent drainage-imbibition displacements and thus cause hysteresis in constitutive relations [56,36,57,58,59]. For instance, core-flooding measurements for (supercritical or gas) CO 2 -water system have indicated that WAinduced alteration in the residual saturation and capillary pressure curves despite the fact that they were measured following a standard procedure, i.e., where "pressure equilibration" is obtained after each increment in pressure [17,19,14,20,45,13]. In these measurements, a steadily change in capillary pressure function over time was observed. More importantly, the capillary pressure deviation from the initial-wet state curve could not be explained by classical scaling arguments. The instability 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 predict relative permeability and capillary pressure under dynamic, long-term WA. One alternative is to use mixed-wet model, e.g. Kjosavik et al. [60] and Lomeland et al. [61], that capture the static heterogeneity of wettability in the relative permeabilities. The main feature of these models is their 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 relative permeabilities. The first class of these models involves a heuristic approach that interpolates between the 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 wetting states are characterized by standard functions, e.g. Brooks-Corey [63] or van Genuchten [64]). The other approach incorporates the effect of the instantaneous WA into the relative permeabilities through the residual saturation directly [26]. To date, only Al-Mutairi et al. [29] have considered the effect of time-dependent WA in both the relative permeabilities and capillary pressure functions explicitly. However, their model does not sufficiently incorporate or upscale the WA processes to core-scale laws.
Appropriate upscaling of the pore-scale time-dependent WA process connected to the capillary pressure function was the subject of our recent work [65]. There, WA dynamics were upscaled by introducing a mechanistic time-dependent CA model at the pore-level that was coupled with a cylindrical bundle-of-tubes model and used to simulate capillary pressure curves for drainage and imbibition displacements. The simulated data was used to formulate and quantify a interpolationbased capillary pressure model at the Darcy scale. The new dynamic model resolves the existing interpolation models used in the studies of reservoir simulation [28,25,3,27,62] by including the dynamics in time and quantifying the pore-scale WA process to the interpolation model in a systematic manner.
One may consider employing a similar approach to [65] and an interpolation-type model to capture the pore-scale underpinnings of WA in the relative permeability behaviors. However, timedependent WA may impact the capillary pressure and relative permeabilities in different ways. As observed in [65], WA has a direct impact on the entry pressure in each pore and reflects it at the Darcy scale. Furthermore, a small change in CA exerts a large impact on the dynamics of the capillary pressure function. However, the relative permeability alteration occurs when the WA affects the pore filling/draining orders of pore-sizes. This may lead to a longer exposure time to observe a relative permeability deviation from the initial-wet state curve. Furthermore, unlike the capillary function, the relative permeability curves are constrained between zero and one for any change of wettability. These features of the relative permeability may impact the modeling approach to upscale the pore-scale WA process to the relative permeability behavior.
To our knowledge, a physically reliable model to characterize a prolonged exposure time-dependent WA induced dynamics in the relative permeability behaviors has not been proposed yet. This paper revises and extends the approach discussed in [65] to develop a relative permeability model that includes pore-level time-dependent WA processes. Section 2 summarizes two possible approaches that can be applied to upscale the impact of time-dependent WA in the relative permeabilities. The fluid-fluid CA change is designed as a function of exposure time to the WA agent at the pore level to measure the WA process. This model is coupled with a pore-scale, triangular bundle-of-tubes, model to simulate time-dependent WA induced relative permeability curves. These curves are presented in Section 3 and are used to evaluate the modeling approaches hypothesized in Section 2.

Modeling and simulation approach
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 wettingstate as: or, the dynamic term can be correlated with the parameters of the standard models where a(·) and b(·) are fitting parameters that change along exposure time, whereas f d α represents the WA induced dynamic component, and the subscript α ∈ {w, n} represents the wetting and non-wetting phases.
In this study, we explore both approaches in Eqs. (1) and (2) to quantify and characterize f d α in the relative permeabilities for a system that undergoes a WA. From Equation (1), we propose an interpolation model following our previous work [65], where the dynamic component is designed to interpolate between two end wetting-state curves. To obtain an interpolation model, the dynamic component in Eq. (1) can be scaled by the difference between the initial and final wetting-state relative permeability curves. The resulting quantity is non-dimensional and referred to as the dynamic coefficient, where the superscript i and f represents relative permeabilities at the initial and final wetting-states respectively. This can be substituted into Eq. (1) to obtain dynamic relative permeability models where ω α , the dynamic coefficient, is responsible for capturing the wettability dynamics at the macroscale. Similar model to Eq. (4) were employed to include the impact of instantaneous WA into the 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 relative permeability function through the model parameters. This can be done by formulating the parameters a(·) and b(·) as a function of exposure time and WA agent in similar fashion as ω α . This approach is motivated by the fact that the parameters in the standard relative permeability models are 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-known relative permeability models such as van Genuchten [64] or Brooks-Corey [63], Purcell [66], the LET model [61] or a model proposed by Kjosavik et al. [60]. For the sake of brevity, we focus on the Brooks-Corey (BC) and LET models in this study. The BC relative permeabilities can be derived by integrating the capillary pressure over the capillary tubes [34,60]. After the integration of the BC capillary pressure, one can obtain relative permeabilities for the water-wet system and for a hydrophobic system (see [60]), here a w , a n , m w , and m n are data fitting parameters in which the subscripts w and n indicate the wetting and non-wetting surfaces. Particularly the m's are known to be tortuosity exponents. In 2005, Lomeland et al. [61] have proposed relative permeability models with three parameters L, E, and T for two-phase flow system. Their correlation models read as: where L α , E α , and T α are data fitting empirical parameters. A detailed description and explanation of the parameters can be found in [61]. The important characteristics of the LET model (7) is its flexibility to predict the relative permeability curves for any type of wettability conditions. We note that 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 in order to evaluate the behaviors of ω α (·), a(·), and b(·) along exposure time to the WA agent. The WA induced k rα -S data can be measured from laboratory experiments. However, this approach is expensive in terms of time. Thus, we follow a theoretical approach in which we simulate timedependent k rα -S data from a pore-scale model.

Pore-scale model description
We employed a triangular bundle-of-tubes to represent the pore-scale model. Note that one can also use a simpler pore-scale model, i.e., a cylindrical bundle-of-tubes model, to characterize the impact of WA on k rα -S relations. However, polygonal pores are advanced in the way that they can represent physical processes such as the establishment of mixed wettability within a single pore. This may lead to different fluid distributions within a pore, and establishment of non-wetting fluid layers in the 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 depicted in wetting (left with pressure P res. l ) phase reservoirs. Once the fluid movement is initiated in the tubes, 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 wetting phase, and change of fluids along the tube respectively, whereas α is half angle and θ m is fluid-fluid CA. Detailed calculation of these areas and the angle β m is discussed below.
Let the boundary pressures difference be defined as and the tubes in the bundle are filled with the wetting phase initially. To displace the wetting phase fluid in the m th tube, the pressure drop has to exceed the local entry pressure [22] ∆P > P c,m .
If condition (9) is satisfied, the non-wetting fluid starts to displace the wetting phase, and the volumetric flow rate in the m th can be approximated by the Lucas-Washburn flow model [68], where, µ nw and µ w are non-wetting and wetting fluid viscosities, respectively, the superscript int stands for fluid-fluid interface, q m = dx int m /dt is the interface velocity, R m is the m th tube inscribed radius, and θ m is the fluid-fluid contact angle at pore m. Here, θ m is a general contact angle representation and can be specifically defined as θ r,m if the interface is receding, θ a,m if the interface is advancing, and θ h,m if the interface hinges in the corner of the pores. The G m in Eq. (10) represents the conductance of the fluids and we follow the work of [67] to pre-compute the conductance. The interface is assumed to be trapped when it reaches the outlet of the tube, thus q m = 0 when the interface reaches the boundaries.
The wettability change and/or gradient in the polygonal pores may create distinct fluid configurations, 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 surface in configurations B through F is to show that part of the surface is exposed to a WA agent at some point and experiences a wettability change. Configuration B occurs when the wettability is altered up to a CA value satisfying θ < π 2 − α, where α = π 6 is the corner half-angle and configuration C otherwise. Configuration F may occur when the non-wetting fluid invades configuration E. This is the case when the receding CA is sufficiently smaller than the previous advancing contact angle. Figure 2: Fluid configurations for primary drainage, imbibition, and secondary drainage (water in light color and CO 2 in a dark color). We adopted the configurations from [72]. The bold lines along the sides indicate altered wettability.
Below, relevant aspects of the polygonal pores for the relative permeabilities measurement will be stressed, which includes entry pressure and the areas covered by the two fluids. The capillary pressure at the pore-level is given by where p n and p w are the non-wetting and wetting phase pressures respectively, r m is radius of arc meniscus (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 the system [69]. For isothermal, constant total volume, constant chemical potentials, and incompressible system, 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 ), 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] the entry 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] to calculate 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 calculated in terms of r m and θ m . To calculate these quantities, we numbered the fluid-fluid interfaces in order from the apex if there exist more than one interface in the corner, and we apply the indicator notation 1, if interface k separates bulk nonwetting and corner water −1, if interface k separates bulk water and corner nonwetting.
The bulk cross-sectional area A k b,m in each tube in the bundle is defined as, where Though we used a general notation for CA θ k m above, θ k m may be replaced by, θ r,m and θ a,m when the interface recedes and advances respectively, and θ k h,m if the interface, separating the bulk and corner fluids, is hinging. The hinging contact angle (if exists) changes with the entry pressure P c,m according to: The fluid area that occupies the corner regions can be estimated by where θ m is an argument to determine the appropriate area. The corner surface covered by water in configuration C is obtained by A c,m (θ 1 h,m ), whereas the non-wetting fluid layer in configuration h,m ). The same approach can be applied if there exist many 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 ≤ π 2 − α is satisfied during the drainage displacement. Otherwise, the non-wetting phase occupies the cross sectional area of the 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 > π 2 + α is satisfied. This implies that a dynamic change of CA can also determine the fluid distribution in a single pore.

Pore-scale time-dependent wettability model
Above, we observed that the entry pressure in Eq. (12), fluid distributions in Eqs. (14)- (17), and a conductance in Eq. (10) are wettability evolution dependent. In this section, we introduce a WA mechanism at the pore level to examine its impact on the saturation distribution. We recall that there 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 assumption that: • The CA of the pore surface is altered through exposure time to the WA agent and the alteration is permanent. That means the wettability is not restored to the original wetting condition when the WA agent is displaced by the other fluid unless the displacing fluid has a composition that gradually restores the original wetting condition.
• The WA becomes quasi-static in time if the WA agent is removed from the pore before the final wetting-state is reached. If the agent is reintroduced at some later point, alteration continues until the final state.
According to the assumptions above, the bulk surface area A b,m is supposed to be altered dynamically in time, whereas the corner surface area A c,m may keep the initial condition. Thus, we introduce a functional form to describe a WA mechanism for any arbitrary tube m as where ∆Θ = θ f m − θ i m , θ f m , θ i m 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 so that ϕ 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 very limited. Furthermore, WA is a complex process, where surface free energies, surface mineralogy, fluid composition and exposure time interact. However, adsorption of the WA agent onto the surface area is a natural process in CA change. Such adsorption type wettability evolution is observed for CA measurements [48,50,49,75,77,76]. These all give an insight to model ϕ in Eq. (18) according to the adsorption of the WA agent and can be given as follows where C is a non-dimensional parameter that controls the speed and extent of alteration from initialwet to final-wet system. The derivation of Eq. (19) can be found in our previous work [65]. The variable χ m is a measure of exposure time and is defined as where T is a pre-specified characteristic time, x int m is the fluid-fluid interface position along the tube length L, A b,m is the pore surface area that covered by the non-wetting fluid, and V p,m is the pore volume. Without losses of generality, the characteristic time T is set to be the time for one complete drainage displacement under static initial wetting condition, which can be pre-computed from Eq. (10).
For a given interface position x int m , one can determine the required time to reach the specified interface position from Equation (10). The obtained time is used in Eq. (20) to calculate the exposure history 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 be different for different pores depending on the local saturation history. This process would gives rise to a non-uniform wetting condition across the bundle until all pores have reached the final wetting-state. Moreover, the bulk surface area A b,m is the only surface that undergoes WA, and the corner 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 because the time to drain is assumed to be fast compared to the exposure time for WA to occur.

Simulation approach
The wettability dynamics described in Eq. (18) are coupled into a triangular bundle-of-tubes model to simulate relative permeability curves according to Algorithm 1. Here, the relative permeability for phase α is calculated as 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 to if ∆P > P c then if ∆P < P c then 21: Calculate S nw and χ

23:
Update θ m from Eqs. (18) and (20) 24: end if 25: Increase ∆P 26: 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 for each tube within a few numbers of drainage-imbibition cycles. When we reduce the increment of each ∆P , the flow rate becomes very slow. This imposes a prolonged exposure to the WA agent and thus χ m grows and eventually results in a large change in CA. For example, each ∆P increment is reduced by three order of magnitude for the last cycle compared to the first cycles to complete the alteration process.
The relative permeabilities-saturation "data points" are obtained in each drainage-imbibition cycle, and the obtained data is presented in the following Section. The generated k rα -S curves are used 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 relation between these parameters and changes in the pore-scale WA model parameter C is studied and examined.

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 wetting reservoirs. Each tube in the bundle is assigned a different radius R, with the radii drawn from a truncated two-parameter Weibull distribution [67] 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 relative permeability curves. In the following section, we present and discuss the simulated relative permeability and related results.

End wetting-state relative permeability
In section 2, we point-out that end wettingstate relative permeabilities are the foundation to characterize 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 distribution and fluid properties in Table 1. Thus, we simulated static k rα -S data by fixing the wettability at pre-specified initial θ i m and final θ f m values in each tube, see Table 1. The simulated curves are plotted 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 fitted parameters for the correlation models are found in Table 2. Both the BC correlations and LET models give an excellent match to the simulated phase relative permeability curves under static conditions.
As a complement, we have done experiments on the sensitivity of model parameters (not shown here) for different pore-size distributions by setting the wetting condition to be static. In this experiment, we found that the parameters L n and T w in the LET model can be unity for any poresize 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 are sensitive 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 model   to a two-parameter (but in each phase) model and we call it reduced LET model which read as The model in Eq. (23) is a two parameter model with only one parameter that varies along the wettability change. This makes the reduced LET model more reliable than its counter part, i.e., the BC model. Note that the parameter λ can also vary with wettability. In this case, a better match with the simulated relative permeability data can be obtained.

Simulated relative permeabilities
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 CA model (18). These data are generated with a pore-scale parameter C = 10 × 10 −5 . Note that the CA change may be halted (temporarily) in the pores if the displacement is to configuration D after imbibition. In this case, the drainage curve may follow the previous imbibition path. However, the imbibition curve may show a deviation from the previous drainage curve if wettability is altered sufficiently. According to the CA distribution in Fig 5 and k rα -S curves in Fig 4, the wettability that ranges from strongly to weakly water-wet does not affect the pore filling and draining orders of the poresizes. As a consequence, the k rα -S relations for drainage/imbibition displacements follow the same saturation path of the initial-wet condition. This implies that a WA induced k rα -S hysteresis may not occur in a straight bundle-of-tubes model for a wettability range of 0 • ≤ θ m ≤ 60 • . Similar k rα -S relations are reported for a pore-network model with a wettability range of θ m = 0 • to 45 • [36]. This is 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 occur in 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-saturation path deviation as observed in Fig 4. We observe that the wetting and non-wetting phase relative permeabilities steadily increase and decrease (see Fig 4) respectively in each subsequent drainage-imbibition displacements when the wettability evolves from hydrophilic to hydrophobic conditions. This occurs because the relatively larger pores start allowing the water (originally wetting phase) to imbibe through before the smaller pores when the medium is altered to intermediate/hydrophobic system. On the contrary, the originally 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 relative permeabilities are observed. However, any additional drainage-imbibition cycle would follow along the 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 dynamic 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. This is one of the challenging features of WA during the quantification of the dynamics in relative permeabilities. To eliminate the hysteresis observed in the k rα -S relation, we projected the simulated relative permeability data onto the temporal domain χ in Fig 7. The temporal domain, χ is the measure of the exposure history in averaged sense which defined as: Unlike k rα -S, k rα -χ is uniquely related but non-monotonically i.e., it raises to one and decreases to zero in time along each drainage-imbibition cycle. The WA process also affects the corner fluid distribution in each drainage/imbibition displacement. Fig 7 shows the evolution of (average) corner wetting and non-wetting saturations. The wetting phase saturation decreases through time while the layer saturation grows. This is because wettability is altered from water-wet to hydrophobic and thus, the originally non-wetting fluid prefers to be in the corners. However, when (all) pores become more hydrophobic, the corner water increases and bulges-out during imbibition. This process may result in a layer collapse and has shown in the fifth imbibition, where the layer saturation increases during larger pores were imbibed and decrease when smaller pores imbibed over exposure time. Here, we have checked the establishment of corner fluids during the calculation of k rα -S relations in each drainage-imbibition displacements.

Dynamic relative permeability model development
The interpolation approach, similar to Eq. (4), was applied successfully to capture time-dependent WA mechanisms in the capillary pressure curves [65]. Thus, it is important to test the potential of the 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, we observe 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. However, the ω α -S curves in Fig 8a have piece-wise functional forms i.e., zero and Langmuir-type function of phase saturation, that are altered with exposure time. This imposes another challenge to find a smooth model to predict the curves along with the saturation. But, one can design a piece-wise function to correlate the ω α − S data. From Fig 8a, we observe that the starting point of the Langmuir functional form is exposure time-dependent along the saturation path. Thus, the ω α − S can be represented as, Parameter Vale Parameter Value where S * w and β are time-dependent. The variable S * w is used to transform the starting point of the Langmuir part ω α along the saturation to zero, and β is used to determine the curvature of the curve. From Fig 8a, we observe that the parameters S * w and β are decreasing functions of exposure time, χ.
We matched the designed model in Eq. (25) with the dynamic coefficient data, ω α − S w , to describe the relations S * w − χ and β − χ. The obtained relations have the form, where a i , b i and c i for i = 1, 2 are dimensionless fitting parameters. These parameters are estimated and given in Table 3 for this particular simulation. The models in Eqs. (25)-(26) are then substituted back into the interpolation model (4) to give the dynamic relative permeability model and is read as, The dynamic model (27) is compared with the simulated relative permeability in Fig 9. According to Fig 9, the designed model predicts beyond the initial wetting state at joint-point of the initial wetting-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 parameters to be calibrated.  Figure 9: Comparison of the interpolation model (27) with the simulated (wetting (a) and nonwetting (b)) relative permeabilities.
The second approach discussed in Section 2 relies on establishing a relation between the model parameters 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 model parameters are dependent on the wetting condition of the porous domain in addition to the pore-size distribution. 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 that only E α is sensitive for wettability change in the reduced LET model, whereas both parameters are sensitive in the BC model (see Table 2). This implies that the BC model is more expensive than the reduced LET model to upscale the effect of WA on relative permeabilities. Thus, we choose the reduced LET model to represent the dynamic relative permeability curves. Here, we rearrange the reduced LET model as where λ and E n are determined from the initial wetting-state correlation. Here, L n is designed to change with CA. Thus, we coupled the reduced LET model (28) with curve fitting tool in MATLAB and matched with the k rα -S data to study the functional dependencies between the L n and χ. The obtained relation is linear and has the form L n (χ, E n ) = a n χ + E n , where E n is as given in Table 2 and a n is dynamic fitting parameter for wetting and non-wetting phase relative permeabilities and determines the slope of the relative permeabilities along exposure time. For this particular simulation this parameter is estimated to be a n = 3.57 for all dynamic drainage-imbibition cycles reported above. Now the dynamic term L n in Eq. (29) can be substituted into the reduced LET model (28) to give the dynamic relative permeabilities: where the parameters E n and λ are pre-determined from the initial wetting-state correlation. The proposed dynamic relative permeability models in Eq. (30) and the simulated k rα -S data are com-pared in Fig 10. From Fig 10, we observe that the proposed dynamic model (30) 30)), the reduced LET model is more efficient and simpler to implement in the Darcy flow than the piece-wise interpolation model. Because, the reduced LET model is smooth along time and saturation. We also note that the number of parameters in the modified 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 multiple parameters that change in each cycle (or hysteresis models). Thus the analysis below will concern only on the modified LET model.

The modified LET model sensitivity to pore-scale model parameter
In this section, we will investigate the response of the upscaled model parameter a n to the change in the CA model parameter C in Eq. (18). The parameter C controls the extent of WA at the pore-level, whereas a n determines the WA induced dynamics in the relative permeabilities. We have simulated 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 by considering different values of C. Then, we correlate the estimated parameter values of a n with the chosen values of C, which can be read as a n = P 1 C + P 2 , where P 1 and P 2 are fitting parameters. The correlation result is plotted in Fig 11, where the parameters are estimated to be P 1 = −18420 and P 2 = 6. Therefore, we can predict the upscaled dynamics of the relative permeabilities directly from the pore-scale WA process. Data Model 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) to complete 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 estimated from the calibration of CA model (18) with experimental data. Validating the underlying CA change model is beyond the scope of this paper. Rather, we consider the CA model as a reasonable basis to perform and analyze the upscaling process.

Applicability of the modified LET model to arbitrary saturation history
The saturation path that used to generate the relative permeability data in Fig 12 can be considered as one of the arbitrary paths in the domain S w × χ. This saturation path is dependent on the history of the exposure time to the WA agent (degree of WA change) and the reversal-point of saturation. This implies that the relative permeability-saturation dynamics may behave differently if one chooses a different saturation path that entails a prolonged exposure time for a fixed saturation profile and/or flow reversal at intermediate saturation.
Here, we simulate as many as possible k rα -S curves that involve different reversal-points and exposure history within the S w ×χ. Note that we use the pore-scale model parameter C = 10×10 −5 . The simulated data is plotted in Fig 12a and 12b for phase relative permeabilities. These arbitrary curves are used to test the potential of the modified LET model in Eq. (30). To do so, we apply the calibrated dynamic relative permeability model (30) to generate the k rα -S-χ surface. The absolute difference between the simulated data and the surface generated by the calibrated model is depicted in Fig 12c. According to the results in Fig 12, we can justify that a single saturation history path is sufficient to calibrate a dynamic model that can be applied to any saturation-time path.

Discussion
In contrast to the dynamic capillary pressure model presented in [65], the interpolation-based approach is poorly correlated with the simulated relative permeability curves. However, further in- 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 the dynamic model with the simulated data vestigation may be needed to re-evaluate the potential of capturing the WA process in the relative permeabilities based on the interpolation approach. We also examined the BC model to upscale the WA induced dynamics in the relative permeabilities though we do not report it here fully. For the sake of brevity, we have highlighted only the response of the BC model to the wettability change in Section 3.3 for end wetting conditions. In general, other models that involve more than two parameters (sensitive to CA change) could be calibrated with reasonable accuracy. But, these parameters need to be adjusted in each drainage-imbibition displacements. This complicates the modeling process and the resulting model may involve many parameters that may impose an extra challenge to analyze 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 existing static-wettability model i.e., the LET model. The original LET model consists of three parameters in each phase that need to be adjusted differently for different wetting conditions. First, we did numerical 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 phase relative permeability models are constant for any pore-size distribution. From this, we reduced the LET model to a model (the reduced model in Eq. (23)) that involves two parameters in each phase.
We 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 relation between E α and the exposure time χ which allows us to come up with a single-valued dynamic relative permeability model that represents any arbitrary drainage-imbibition cycle. We note that the model involves 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 these values, the parameter a n is the only parameter that controls the WA induced dynamics in the relative permeabilities for both phases.
The proposed model is at the Darcy scale, that allows for a change in relative permeability as a function of averaged variables such as saturation (S w ) and exposure time to a WA agent (χ). This implies that the relative permeability in a grid block that is exposed to the WA agent may change over time even for constant saturation profile. However, if the grid block is not exposed to the WA agent, χ is zero for this particular grid block. In this case, the dynamic model predicts the initial wetting-state curve. However, the developed model may continue further after the end wetting-state curve has attained which is in contrast to the interpolation-type model, where prolonged exposure does not contribute to the dynamics once the final wetting curve is met, see [65]. Thus it is important to propose a strategy to ensure that the relative permeability dynamics do not to cross the end-state curve. 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 attained when a n χ + E n = E f n is satisfied. From this, we can estimate the exposure time needed to reach the final wetting state (say χ max ) such that the relative permeability is represented by the end wetting-state curve. After knowing this we can set the dynamic variable as 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 end of the curve for the wetting and non-wetting phases respectively.
Previous studies [25,3,27,28] represent the impact of instantaneous WA on relative permeabilities by an interpolation model which matched directly to core-scale data in a heuristic manner. This study revealed that the interpolation model is not the best approach to upscale the pore-scale WA process. Rather, we have shown the potential of a modified LET model to capture the underlying WA process at the pore-scale represented by a triangular bundle-of-tubes. The proposed models are smooth and simple to use for practical applications. Most importantly, the models are designed to eliminate the relative permeabilities hysteresis induced by CA change during drainage and imbibition displacements. Similar to the developments in [65], we have quantified the link between the 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 results, 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-level can be used to predict the macroscale dynamics without performing pore-scale simulations. This is an important and valuable generalization for making use of experimental data to inform core-scale relative permeability-saturation relations.

Conclusion
In this paper, we developed a dynamic relative permeability model that includes the pore-scale underpinnings of WA in the relative permeability-saturation relationships at the Darcy scale. We found that the developed model (i.e., the modified LET model in Eq (28)) is simple to use and can predict WA induced changes in the relative permeabilities. The modified LET model shows a good agreement with the simulated relative permeability data. Furthermore, this model is independent of the saturation-time paths generated by any drainage-imbibition cycles. More importantly, the WA dynamics in the relative permeabilities is controlled by a single-valued parameter that has a clear relationship with the time-dependent CA change model parameter.