Dynamic wettability alteration in immiscible two-phase flow in porous media: Effect on transport properties and critical slowing down

The change in contact angles due to the injection of low salinity water or any other wettability altering agent in an oil-rich porous medium is modeled by a network model of disordered pores transporting two immiscible fluids. We introduce a dynamic wettability altering mechanism, where the time dependent wetting property of each pore is determined by the cumulative flow of water through it. Simulations are performed to reach steady-state for different possible alterations in the wetting angle ($\theta$). We find that deviation from oil-wet conditions re-mobilizes the stuck clusters and increases the oil fractional flow. However, the rate of increase in the fractional flow depends strongly on $\theta$ and as $\theta\to 90^\circ$, a critical angle, the system shows critical slowing down which is characterized by two dynamic critical exponents.


I. INTRODUCTION
The world's primary energy demand is predicted to increase by one-third between 2011 and 2035, where 82% of it comes from fossil fuels [1]. In this scenario, the fact that some 20 to 60 percent of the oil remains unrecovered in a reservoir after the production is declared unprofitable, is a challenge of increasing importance [2]. The main reason for this loss is the formation of oil clusters trapped in water and held in place by capillary forces, which in turn are controlled by the wetting properties of the reservoir fluids with respect to the matrix rock. The production from complex oil reserves that today are considered immobile or too slow compared to the cost is therefore an important area of research. In this context, the role of formation wettability is a focus area within the field of Enhanced Oil Recovery (EOR) [3]. Different reservoir rocks have widely different wetting characteristics [4]. Wettability may vary at the pore level from strongly oil wet through intermediate wetting to strongly water wet. Carbonate reservoirs contain more than half of the world's conventional oil reserves, but the oil recovery factor is very low compared to sandstone reservoirs [5]. This is due to the complex structure, formation heterogeneity and more chemically active wettability characteristics of the carbonate reservoirs, which leads to uncertainty in the fluid flow and oil recovery [6]. Sandstone is strongly water wet before oil migrates from a source rock into the reservoir. When oil enters a pore, it displaces the water which leaves behind a water film sandwiched between the oil and rock surface. This happens as a result of balancing van der Waals and electric double layer forces, capillary pressure and grain curva- * vegard.flovik@ntnu.no † santanu.sinha@ntnu.no ‡ alex.hansen@ntnu.no ture [7]. A permanent wettability alteration is then believed to take place by adsorption of asphaltenes from the crude oil to the rock, and leads to high but slow recovery through continuous oil films [8,9]. As the oil saturation drops, these films can become discontinuous, leaving immobile oil clusters held in place by capillary forces. Wettability is therefore an important petrophysical property which plays a key role in the fluid transport properties of both conventional and unconventional reservoirs [10], and there is great potential to improve the oil recovery efficiency by altering the wetting properties. Main factors which can alter the pore wettability are; lowering the salinity [11], adding water-soluble surfactants [12,13] or adding oil-soluble organic acids or bases [14]. Increasing the reservoir temperature also increases water-wetness [4,15]. There are some correlations with the wetting behavior to the electrostatic forces between the rock and oil surfaces [16], but there is no consensus on the dominating microscopic mechanism behind the wettability alteration. It is known from laboratory experiments and field tests that a drift from strongly oilwet to water-wet or intermediate-wet conditions can significantly improve the oil recovery efficiency [14]. The amount of change in the wetting angle is a key factor here [17,18] which not only decides the increase in oil flow but also the speed of the process. An improper change in the wetting angle can also make the recovery very slow and not profitable.
Given there is a certain change in the wetting angle due to a brine, the next important factor is the flow pathways in the matrix rock which transports the oil and brine. One cannot expect any change in the wetting angle of a pore if there is no flow of the brine through it. The flow pathways depend on several different factors: the porous network itself, oil saturation, capillary number and also on the present wettability conditions. A change in the wettability will cause a perturbation in the flow distribution of the system. This will in turn again affect the wettability change through the altered flow pathways, causing further changes in the flow distribution. The dynamics of wettability alterations is therefore controlled by a strongly correlated process.
There are some studies of wettability alterations in two-phase flow by equilibrium-based network models [19] for capillary dominated regimes where viscous forces are negligible. Moreover, investigation of the time-scale of dynamics lacks attention in such models. In this article, we present a detailed study of wettability alterations in two-phase flow considering a network model of disordered pores transporting two immiscible fluids where a dynamic wettability alteration mechanism, correlated with the flow-pathways, is implemented. We will focus on the transport properties due to the change in the wettability as well as on the time scale of the dynamics.
We study in the following the effect of wettability changes on immiscible two-phase flow based on a network model [21][22][23]. In section II, we present the model and how we adapt it to incorporate the dynamic wettability changes. In Section III, we present our results. Initially, we let the two phases settle into a steady state where the averages of the macroscopic flow parameters no longer evolve. At some point, we then introduce the wettability altering agent, so that it starts changing the wetting angle. The wetting angle alteration depends on the cumulative volume of the wettability altering fluid that has flowed past a given pore. This induces transient behavior in the macroscopic flow properties and we measure the time it takes to settle back into a new steady state. We find that there is a critical point at a wetting angle of 90 degrees and we measure its dynamical critical exponents; the exponents are different whether one approaches the critical point from smaller or larger angles. In Section IV we summarize and conclude.

