Digital rock workflow to calculate wettability distribution in a reservoir rock

. Wettability has a strong influence on multi-phase flow behavior through reservoir rock. Reservoir rocks tend to have spatially varying wettability. Prior to contact with oil, rocks are almost always naturally water-wet. As oil invades the pore-space over geologic time, the initial water-wet state may be altered in certain locations due to adhesion of substances within the oil phase to the grains. Mechanisms of wettability alteration depend on various properties such as pressure, temperature, mineral chemistry, surface roughness and fluid composition. In this study wettability alteration in a reservoir rock is studied through direct simulation using multiphase Lattice Boltzmann method where the computational grid is constructed from segmented micro-CT images of the rock sample. The pore-grain interface is defined by a triangulated surface mesh for accurate fluxes near boundary and local curvature calculation. A capillary pressure drainage simulation is conducted in a water-wet Berea sandstone sample initially filled with water. When oil invades the pore space as the capillary pressure is increased, a fraction of the pore-grain surface is altered towards an oil-wet condition, as determined by a novel wettability alteration process. This process calculates local curvature at every surface element of the rock, obtains local capillary pressure from the simulation and assumes a disjoining pressure to determine water-film breakage at every location of the pore-grain surface. As a result, a spatially varying rock wettability is created. Using this new wettability distribution, the simulation is continued to allow the fluid phases to redistribute accordingly. The process is iteratively carried out until both fluid saturation and wettability distribution converged at a given applied capillary pressure. Afterwards, the pressure is ramped up to the next stage and the process is repeated again. It has been found that the wettability alteration is a slow dynamic process where the non-wetting phase can gradually invade finer pore space as the surrounding grain wettability is altered. In this study, it has also been found that wettability alteration of the reservoir rock produces lower connate water saturation during primary drainage compared to the simulation results without alteration. The resulting spatially varying wettability distribution from primary drainage is used for a subsequent water flooding simulation to calculate water-oil relative permeability curves. The methodology presented in this work can be leveraged to better understand and predict an improved mixed wetting conditions found in the reservoir rocks which is needed for more accurate displacement tests such as relative permeability simulations.


