The Problem of Stability of Gas-Condensate Mixture at Pore-Scale: The Study by Density Functional Hydrodynamics

. The method of the density functional hydrodynamics (DFH) is used to model compositional gas-condensate systems in natural cores at pore-scale. In previous publications, it has been demonstrated by the authors that DFH covers many diverse multiphase pore-scale phenomena, including fluid transport in RCA and SCAL measurements and complex EOR processes. The pore-scale modeling of multiphase flow scenarios is performed by means of the direct hydrodynamic (DHD) simulator, which is a numerical implementation of the DFH. In the present work, we consider the problem of pore-scale numerical modeling of three-phase system: residual water, hydrocarbon gas and hydrocarbon liquid with phase transitions between the two latter phases. Such situations happen in case of gas-condensate or volatile oil deposits, in oil deposits with gas caps or in EOR methods with gas injection. The corresponding field development modeling by the conventional reservoir simulators rely on phase permeabilities and capillary pressures, which are provided by laboratory core analysis experiments. But the problem with gas-liquid hydrocarbon mixtures is that in laboratory procedures it may be difficult or even impossible to achieve full thermodynamic equilibrium between phases as it must be under the reservoir conditions of the initial reservoir state. However, reaching the said equilibrium is quite possible in numerical simulation. In this work, the gas-liquid mixture, after being injected into core sample, would slowly undergo the rearrangement of the phases and chemical components in pores converging to the minimum of the Helmholtz energy functional. This process is adequately described by DFH with consequent impact on phase permeabilities and capillary pressure. We give pore-scale numerical examples of the described phenomena in a micro-CT porous rock model for a realistic gas-condensate mixture with quantitative characterization of phase transition kinetic effects.


Introduction
Core analysis methods based on using digital rock models (DRMs) gradually become more and more accepted as a complimentary tool to traditional lab core analysis.Although not free from its own difficulties, digital rock approaches can help to obtain data and reduce uncertainties in cases where traditional lab-based core analysis becomes especially challenging to perform, e.g., poorly consolidated cores, low permeable cores, complex fluids like live oils and gas-condensates.Currently there are known many alternative methods used for numerical simulation of dynamic processes on DRMs and we refer the reader to reviews in [1][2][3][4][5].The modeling method developed by the authors is called density functional hydrodynamics (DFH).
In previous publications it has been demonstrated that various multiphase compositional problems can be modeled in the frame of the density functional hydrodynamics (DFH) [7][8][9][10][11][12][13][14][15][16][17], including practical digital rock applications [18][19][20][21][22][23][24][25][26].In the present work we discuss the pore-scale modeling of gas-condensate mixtures in digital rock models obtained by X-ray microCT.The following are specific features of this class of problems: a) the thermodynamic properties of hydrocarbon mixtures are related to the multicomponent nature of such mixtures, so it is important to retain in the modeling the influence of different components; b) gas and liquid phases undergo chemical component exchange, which can be slow in comparison to the transport processes observed at the laboratory scale.The latter problem can affect the applicability of laboratory results to the gascondensate reservoir modeling.Namely in the reservoir far from the production and injection wells the transport processes in pores are sufficiently slow, so the gascondensate mixture at any moment can be close to the minimum of the Helmholtz energy functional, while under laboratory protocols the state of the mixture can be in some metastable state far from the true equilibrium state.In order to demonstrate this problem, we perform the simulation of transport of gas-condensate mixtures in two different scenarios: a) in accordance with laboratory protocols; b) using equilibrium gas-condensate states as initial input for the computations of phase permeabilities and capillary pressures.The results highlight the possible difference between two approaches, which in our opinion proves the unique value of digital rock analysis in this class of problems.
The summation over repeated indices is implied everywhere.The indices , , 1, 2,3 abc = are related to Cartesian coordinates a x , the indices , , 1,..., i j k M = are related to chemical components in the fluid mixture.We consider isothermal processes, so the temperature is assumed to be fixed, and dependence of certain variables on temperature is omitted.We use short symbols for partial derivatives:

Reminder of density functional hydrodynamics
Here we provide the reminder of some of the basic concepts of DFH, which are necessary for the modeling of gas-condensate mixtures.A detailed exposition of the DFH can be found in [7][8][9][10].
We consider continuum mechanics description of a mixture of M chemical components present inside a spatial region D with overall volume . Like elsewhere in continuum mechanics, the small volume limit is understood as the convergent procedure with D V being small, but still large in comparison to the molecular volume.The molar densities ) where D  is the boundary surface of the region D (when the region is finite), The model in Eqs. ( 1) and ( 2) is adequate for description of many important phenomena involving multiphase multicomponent mixtures.Up to now, it was successfully used to simulate multiphase multicomponent phenomena with or without phase transitions, surfactants, and mixtures with solid phases such as gas hydrates or solid particles [3][4][5][6][7][8][9][10][11].
The multiphase compositional transport in case of isothermal flow of fluids with Newtonian rheology is governed by the following equations written in a form of conservation laws for chemical components and momentum together with the relations ) and subject to the boundary conditions Since there exists the identity () If both phases A and B are in contact with solid surface and ** ( , the phase with lower surface Helmholtz energy density will form a layer in contact with solid surface.

Modeling of gas-condensate mixtures by density functional hydrodynamics
Our aim is to apply DFH for the numerical modeling of gas-condensate mixtures inside of the pore system of natural rock sample using 3D microCT X-ray models of the rock microstructure.Previously we have demonstrated how the DFH method can be used to model condensate banking [15] though in less rigorous mixture composition formulation than the one proposed in this work.
Here we use a Cartesian grid with cubic cells in one-to-one correspondence with 3D voxels of segmented microCT images.The numerical realization of DFH equations is described in previous publications [10,13].Here we discuss specific numerical features associated with gas-condensate modeling problem.
First, it is necessary to specify the bulk Helmholtz energy density , which represents bulk mixture thermodynamics including the phase behavior.Currently it is widely accepted, that the thermodynamic behavior of hydrocarbon mixtures is adequately described by Peng-Robinson equation of state (PR EOS) [27], so it is natural to rely on the corresponding analytic expression of Helmholtz energy.The main computational problem with gas-condensate mixtures is that the large number M of chemical components leads to prohibitively large computational time.The way to remedy this problem is to perform computations in terms of smaller number of pseudocomponents, which are artificially constructed from the initial set of chemical components.In our case the transition from the true components to pseudocomponents (lumping) must be compatible with the DFH equations ( 4)-( 13).Now we shall describe the mathematical procedure, which helps to reduce the number of effective components.
Let us consider linear transformation of molar densities with some invertible matrix i   (here and below index 1,..., M  = , and the summation is assumed over the repeated indices) as follows We define an auxiliary vector g  with all components equal to 1. Since the matrix i   is invertible there exists the inverse matrix . We assume that the matrix i   satisfies the following additional condition This condition means that the new overall molar density '' coincides with the old one n .Therefore, there exists the linear relation between the new and the old concentrations, as follows with the normalizing condition for concentrations .So far, the transformation ( 16) is purely formal and has no physical sense, but it can help to reduce the number of effective components.Indeed, let us consider the phase transitions in the narrow pressure range, when the mixture with overall composition ' ' In this case the hydrodynamic processes for the considered mixture, i.e., phase transitions, convective transport, and mixing, negligibly change the quantities (19).Then one can assume that there is a fixed quantity The problem is to find the matrix i   , which provides as many approximations like (19) as possible.Here is the mathematical procedure for the construction of this matrix.
In digital rock transport problems, the variations of pressure are small because of the small size of microCT rock models, so the matrix i   is constructed for certain fixed pressure.At the same time variations of local composition can be significant, so we consider variations in the range with fixed positive parameter h .Note that the concentration vector i c belongs to the hyperplane small.Thus, we described a new lumping technique that is consistent with the DFH.The technique is purely mathematical and can be classified as a space reduction approach.
In this work, we chose the definition of i   with 1 h = .After defining the new pseudocomponents, it is necessary to calculate the PR EOS Helmholtz energy function in terms of these new variables.To achieve further reduction in computational time, we approximated this function with Padé-like rational approximation, which is used in numerical solution of the DFH equations.Besides being a thermodynamic description of bulk fluid, this final version of the Helmholtz energy function is used to fit the matrix ij  to experimental IFT data.
It is necessary to express the surface Helmholtz energy density * f and transport coefficients ij D through the new variables.We use a linear approximation for * f , which is compatible with Young's equation as it has been discussed in Section 2.1 (see also [16]).It is necessary to note that in digital rock problems the surface Helmholtz energy function * f must also compensate the difference between actual specific surface and smaller specific surface in microCT rock models.The transport matrix ij D is initially specified using known data on diffusion [28] and then it is transformed in accordance with Eq. ( 16) using the relation '

Numerical results
In interpretation of the numerical results for final static states it is necessary to remember that in accordance with equilibrium conditions ( 14) the chemical potentials for gas and liquid are equal and constant everywhere except for interfacial regions and regions in vicinity of the rock surface: Among other things this means that the composition of gas phase is the same everywhere where the conditions ( 23) take place, and the same statement is true for liquid phase.But variations of composition in grid cells, which cover the interfacial regions or are adjacent to the rock surface are permissible.
In this work, we studied a rich gas-condensate system corresponding to a gas-condensate in southern Russia.The reservoir pressure is 50.6 MPa and temperature is 82 C. The dew point is at 48.7 MPa.A somewhat simplified composition of the reservoir gas, i.e., above dew point, is listed in Table 1.21) are presented at Fig. 1.One can see that if maximum relative error at 1 h = (that is at 100% deviation from the original mixture composition) is set to be 5%, then the reasonable choice is 4 m = , because for 4 m  the error grows rather drastically and goes well above 5%.Thus, we settle for a 4-pseudocomponent space that will be used in numerical simulation of gascondensate flow on a digital rock model.
A core sample used in this work is a lowpermeable, but relatively well sorted carbonate rock.The sample was X-ray scanned at 0.82 um/voxel spatial resolution.Image segmentation was based on an extended version of the Indicator Kriging approach [29].An example of a 2D cross-section of the sample as well as a 3D view of the DRM used for numerical simulation are shown in Fig. 2. The porosity of the DRM is 0.175 and permeability is 1.8 mD.The DRM resolved porosity is somewhat lower than the 8 mm diameter core plug measured porosity, which is 0.225.The core plug measured permeability is 2.5 mD.This is a rather typical situation, especially in low-permeable rocks, when a considerable portion of porosity appears below the X-ray microCT resolution.For such cases, we have developed a multiscale workflow based on resolving a core plug at a range of different resolutions and using high-resolution models for determining effective properties that are then used in a combined multiscale low-resolution model [17].That workflow is fully compatible with the compositional gas-condensate modeling in the frame of the space reduction approach, but we do not use it in the present work to avoid unnecessary complication and blurring the focus of this paper.
As the first step in numerical simulation on the DRM, we have modeled initial saturation of the reservoir gas and connate water (Fig. 3) by solving a full thermodynamic equilibrium problem (global minimum of Helmholtz free energy) following the methodology described in [25].Spatial distribution of wettability was established by surface Helmholtz energy density field ( * f ) with correlated stochastic noise representing expected natural variability in rock surface properties.The connate water saturation is Sw=0.18.According to the mixed-wet nature of the core connate water is primarily accumulated within small pores.As the next step of numerical simulation, we have evaluated gas-condensate relative permeabilities at connate water.In order to obtain different condensate saturations, the simulations were conducted at different pressures below dew point.Also, two cases distinct in the simulation workflow have been modeled.
In Case 1, gas-condensate mixture corresponding to selected pressure below dew point was flowing through the DRM until stabilization of phases saturations within the model and then gas-condensate phase permeabilities were evaluated for each saturation point similarly to how we usually do this in case of immiscible phases [23].In Case 2, after stabilization of phases saturation, we have waited until thermodynamic equilibrium takes place and pseudocomponents redistribute to minimize global Helmholtz free energy.In each case, the injection was arranged at a constant rate corresponding to a low capillary number flow.The mixture was injected at the bottom side of the model and arriving fluids were removed at the top side to sustain constant overall pressure.The simulations were carried out using the following parameters (inferred for the considered deposit at 42 MPa): sc 0.92 mPa×s As can be seen from the pictures, the heavy pseudocomponents (shown in red and yellow colors in Figs. 6 and 7) tend to group closer to connate water and occupy smaller pores, while bigger pores are mostly occupied by the lighter pseudocomponents (shown in green and aqua colors in Figs. 6 and 7).This is the consequence of the fact that condensate, which is dominated by heavy pseudocomponents, has lower interfacial tension with water than gas that is mostly formed by light pseudocomponents.The said details are seen clearer from the 2D cross-sections in Figs. 6 and 7.During the flow gas-condensate mixture is not necessarily in equilibrium state everywhere (here it is worth stressing that the DFH equations are compositional and allow for non-equilibrium flows; thermodynamic equilibrium is produced as a result of global minimization of the Helmholtz free energy and the minimization process is coupled with the transport phenomena).Therefore, it is necessary to establish definition of phases for arbitrary mixture composition [10,12].This is done by the bulk Helmholtz energy analysis, which has minima corresponding to the equilibrium phases.If mixture composition appears within convex part of the Helmholtz energy around a phase then this composition corresponds to this phase.If, however, the composition appears outside any convex regions then such composition is unstable and does not belong to any of the phases.
A comparison between gas-condensate relative permeabilities is presented in Fig. 8.The gas and condensate saturations have been normalized relative to the total volume occupied by the hydrocarbon phases (connate water was excluded).It is noticeable that in Case 2 the two-phase mobility region is somewhat narrower than in Case 1.For example, condensate mobility threshold, Sc0, is equal to 0.28 in Case 2, while in Case 1 it is 0.23.This is explained by the fact that after fluid injection to saturation stabilization, the nonwetting phase tends to form a cluster connecting the injection side of the model with the opposite side.But after further stabilization towards thermodynamic equilibrium that cluster can deform, become less connecting, or can even break entirely.It is only because of low interfacial tension between gas and condensate that the difference in relative permeabilities between Cases 1 and 2 appeared to be rather modest; it might have been far greater have the phases had bigger IFT.