II. MODEL
We model the porous medium by a network of tubes (or links) oriented at 45 degrees relative to the overall flow direction. The links contain volumes contributed from both the pore and the throat, which then intersect at volume-less vertices (or nodes). Any disorder can be introduced in the model by a proper random number distribution for the radius r of each link, and we choose a uniform distribution in the range [0.1l, 0.4l] here, where l is the length of each tube. The network transports two immiscible fluids (we name them as oil and water), one of which is more wetting than the other with respect to the pore surface. The pores are assumed to be in between particles, and the pore shape is thus approximated to be hour-glass shaped, which introduces capillary effects in the system. The model is illustrated in figure 1.
Due to the hour glass shape of the pore, the capillary pressure at a menisci separating the two fluids is not constant, and depends on the position x of the menisci inside the pore. The capillary pressure p c (x) at position x inside the ith pore is then calculated from a modified form of the Young Laplace equation [20,21], where γ is the interfacial tension between the fluids and ϑ i is the wetting angle for that pore. As an interface moves in time, p c (x) changes. The capillary pressure is zero at the two ends (x = 0 and l) and it is maximum at the narrowest part of the pore. It makes the model closer to the dynamics of drainage dominated flow, where the film flow can be neglected. When there are multiple menisci in a pore, the total capillary pressure inside the ith pore is obtained from the vector sum ( i p c (x)) over all the menisci in that pore. The flow is driven by maintaining a constant total flow rate Q throughout the network, which introduces a global pressure drop. The instantaneous local flow rate q i inside the ith link between two nodes with pressures p 1 and p 2 follows the Washburn equation of capillary flow [24], where ∆p i = p 2 − p 1 . k i = r 2 i /8 is the permeability of cylindrical tubes. Any other cross-sectional shape will only lead to an additional overall geometrical factor.
is the volume averaged viscosity of the two phases inside the link, which is a function of the oil saturation s i in that link. Here µ o and µ w are the viscosities of oil and water respectively.
The flow equations for the tube network are solved using a conjugate gradient method [25]. These are the Kirchhoff equations balancing the flow, where the net fluid flux through every node should be zero at each time step, combined with the constitutive equation relating flux and pressure drop across each tube. The system of equations is then integrated in time using an explicit Euler scheme with a discrete time step and all the menisci positions are changed accordingly. Inside the ith tube, all menisci move with a speed determined by q i . Physical processes like bubble snap-off and coalescence are introduced in the model. Due to these, bubbles can be formed or merged, which changes the number of menisci inside a link with time. When a meniscus reaches the end of a tube, new menisci are formed in the neighboring tubes depending upon the flow rates. In each link, a maximum number of menisci is allowed to form. When this number is exceeded, the two nearest menisci are merged, keeping the volume of each fluid conserved. In our simulations, we considered a maximum of four menisci inside one pore which can be tuned depending on the experimental observations. Further details regarding the menisci movement can be found in Reference [22].
The simulations are started with an initial random distribution of two fluids in a pure oil-wet network. Biperiodic boundary conditions are implemented in the system, which effectively makes flow on a torus surface. The flow can therefore go on for infinite time, keeping the saturation constant and the system eventually reaches to a steady state. In the steady state, both drainage and imbibition take place simultaneously and fluid clusters are created, merged and broken into small clusters. One can consider this as the secondary recovery stage. Once the system reaches the steady-state in a oil-wet network, the dynamic wettability alterations are implemented, which may be considered as the tertiary recovery stage or EOR. In the following we discuss this in detail.