Introduction
Two-phase relative permeability describes how easily one fluid can be moved through a porous medium in the presence of another fluid.It is an important characteristic of hydrocarbon reservoir rocks and a crucial input to oil and gas reservoir simulation modeling.Relative permeability, and multiphase flow through porous media in general, is dependent on various characteristics of the fluid-fluid-rock system, including rock surface properties, physical properties of each fluid, and flow conditions.They can be expressed as dimensionless numbers such as viscosity and density ratio, capillary number, wettability and geometric state function [1].Wettability is a measure of a fluid's tendency to adhere to a solid surface in the presence of other immiscible fluids.Macroscopically, it can be measured by the contact angle at the three-phase contact line and classified as water-wet (minimum 0°, maximum 60 to 75°), neutral-wet (minimum 60 to 75°, maximum 105 to 120°) and oil-wet (minimum 105 to 120°, maximum 180°) in the petroleum engineering context.Microscopically, it can be related to the molecular energy difference between the energy of the combined system and the combined bulk energies of the individual components.The physics of multiphase flow at pore scale is very different for water-wet and oil-wet rocks and its impact can be found in capillary pressure or relative permeability curves [2][3][4][5][6], both critical inputs for reservoir simulation forecasting.Thus, an accurate wettability condition is an important input parameter for multiphase fluid displacement simulations and impacts the overall productivity of the hydrocarbon reservoir [7].
Reservoir rocks, unlike man-made materials, tend to have spatially varying wettability, i.e., the contact angle varies from location to location on the surface of the grains in contact with the fluids in the pore space of a rock [8][9][10].The contact angle distribution is a result of the mineral composition and the geological history of a hydrocarbonbearing rock as well as surface texture, chemical composition of fluids (e.g., water, oil) in contact with the different crystalline planes of the grain minerals present in the rock, etc. Prior to contact with crude oil, rocks are almost always naturally water-wet.As oil invades the pore-space over geological time, the initial water-wet character may be altered in certain locations due to adhesion of surface-active materials from the oil phase.Mechanisms of wettability alteration depend on various local system properties such as pressure, temperature, mineral type, and fluids composition.
When rock samples are extracted from oil/gas reservoirs during the drilling process, they are often contaminated with mud and other fluids and chemicals.When those rock samples are subject to laboratory experiments to obtain two-phase relative permeability curves or capillary pressure curves, a common lab preparation procedure is to first clean the rock of fluids using solvents and then attempt to restore the rock's natural wettability.This process of attempting to reestablish fluid phase distribution and wettability distribution (e.g.contact angles) representative of the subsurface reservoir conditions is referred to as "aging".A typical aging procedure would involve cleaning the rock sample of original fluids and contaminants and treating the rock sample with chemicals that induce a water-wet condition.Next, the sample is saturated with brine, and in a "drainage" procedure oil is pushed in (hence brine is pushed out) using estimated reservoir conditions of temperature and pressure.The rock sample is then "aged" for a period of time, e.g. 4 weeks, allowing wettability alteration to occur, presumably in a fashion similar enough to what had occurred in the real subsurface rock formation so that a realistic wetting condition is restored.
With the advancement of imaging and computational fluid dynamics, it is now possible to perform a relative permeability or capillary pressure analysis numerically for rock samples.But the accuracy of these simulations also depends on the spatial distribution of the wetting condition of the grains among many other inputs [5].Using x-ray microtomography (micro-CT), it is possible to measure in-situ contact angles [11][12][13][14][15][16] in an aged rock sample from laboratory procedure and use as input for numerical relative permeability analysis.Observations of the relative fluid coverage of rock mineral surfaces can also be used to create a 3D map of the wetting state of the rock [8][9].However, establishing an initial fluid distribution and wetting conditions using numerical simulation remain challenging.It includes understanding the role of thin films in porous media.The type and thickness of films coating pore walls determines a reservoir rock wettability and whether or not a reservoir rock can be altered from its initial state of wettability [5,17].Thin film physics also depends on the grain surface mineralogy and the fluid pair.It is interrelated with the pore shape, especially pore curvature and surface roughness.As a result, the solution scope to this problem spans from molecular to microscopic scales, where molecular interactions determine the Derjaguin-Landau-Verwey-Overbeek (DLVO) forces [18,19], and interfacial tension, while microscopic pore topology determines candidate surfaces for wettability alteration.
In this work, we describe a numerical Digital Rock workflow starting from pore scale 3D image acquisition to relative permeability fluid flow simulation, including a novel wettability alteration process.Wettability alteration is carried out in a digitally reconstructed 3D pore-grain model, paying attention to the micro scale interactions.The resulting spatially varying wettability distribution is used for relative permeability simulation.The next section describes the methodology, including image processing, an overview of the numerical modeling approach, and detailed description of the wettability alteration process.Afterwards, wettability alteration and relative permeability results for a Berea sandstone can be found in the results section, which is followed by conclusions to shed some light on future research in this regard.

Sample imaging
This study uses a Digital Rock workflow [15] that starts with the acquisition of a 3D image of a Berea Sandstone sample using an x-ray micro-tomography (micro-CT).Fig 1(a) shows a 2D slice of the 3D scans captured perpendicular to the fluid displacement direction.Emphasis was given to acquiring images that are free of x-ray imaging artefacts such as beamhardening, shadowing, ring artefacts etc. Next, a segmentation workflow is applied to that image to identify the pores and grains as shown in black and white respectively in Fig 1(b).This segmentation uses differences in grey values and intensity gradients to identify voxels along the solid-pore boundary and produce the resolved pore structure to be used in the flow simulation.Further reference for image acquisition and processing can be found in [20].The stack of micro-CT images of a Berea sandstone was captured at 2.0 per voxel resolution to ensure good mesh quality for subsequent fluid flow simulation.At this resolution, the resulting 3D image shows sufficient pore connectivity in all three directions.The pore bottleneck radii, which is defined by the radius of the largest rigid sphere that can be percolated across the 3d rock pore space without getting stuck, were calculated to be 3.0, 3.6 and 3.2 voxels in the x, y and z directions respectively.Note that the bottleneck radius needs to be large enough to ensure sufficient accuracy of the fluid flow simulation and that requirement is solver specific.In a previously published work [20], it has been demonstrated that the Lattice Boltzmann method (LBM) solver used for this study can provide accurate results at this resolution.
Attention was also paid to ensure that the field of view of the 3D image is large enough to capture a representative elementary volume.A size growing analysis was performed, where the entire analysis domain is divided into incrementally bigger, overlapping, and centered subdomains of equal aspect ratio.For each of the sub-domains, porosity and single-phase permeability values are calculated as shown in Fig. 2(a).Additionally, a coefficient of variation of permeability as a function of sub-domain size is calculated in Fig. 2(b) [21].This chart represents likelihood that a certain subdomain permeability value will be one standard deviation of the mean value.Based on these results, we chose a subdomain of 800 voxels (corresponding to 1.6 mm) cubic size for subsequent multiphase simulation.