Conclusion
In this work, we developed a compositional mixture lumping technique that is fully compatible with the DFH.This new lumping technique is a purely mathematical procedure that can be classified as a space reduction approach.As a result of its application, a pseudocomponent molar density space is built that has fewer dimensions than the space of the original molar densities.However, the inverse transformation exists that allow to reconstruct the original components from the new pseudocomponents.The error of such transformation can be estimated beforehand and can be used to select the optimal number of pseudocomponents.
Starting from an 11-component gas-condensate system we have arrived to a 4-pseudocomponent reduced system yielding maximum error of 5% within a subspace specified by a 100% deviation from the original mixture composition.Using this reduced system, we have conducted digital rock simulations by the DFH and obtained gas-condensate relative permeabilities.The simulations have been carried out in two different statements -a classical one similar to a standard laboratory protocol, and a new one attempting to mimic thermodynamic equilibrium conditions that can form over prolonged periods of time.It was observed that the second statement provides somewhat poorer overall mobility of phases as compared to the classical formulation disregarding thermodynamic equilibrium.

I
flow rate of molecules through a small area inside the mixture, one can define the component can be represented as a combination of transport term to derive the expression for interfacial tension (IFT), which helps to tune the matrix ij

.
The transformation(16) can be applied to the DFH equations (4)-(13).The mass velocity a v and mass density '' ii m n m n   == are not transformed given we transform molar masses ' ii mm   = .The representation of the DFH problem (4)-(13) in new variables is evidently equivalent to the old representation, and any solution in terms of ', a nv  can be used to compute solution in terms of , a nv  with corresponding reduction in system of Eqs.(4).If there are similar situations with other concentrations ' c  , the number of effective component densities ' n  can be reduced even further.
).Using PR EOS for any vector i c we can calculate corresponding gas and condensate concentrations () .If by chance we have one-phase system, we put both compositions equal to i c .Let us define by induction the sequence of vectors the error function (21) is solved numerically.Using this solution, which is dependent on parameter h , it is possible to define the matrix i be interpreted as new effective pseudocomponent densities.The cumulative error of this approximation is characterized by the function m J .The analysis of this function can help to determine the best minimum number m , when the error is still sufficiently

Fig. 2 .
Fig. 2. (a) Micro-CT grayscale cross-section of the carbonate core sample at 0.82 m resolution, (b) 3D view of the segmented pore space of a cubic DRM with dimensions 500 3 at 0.82 m resolution.

Fig. 3 .
Fig. 3. 3D view of the distribution of the reservoir gas shown in green and connate water shown in semitransparent blue within the pore space of the DRM.
are the condensate-gas, condensate-water, and gas-water interfacial tensions, respectively.Some examples of phases saturation and pseudocomponent distributions obtained during simulations are presented in Figs. 4 through 7. Fig. 4 compares gas, condensate and connate water saturations in Cases 1 and 2 obtained after injection of gascondensate mixture at 42 MPa pressure.Although the distributions are generally similar, there are visible differences in details.Fig. 6 compares distributions of pseudocomponents that correspond to the saturations from Fig. 4. Similarly, Fig. 5 and 7 show distribution of phases and pseudocomponents, respectively, for lower value of pressure equal to 33 MPa.
m J of Eq. (