Dynamic wettability alteration
We now introduce a dynamic wettability alteration mechanism to simulate any wetting angle change, decided by the oil-brine-rock combination and the distribution of the flow channels in the system. In a previous study [23], a simplified static wettability alteration mechanism was studied, where the alteration probability was considered equal for all pores without any correlation with the flow of brine inside a pore. However, for wettability alterations to occur, the wettability altering agent (e.g. low-salinity water or surfactant) needs to be in contact with the pore walls. Thus, the wettability alteration should follow the fluid flow pathways and any change in the wetting angle inside a pore should depend on the cumulative volumetric flow of brine through that pore. This claim is rather trivial, as one can not expect any wettability change in a pore if the altering agent is not present. This means that if a certain pore is flooded by large amounts of brine, the wetting angle should change more in that pore than the one which had very little water flooded through. This is implemented in the model by measuring the cumulative volumetric flux V i (t) in each individual pore with time t, where t 0 is the time when the injection of low salinity water is initiated, ∆t is the time interval between two simulation steps and (1 − s(t)) is the water saturation. V i (t) is then used to change the wetting angle for each tube continuously, updated at every time step after t = t 0 . The wetting angle ϑ i of the ith pore can change continuously from 180 • to 0 • as V i (t) changes from 0 to ∞. Correspondingly, the cos ϑ i term in Eq. (1) will change from −1 to 1 continuously. This continuous change of the wetting angle with the variation of V i (t) is modelled by a function G i (t) given by, which replaces the cos ϑ i term in Eq. (1). The pre-factor 2/π is a normalization constant to set the range of the function, and the parameter C adjusts the slope during the transition between oil wet and water wet. This leads to the change in the wetting angle as a function of V i (t) as shown in Fig. 2. As our model does not include film flow, the wetting angle should not reach either 0 or 180 degrees. Rather, a starting point of ϑ i ≈ 165 • was used by choosing C = 20 in our simulations as shown in Fig. 2. Next, as a larger pore will need more brine to be flooded in order to have a similar change in the wetting angle than a smaller pore, a threshold value V th i is introduced, which is proportional to the volume of that pore, At V i (t) = V th i , the wetting angle reaches to 90 • in that pore and p c (x) essentially becomes zero. Here η is a proportionality constant which decides how many pore volumes of water is needed to reach V i (t) = V th pore. This parameter can possibly be adjusted against future experimental results, but is considered as a tuning parameter in this study. The expression for the capillary pressure at a menisci from Eq. 1 then takes the form, The maximum amount of wetting angle that can be changed depends on the combination of brine, oil and rock properties [17,18]. We therefore set a cut-off θ in the wetting angle change, such that any pore that has reached to a wetting angle ϑ i = θ, can not be changed further. The model thus includes all the essential ingredients of wettability alteration study -it is a time dependent model where the wettability alteration is correlated with the flow pathways of the brine, and can be used to study any oil-brine-rock combination decided by θ.