Numerical modelling approach
The multiphase fluid flow numerical simulation is conducted in the segmented 3D pore-space voxels of the selected rock domain using a multiphase LBM solver.The LBM is based on mesoscopic kinetic theory and solves a discrete form of the Boltzmann transport equations instead of Navier-Stokes equations [22][23][24][25][26].In this method, most of the operations can be performed locally on a cubic lattice, resulting in highly parallel computational performance [27].This simulation performance quality makes LBM an ideal candidate for porous media applications where typically a large domain with complex pore network needs to be simulated.The LBM has been successfully used for various porous media applications such as in hydrocarbon bearing rocks [20,28], gas diffusion layers in fuel cells [29,30], porous electrodes in lithium-ion batteries [31,32] etc.
The LBM implementation used in the current study is based on the multiphase framework presented by Shan and Chen [33] that later went through various improvements for simulating porous media applications with sufficiently high accuracy at low resolution and small viscosity [27,34].This is a diffuse interface method where different particle species are used to generate separated fluid phases, and the interaction between the different particle species determines the inter-component interfacial tension.Thus, the interfaces between different components are automatically determined once the species' interactions are defined.As a result, the complexity of explicit interface tracking such as in traditional volume-of-fluid or level-set approaches can be avoided.
In order to capture correctly the physics of multicomponent fluid flow, even at low resolution, the current LBM implementation uses, in addition to the voxel elements, triangulated surface elements also known as surfels, to accurately capture the pore-grain boundary as demonstrated in [35][36][37].The use of surfels to define the exact wall location results in more accurate fluxes near boundaries, and increased accuracy at lower numerical resolution compared with operating directly on the voxel grid of the micro-CT scan.Usage of surfels also improves the accuracy of local curvature and curvature-based pressure calculation as described in section 2.3.3.

Capillary pressure displacement
A capillary pressure driven primary drainage (oil displacing water) simulation is performed in the selected rock domain to establish a fluid phase (oil/water) distribution for subsequent water flooding simulation.Initially, the simulated domain pore space is fully saturated with water and a uniformly distributed strong water-wet wetting condition is applied to all pore-grain boundary surfaces.In all previous Digital Rock multiphase flow studies, the initial wettability condition is maintained constant along the simulation.In the current study, the wetting condition is altered during this drainage simulation using a process described in section 2.3.5 and compared with a traditional no-alteration scenario.
As shown in Fig. 3, the rock sample is placed inside a container box.The container has solid walls in the x and y direction, while in z direction (flow direction), oil and water buffers are placed in the opposite sides of the rock domain to control the capillary pressure.In addition, a numerical membrane acting as a porous plate is placed in between the rock domain and water buffer.This numerical membrane behaves as a volume-free and massless barrier that prevents oil from entering into the water buffer.Water on the other hand, can flow through the membrane without any resistance.Thus, the usage of this numerical membrane serves as a capillary pressure barrier for the non-wetting phase, similar to the porous plate in capillary pressure lab analysis [38].It is worth mentioning that this numerical membrane is a unique feature of the current LBM implementation we use, and that can practically replace the modeling of a finely resolved porous plate which can be computationally expensive and can potentially raise numerical issues.
During primary drainage, the pressure at the oil buffer is gradually ramped up so that oil starts to invade the pore space filled with water.At each pressure level the simulation continues until the saturation change is within a preset convergence value and a static equilibrium is achieved.

Relative permeability displacement
Using the initial fluid distribution, obtained from the previous capillary pressure drainage simulation, a multiphase unsteady-state kr relative permeability fluid flow simulation is carried out.In this unsteady-state kr method, a periodic boundary conditions and a driving force are applied in the flow direction, while no-flow boundary conditions are used in the directions perpendicular to flow.Further details of this procedure can be found in [20].

Curvature calculation
The pore-grain boundary, which consists of triangular surface elements, is used for both mean curvature calculation and to define wall boundary conditions for fluid flow.As described in [39], from a theoretical point of view, triangular meshes do not have any curvature at all, since all faces are flat and the curvature is not properly defined along edges and at vertices because the surface is not differentiable.One possible alternative is to think of a triangular mesh as a piecewise linear approximation of an unknown smooth surface and calculate mean curvature at the vertex.In the current study, we use an open-source, freely available tool "Visualization Toolkit" (VTK) [40] to calculate the mean curvature at every vertex: In the current workflow, since we need to calculate the curvature for every voxel in the computational domain, the curvature at the triangular surface element is calculated by: where, $ and % represent the triangular surface element and one of the vertex in that element respectively.Since the surface elements do not cross over from one cubic lattice voxel to another, curvature for every voxel can be calculated by area weighted averaging: where, + 123 | | represent area and curvature of a triangular surface element while 2 represents number of elements in a given voxel.

Force balance
The basic mechanism that drives wettability alteration is related to achieving a local force balance in the fluid-fluidsolid interfaces.It involves three forces: 1) Capillary pressure 4 5 , i.e. the pressure difference between non-wetting (oil) and wetting (water) fluids.Since the capillary pressure measurement and subsequent wettability alteration takes place when both of the fluid phases achieve static equilibrium, the difference between inlet oil pressure and outlet water pressure can also be considered as the capillary pressure of the system.2) Curvature based pressure 4 is the pressure caused by solid geometry forcing the shape of the oil/water interface into the shape of the pore geometry: where, 9 is the interfacial tension coefficient.3) Disjoining pressure Π is the pressure that prevents a thin water film from rupturing.As explained in [17], when oil invades a water-wet pore space, a thin water film is developed which coats and adheres to the solid surface.Disjoining pressure is a function of film thickness, h as shown in the schematic in Fig. 5. Positive disjoining pressure can be seen as the repulsive stabilizing force that keeps the oil away from the water-wet interface.Depending on the film thickness, either one of the following two repulsive forces can contribute to positive disjoining pressure: The first one is the electrostatic force that originates from the electrical double layers present at the solid/water and water/oil interfaces.The second one is a strong hydration force that plays dominant role when film thickness approaches molecular dimensions.On the other hand, dispersive van der Walls forces, are usually attractive and destabilizes thin water films contributing to the negative disjoining pressure.
The shape of the disjoining pressure isotherm, Fig. 5, suggests the existence of a local maximum disjoining pressure.When the applied capillary pressure exceeds the sum of the local curvature-based pressure and maximum disjoining pressure, the local thin water film can be considered ruptured making the adjacent grain surface a candidate for wettability alteration: It should be emphasized that Eq. ( 6) i.e. the augmented Young-Laplace equation, considers both disjoining pressure and local pore shape for wettability alteration [5].To illustrate the local pore shape dependency, consider a non-wetting fluid flowing through two separate micro channels, one with circular cross-section and other with star shaped crosssection.If they have the same bottleneck radius, as per Eq. ( 6), the thin film in a star shaped pore is more likely to be ruptured, even though the disjoining pressure and applied capillary pressure in star shaped pore is the same as the circular pore.Also, disjoining pressure is a necessary ingedient to explain capillary equilibrium for saturation high enough where only a smaller fraction of interfaces are really "free" interfaces but the majority of the oil-water interface follows the shape of the solid.In some instances the interfaces can even be flat, i.e. zero curvature.In order to sustain a force balance (Eq.6), the phase pressure difference has to be taken up by the disjoining pressure.1) The simulation is continued until static equilibrium is achieved for a certain capillary pressure 4 5 .
2) The surface elements in contact with the nonwetting phase are identified as shown by the red lines in Fig. 6(a) 3) The curvature-based pressure, 4 is also calculated using Eq. ( 5) at every voxel containing the poregrain surface boundary as described in section 2.3.3 and 2.3.4.4) Then Eq. ( 6) is used to determine whether a local water-film breakage may occur as represented by the red lines in Fig 6(b).5) Afterwards, a subset of surface elements from step 2) and 4) is identified as the candidate location for wettability alteration as shown by the red lines in Fig 6 (c).6) The wettability alteration can be a slow, dynamic process that occurs over many small discrete time intervals.In order to facilitate further change in fluid saturation, the simulation is continued with the updated wetting condition maintaining same capillary pressure until static equilibrium is achieved.7) Afterwards the pressure is ramped up to the next value and the above steps are repeated again.Thus, it is ensured that the fluid saturation is converged not only for certain applied pressure but also due to the corresponding change in wetting condition.It should be noted that the above procedure takes disjoining pressure for different minerals and fluid pair as an input.In the current study, we assumed a positive, repulsive disjoining pressure of 5000 Pa.As a next step, the workflow can be extended to calculate disjoining pressure using Molecular Dynamics simulation or from the experiments.