III. RESULTS
Simulations are started with a random initial distribution of oil and water in an oil-wet network, where θ ≈ 165 • for all links. First, the oil-wet system is evolved to a steady state before any wettability alteration is started. This will allow us to compare the change in the steady-state fractional flow of oil (F ) with a change in the wetting angle. The oil fractional flow (F ) is defined as the ratio of the oil flow-rate (Q oil ) to the total flow-rate (Q) given by, F = Q oil /Q. The flow rate (Q) is kept constant throughout the simulation, which sets the capillary number Ca = µ eff Q/(γA), where A is the crosssectional area of the network. A network of 40 × 40 links are considered, which is sufficient to be in the asymptotic limit for the range of parameters used [22]. An average over 5 different realizations of the network has been taken for each simulation. As the simulation continues, both drainage and imbibition take place simultaneously due to bi-periodic boundary conditions and the system eventually evolves to a steady state, with a distribution of water and oil clusters in the system. In Fig. 3, F is plotted against the number of pore volumes passed (N ) through the network. As we run the system with constant flow-rate, N is directly proportional to the time t, N = tQ/v where v is the total volume of the network. The initial 200 pore volumes are for an oil-wet network, where it reaches to a steady state with F ≈ 0.235. We then initiate the dynamic wettability alteration which resembles the flow of a wettability altering brine and F starts to drift. Here we run simulations for different values of η, defined in Eq. 5, and the results are plotted in different colors. θ = 0 • in these simulations, which means any pore can change to pure water-wet depending upon the flow of brine through it. One can see that F approaches to a new steady-state with F ≈ 0.308 due to the wettability alteration. Here different values of η make the simulation faster or slower to reach the same steady state, where a higher value of η makes the system slower to reach the steady state. The initialization of steady state measured in terms of the pore-volumes (N ss ) is defined as the instant when the average fractional flow stops changing with time anymore and essentially stays within its fluctuation. In the inset of Fig. 3, N ss v is plotted with η and a simple linear dependency is observed. In order to save computational time, we therefore use η = 10 in all our following simulations.
How the two fluids and the local flow-rates are distributed in the network in the two steady-states before and after the wettability alteration are shown in Fig.  4. The network size is 64 × 64 links here with an oilsaturation S = 0.3 and the capillary number Ca = 10 −2 . All the links are hour-glass shaped in the actual simulation with disorder in radii, but shown as a regular network for simplicity in drawing. The upper row shows the distribution of oil bubbles drawn in black. The left column (a) shows the steady state in a oil-wet network and the right column (b) shows the steady-state after the wettability alteration is initiated with maximum possible wetting angle change θ = 0 • . A closer look in these bubble distributions shows more clustered oil bubbles in (a) than in (b) where they are more fragmented. A more interesting picture can be seen in the local flow-rate distribution in the bottom row, which shows a more dynamic scenario. The left (c) and right (d) figures are for the same time-steps before and after wettability alteration as in (a) and (b) respectively. Here the local flow-rates in each pore, normalized in between 0 to 1, are shown in gray scale. Interestingly, in the oil-wet system (c), the flow is dominated in a few channels (black lines) where the flow-rates are orders of magnitude higher than the rest of the system. Other than those channels, the system has negligible flow, indicated by white patches which means the fluids are effectively stuck in all those areas. This situation happens when the difference in the saturation of the two fluids is large, where the phase with higher saturation (water here) tries to percolate in paths dominated by a single phase with less number of interfaces. This is not favorable in oil-recovery, as it leaves immobile fluid in the reservoir. In the flow distribution after the wettability alteration (d), the flow is more homogeneous and distributed over the whole system, indicating higher mobility of the fluids. However, one should remember that when the wettability alteration is started in a sys- tem shown in (c), the wettability alteration starts taking place only in those pores with active flow. But then it perturbs the flow-field and starts new flow paths and eventually the system drifts towards a more homogeneous flow with time, as shown in (d).
We now present the results when the wetting angle of any pore can change all the way down to zero degree (θ = 0 • ). In Fig. 5 the steady-state oil fractional-flow in an initial oil-wet system (F ) is compared with that in the steady-state after wettability alteration (F ′ ). Results are plotted as a function of S for two different capillary numbers, (a) Ca = 10 −1 and (b) Ca = 10 −2 . The diagonal line corresponds to F = S. A miscible fluid mixture would follow this line, and it is interesting to use this as a reference to compare how the data points lie above or below this line. If both the fluids flow equally in the system then F will be equal to S. But in Fig. 5, F stays below the F = S line for low oil saturation values in the oil-wet system, which means that the flow of oil is less than that of water and there are stuck region of oil. A significant increase in F can be observed here due to the wettability alteration for the full range of oil-saturation. Moreover, increase in F is higher for the lower capillary number, indicating that wettability alteration is very significant in the case of oil recovery, as Ca can go as low as 10 −6 in the reservoir pores. Fractional flow also obeys the symmetry relation F ′ (S) = 1 − F (1 − S) [23] which implies that, if the wetting angle of any pore is allowed to change all way down to zero degree (θ = 0 • ), the system will eventually become pure water-wet with time.
As noted earlier, the maximum change in wetting angle for a system, depends on the properties of the reservoir rock, crude oil and brine, and also on the temperature. Existing wettability alteration procedures generally turns the oil-wet system into intermediate wet, rather than to pure water-wet. Some examples of the change in the wetting angle for different rock materials and brine can be found in [17,18]. In our simulation this is taken care of by the parameter θ, which decides the maximum change in the wetting angle ϑ i for any pore. One should remember that, we are not forcibly changing the wetting angles ϑ i to θ, rather the change in ϑ i is decided independently for individual pores by the amount of brine passed through it (by Eqs. 3 to 6), and there is a max- imum allowed change in any ϑ i . As before, simulations are started with a pure oil-wet system to reach a steady state and then wettability alteration is started and simulation continues until the system reaches to a steady state again. Independent simulations have been performed for different values of θ. The proportionate change in the oil fractional-flow due this wettability change from the oil-wet system, ∆F/F = (F ′ − F )/F is measured for different simulations and plotted in Fig. 6 against θ.
There are a few things to notice. First, as one can immediately see, fractional flow increases with the decrease of oil-wetness, θ → 0 • . The maximum increase in F is higher for lower Ca, about 86% for Ca = 10 −2 and about 32% for Ca = 10 −1 . This is because the change in wetting angle affects the capillary pressures at the menisci, so the change in F is larger when the capillary forces are higher. Secondly, the major change in F happens in the intermediate wetting regime, upto θ ≈ 60 • , and then it becomes almost flat afterwards. Moreover, this increase in F is more rapid for lower value of Ca. All these facts points towards an optimal range of wetting angle change to increase the oil flow. This is an important observation for practical reasons, as it is not necessary to change the wetting angle further. Thirdly, there is a discontinuity in the curve exactly at θ = 90 • , as we will discuss later. Increase in the oil fractional-flow with the increase in the water-wetness may seem to be obvious and reported by many experiments and field tests. But, the most important concern for the oil industry is the rate of increase, or the time required to achieve a significant increase in the oil production. If the increment in oil flow is very slow compared to the cost of the process, then the oil recovery is declared as not profitable and the reservoir may be considered to be abandoned. As per our knowledge, there are very few systematic studies reported in the literature predicting the time scale to change the oil flow due to the wettability change by two-phase flow of brine and oil in a porous media. We observe that, due to the correlations between the flow paths and the wetting angle change, the time scale of the process varies dramatically with θ. This is illustrated in Fig. 7, where F is plotted as a function of pore volumes (N ) of fluid passed through the system. The initial 400 pore volumes are for an oil-wet system and then results of few different simulations with θ = 85, 88, 89, 90, 91, 92 and 94 degrees are plotted. Interestingly, the rate at which the system reaches a new steady-state, varies significantly depending on the value of θ. For example, after the wettability alteration is started, it needs to flow less than 100 pore volumes to reach the new steady state for θ = 94 • whereas more than 300 pore volumes are needed to reach a steady state for θ = 91 • . Therefore, even if the final steady-state fractional flow is higher for θ = 91 • than for θ = 94 • , it might not be profitable to alter the wetting angles to 91 • because of the slow increase in F . In general, the process becomes slower and slower as θ → 90 • from both sides. Such kind of slow increase in oil recovery as θ → 90 • is also observed in experiments [17,18]. This slowing down of the process is an combined effect of two factors. First, the fact that wettability only can change in the pores where there is flow of brine and the second is the value of θ. All the pores were initially oilwet (ϑ i ≈ 165 • ) and when it reaches the steady state, the flow finds the high mobility pathways depending on the mobility factor of the pores and the capillary pressures at the menisci. When the wettability alteration is started, the wetting angles of the existing flow pathways start decreasing. As a result, capillary pressures at menisci in those channels first decreases as ϑ i → 90 • and then it increase afterwards as ϑ i → 0 • . This creates a perturbation in the global pressure field and correspondingly viscous pressures start changing with time which changes the flow field. However, capillary pressures at the zero-flow regimes are now higher than the high-flow regimes which makes it difficult to invade the zero-flow regimes causing a slower change in the flow field as ϑ i approaches 90 • . An interesting feature is observed exactly at θ = 90 • , where the average fractional flow does not change at all after the wettability change. At exactly θ = 90 • , capillary pressures in all the pores in the existing flow pathways essentially become zero, making them the lowest resistive channels with zero capillary barriers. As a result, the fluids keep flowing in the existing channels forever and the system stays in the same steadystate. The time taken to reach another new steady-state is therefore infinite at θ = 90 • and it therefore is a critical point for the system.
We now measure the steady-state initialization time τ , defined as the moment when the average of the fractional flow stops changing with time and becomes horizontal with the x axis. This is shown In Fig. 7, where the initialization of steady states is marked by crosses on the respective plots. As the simulations are performed with constant Q, τ is proportional to the fluid volume passed through the system and therefore we measure τ in the units of N . τ for different simulations with different values of the maximum wetting angle (θ) is plotted in Fig.  8 for (a) Ca = 10 −1 with S = 0.3 and (b) Ca = 10 −2 with S = 0.4. One can see that τ diverges rapidly as θ approaches θ c = 90 • from both sides, θ > 90 • and < 90 • . This divergence of the steady-state time τ as θ → θ c indicates the critical slowing down of the dynamics, which is a characteristics of critical phenomena. The critical slowing down is the outcome of the divergence of correlations at the critical point and can be characterized by a dynamic critical exponent z defined as τ ∼ ξ z , where ξ is the correlation length [26]. As θ → θ c , the correlation length ξ diverges as |θ−θ c | −ν and therefore τ ∼ |θ−θ c | −α , where α = zν. In Fig. 9, τ is plotted versus |θ−θ c | in loglog scale which gives two different slopes for θ > θ c and θ < θ c . We then find the value of the dynamic exponents α as α = 1.2 ± 0.1 for θ > θ c and α = 1.0 ± 0.1 for θ < θ c . However, they are the same within error bar for different capillary numbers and saturations (Fig. 9). The value of the dynamic critical exponents depend on the underlying dynamics and on the model [27]. In this case, wettability alteration was started from an oil-wet system with ϑ i = 165 • for all the pores. So for the simulations with θ < 90 • , the wetting angles cross the critical point (90 • ) when the capillary forces change directions. This might cause the system to mobilize the clusters somewhat faster than for θ > 90 • when the capillary forces does not change any direction. As a result, α becomes smaller for θ < 90 • than for θ > 90 • .

IV. CONCLUSIONS
In this article we have presented a detailed computational study of wettability alterations in two-phase flow in porous media, where the change in the wetting angle in a pore is controlled by the volumetric flow of the altering agent through it. When the wetting angles are allowed to alter towards water-wetness, the stuck oil clusters start to mobilize and oil-fractional flow increases. However, due to the correlations in the wetting angle change with the flow pathways, the time-scale of the dynamics strongly depends on the maximum allowed change in the wetting angle. We find that, as the final wetting angle is chosen closer to 90 • , the system shows a critical slowing down in the dynamics. This critical slowing down is characterized by two dynamic critical exponents. The critical point we are dealing with is an equilibrium critical point as the system is in steady state. The dynamical critical exponents measure how long it takes to go from one steady state to a new one. To our knowledge, this is the first example of there being different values for the exponents on either side of the critical point. Our findings are in agreement with experimental observations reported in literature, and are extremely important for application purposes like oil recovery, where the time scale of the process is a key issue.