Wettability alteration in a Berea sandstone
Fig. 7 shows capillary pressure drainage simulation results for an oil/water system in a strongly water-wet Berea sandstone initially filled with water.The black and red measurement points in Fig. 7, 3 rd quadrants represent equilibrium states after pressure and wettability adaption respectively.At the beginning of the simulation, the water saturation changes are small in response to the increase in oil pressure.But as the pressure is ramped up enough to overcome the capillary pressure imposed by the local pore throat radii, oil rapidly invades the pore space as can be seen in between 90% and 70% of the water saturation, which is a phenomena called "Haines Jump" [41].As reported in [42], the Reynolds number calculated during Haines Jump are estimated to be in the order of one, indicating laminar flow, but the inertial forces are still relevant.The current simulation uses a Lattice Boltzmann formulation that represents the full transient Navier-Stokes equation which is essential to capture this phenomenon [18].However, not all the pore space is filled during Haines jump, especially during a slow convergence to static equilibrium.Even though Haines Jump process is far from quasi-static, during the subsequent slow convergence stage the phase pressure difference is not significantly large enough for the non-wetting phase to rapidly invade and the flow of oil can be described as quasi static and capillary dominant.In such a scenario, oil droplets can be snapped off when passing through a narrow pore-throat, similar to "Roof snap off" phenomena as described in [43,44] which can also be seen in current simulation as illustrated in Fig 8 .As can be seen in Fig. 7(c), oil breakthrough occurs for this rock sample at around 60% of water saturation, when the numerical porous plate stops further advancement of oil into the water buffer.Thus, it allows the oil pressure to build up, so that the oil non-wetting phase can intrude smaller pores too.As a result, having porous plate also helps achieving higher pressure-saturation points as shown in Fig. 7(d).For this simulation we carry out the wettability alteration after every four pressure steps by following the methodology described in section 2.3.5.In Fig 7, the 4 th quadrants show the water wet (blue) and oil wet (red) surfaces at different stages of the simulation.Only a small fraction of the surface area in contact with the oil non-wetting phase is subjected to alteration at low capillary pressures, because the applied capillary pressure is not large enough to overcome the repulsive disjoining pressure for thin film rupture.In contrast, alteration at higher capillary pressures creates larger fraction of oil-wet surface.By the end of simulation, almost 70% of the pore-grain surface is altered from strongly water-wet to oil-wet condition.Fig. 9 shows a comparison between two capillary pressure curves with and without wettability alteration.Drainage with wettability alteration exhibits smaller capillary entry pressure.Fig. 9 also confirms that the inclusion of the wettability alteration mechanism in the multiphase flow simulation leads to an overall lower water saturation in the drainage displacement.It can be hypothesized that for a given topology, the degree of observed difference between with and without alteration depends on the rock mineralogy.Different mineral and fluid pair can results in much higher or lower or even negative disjoining pressure which can potentially make the wettability alteration effect more pronounced.This phenomenon can be explained by the slow dynamic impact of wettability alteration that allows non-wetting fluid to intrude into very small pores, as the wettability can be changed in such a way that capillary pressure does not need to be overcome.In other words, the iterative process between wettability alteration, the consequential change of capillary pressure and establishing a new capillary equilibrium from the updated capillary pressure allows the contact line to move into smaller pores than without wettability alteration.Fig. 10 illustrates an example of a local pore scale displacement where the flow of oil is blocked by a narrow water-wet porethroat at certain applied capillary pressure.Next, the mechanism of wettability alteration of the pore-grain surface is carried out and certain portion of the surface in contact with oil is altered to more oil-wet condition, as shown in Fig 10(b).As the simulation is resumed with the updated wetting condition, but still at the same applied capillary pressure, the new oil wet surface helps oil to push through the small porethroat and displace more water as shown in Fig 10(c From this study, it is clear that wettability alteration has an impact on the overall displacment, but the alteration has been modeled as a discrete event happening every after certain pressure interval.We plan to investigate further the impact of discrete pressure interval of alteration as a next step.

Relative permeability analysis
The initial fluid and wettability distributions obtained from the capillary pressure drainage simulation are used as inputs for a subsequent relative permeability simulation where water displaces oil (imbibition).Fig. 11 shows the relative permeability simulation results in the Berea Sandstone model, for both altered and un-altered wettability conditions.Wettability alteration creates a mixed wet condition that allows oil to penetrate deeper into the pore space.As a result, the initial water saturation is lower than the un-altered wettability case (water-wet).For the mixed wet case, a thin oil film develops onto the grains where wettability alteration took place and the flow of water during water-flooding kr simulation may not be in direct contact with the grain wall.As a result, water relative permeability is higher for the altered wettability than the unaltered case where there is no oil film.The film also creates persistent oil connectivity even at higher water saturation which leads to more oil recovery for the alteration induced mixed wet case.

Conclusion
A numerical wettability alteration methodology based on a provided disjoining pressure isotherm has been established and implemented using a multiphase Lattice Boltzmann Method.In particular -1.An iterative process for wettability alteration has been developed based on force balance between curvature-based pressure, applied pressures and disjoining pressure.
2. This work also offers a secondary oil displacement mechanism based on slow dynamic wettability alteration process.This theory can be helpful to understand how oil can get inside into very fine pore space even though capillary pressure is not sufficient.
3. The effect of wettability alteration on overall oil production has also been illustrated by performing a subsequent relative permeability simulation.
In the future, this work can be extended to encompass molecular models of thin film physics, mineralogy and fluid pairs.Using Molecular Dynamics simulation, it is feasible to calculate interfacial tension and wettability on first principles basis as demonstrated in [45][46][47].Current Lattice Boltzmann workflow to calculate spatial wettability distribution can be integrated with a Molecular Dynamics workflow to obtain interfacial tension, contact angle and disjoining pressure isotherm inputs.

Fig. 1 .
Fig. 1.Two-dimensional section of 3D scans for the Berea Sandstone used in this work: (a) original grey scale, (b) segmented image.The stack of micro-CT images of a Berea sandstone was captured at 2.0 per voxel resolution to ensure good mesh quality for subsequent fluid flow simulation.At this resolution, the resulting 3D image shows sufficient pore connectivity in all three directions.The pore bottleneck radii, which is defined by the radius of the largest rigid sphere that can be percolated across the 3d rock pore space without getting stuck, were calculated to be 3.0, 3.6 and 3.2 voxels in the x, y and z directions respectively.Note that the bottleneck radius needs to be large enough to ensure sufficient accuracy of the fluid flow simulation and that requirement is solver specific.In a previously published work[20], it has been demonstrated that the Lattice Boltzmann method (LBM) solver used for this study can provide accurate results at this resolution.

Fig. 2 .
Representative elementary volume calculation: a) size growing analysis for porosity and permeability b) Permeability coefficient of variation, defined as the ratio of standard deviation to mean value.

Fig. 3 .
Fig. 3. Schematic of a capillary pressure drainage simulation setup mimicking a porous-plate method.

Fig. 4 .
Fig. 4. Curvature calculation: cubic lattice voxel with more than one triangular surface elements inside.| | | | 1 4 ‖ ⃗ ‖ | | 1 where, | |, ⃗ and represent absolute mean curvature, edge shared by two triangular surface elements and dihedral angle between the normals of the adjacent triangles at an edge ⃗ , respectively, as shown in Fig. 4 .This integral value is then normalized by the surrounding area of a vertex to calculate the curvature at the vertex.

Fig. 6 .
Surface selection (red) for wettability alteration: (a) selection is adjacent to a non-wetting phase, (b) selection satisfies 4 5 < 4 Π AB' condition and (c) selection is subset of (a) and (b).The wettability alteration workflow starts by creating a digital representation of the pore-grain surface and setting up the scenario for capillary pressure drainage simulation as described earlier.Afterwards the following steps are implemented to change the local wettability within a multiphase fluid flow simulation:

E3SFig. 7 .
Fig. 7. Various stages of wettability alteration simulation.The 1 st quadrant in each figure shows an iso surfaces of oil (beige) and a volume rendering of water (blue), 2 nd quadrant shows the water saturation profile along the flow direction, 3 rd quadrant shows capillary pressure versus water saturation chart and 4 th quadrant shows oil wet (red) and water wet (blue) surface within the pore space in 3D.

Fig. 10 .
Fig. 10.A pore scale event showing the slow dynamics of wettability alteration: (a) flow of oil is blocked by a narrow pore-throat, (b) local surface wettability alteration, and (c) updated fluid distribution in response to the wettability alteration.
). Wettability alteration shown in 10(b) is an outcome of procedure described in 2.3.5 that implements a steady-state solution of the augmented Young-Laplace equation.Once the wettability of the candidate surface is altered, the transient fluid flow simulation is resumed and the slow dynamic impact of alteration such as in Fig10(c) is simulated in time domain.