Density Functional Hydrodynamics in Multiscale Pore Systems: Chemical Potential Drive

. We use the method of density functional hydrodynamics (DFH) to model compositional multiphase flows in natural cores at the pore-scale. In previous publications the authors demonstrated that DFH covers many diverse pore-scale phenomena, starting from those inherent in RCA and SCAL measurements, and extending to much more complex EOR processes. We perform the pore-scale modelling of multiphase flow scenarios 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 numerical modelling of fluid transport in pore systems with voids and channels when the range of pore sizes exceed several orders of magnitude. Such situations are well known for carbonate reservoirs, where narrow pore channels of micrometer range can coexist and interconnect with vugs of millimeter or centimeter range. In such multiscale systems one cannot use the standard DFH approach for pore-scale modeling, primarily because the needed increase in scanning resolution that is required to resolve small pores adequately, leads to a field of view reduction that compromises the representation of large pores. In order to address this challenge, we suggest a novel approach, in which transport in small-size pores is described by an upscaled effective model, while the transport in large pores is still described by the DFH. The upscaled effective model is derived from the exact DFH equations using asymptotic expansion in respect to small-size characterization parameter. This effective model retains the properties of DFH like chemical and multiphase transport, thus making it applicable to the same range of phenomena as DFH itself. The model is based on the concept that the transport is driven by gradients of chemical potentials of the components present in the mixture. This is a significant generalization of the Darcy transport model since the proposed new model incorporates diffusion transport in addition to the usual pressure-driven transport. In the present work we provide several multiphase transport numerical examples including: a) upscaling to chemical potential drive (CPD) model, b) combined modeling of large pores by DFH and small pores by CPD.


Introduction
In recent years we witnessed an increasing level of capability and acceptance of digital rock as a complementary tool in core analysis. This tendency is following a gradual increase in both resolution of rock imaging techniques and the power of high-performance computing that are the key to unlock the potential in the core analysis by digital rock. As an answer to the growth in the hardware possibilities various modeling methods applicable in digital rock also evolve; see reviews in [1][2][3][4][5][6]. One of such methods, which the authors are actively developing, is the density functional hydrodynamics (DFH).
DFH is essentially a pore scale multi-phase, multicompositional approach. As it has been demonstrated in previous publications [7][8][9][10], various multiphase compositional problems can be described in the frame of the DFH. The DFH method has found many applications in direct pore-scale modeling of various hydrodynamic and petrophysical phenomena on digital rock models obtained by X-ray micro-CT and SEM [11][12][13][14][15][16][17][18][19][20][21]. The practical usefulness of a standard workflow is very much dependent on the possibility to resolve necessary pore structure adequately, as well as on the rock sample being representative for the considered formation. There are specific ways to meet both of these requirements, but in this work, we discuss the resolution problem only.
In many cases the natural rocks are characterized by a wide range of pore sizes making the choice of X-ray micro-CT resolution very difficult if not altogether impossible. Indeed, if pores of both micrometer and millimeter range coexist within one piece of rock then, at present, it is not always possible to develop 3D model of rock microstructure, which can represent all pores. In situations like this we propose a combined approach, when several 3D digital rock models having different resolutions are integrated into one synthetic model. Combined multiscale approach is also being developed in the frame of pore-network modelling [1], but we follow image-based approach as it is inherent in DFH. High-resolution sub-models with explicitly resolved small pores are used for direct modeling of transport processes by DFH to obtain effective flow characteristics. Then these results are incorporated into low-resolution model with explicit imaging of large pores while small pores are not necessarily resolved. At this stage the hydrodynamics in large pores is still simulated directly by DFH, while the transport in matrix is simulated using the previously computed effective transport properties. All these operations are done by the DFH pore-scale simulator called DHD (Schlumberger).
In effect, the described procedure represents upscaling from pore-scale DFH equations to macroscopic porous medium equations. As concerned with the single-phase flow, our idea is similar to the one described in [22,23]. However, in addition to being applicable in general multiphase compositional scenarios, the rigorous upscaling starting with the DFH equations leads to the transport equations, which constitute the significant extension of the Darcy model. The driving force for the fluid flow happens to be gradients of chemical potentials of fluid chemical components that gives the model its name Chemical Potential Drive (CPD). Because change in pressure is directly related to change in chemical potential, this new model is consistent with pressure-drive Darcy approach as will be shown below. On the other hand, the CPD approach covers certain physical phenomena outside the scope of the Darcy model. Here are two examples in relation to the latter assertion. First, there are cases, when definition of pressure in confined fluid is problematic, e.g., when gas mean free path length is comparable to pore size (high Knudsen number), or liquid in pores with significant disjoining pressure; at the same time, the definition of chemical potential is still correct in both said instances. Second, there is an easy to derive exact solution for two-phase water-oil transport velocities in a hydrophilic circular capillary where the water and oil velocities w o , u u are defined as phase fluxes per capillary cross-section, sw so ,   are the shear viscosities of water and oil, respectively, R is the capillary radius, r is the radius of the central region occupied by oil, and G is the pressure gradient. An important observation following from Eqs. (1), (2) is that the transport of oil is influenced by viscosity of water, because oil, being the central phase, is carried by water. This cross-phase friction is in direct contradiction with the routine concept of phase permeabilities indicating the possible role of cross-terms in phase permeability matrix [24]. At the same time the cross-component influence is inherent in CPD, as it will be demonstrated below.
Limitations of the Darcy law are discussed in the literature and deviations from it are observed in both high-accuracy 4D visualization of pore-scale fluid dynamics as well as in numerical simulations of porescale flows [25]. One may expect more deviations when dealing with complex fluids. In Sec. II we give a brief reminder of the DFH equations, as well as a derivation of upscaling from DFH to CPD. Numerical demonstration of combined DFH and CPD models is presented in Sec. III. The overall summary and discussion of results is in Sec. IV.
The summation over repeated indices is implied everywhere. The indices , , 1, 2, 3 a b c  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:

Density functional hydrodynamics
Here we provide only some of the basic definitions necessary to evolve the step from the DFH equations to the CPD ones. 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 having volume elsewhere in continuum mechanics, the small volume limit is understood as the convergent procedure with D V being small, but still larger than the molecular volume.
By counting the flow rate of molecules through a small area inside the mixture, one can define the component flux where by definition diffusion flux does not influence net mass transfer We assume the existence of the Helmholtz energy functional where D  is the boundary surface of the region D (when the region is finite), is the bulk Helmholtz energy density of homogeneous mixture, i j  is the positive-definite symmetric matrix, and is the surface Helmholtz energy density, which is not equal to zero if D  is a contact surface with some immobile solid. It is convenient to recollect certain thermodynamic equations involving Helmholtz energy where i  is the chemical potential of the i th component, and p is the hydrostatic pressure. The model in Eqs. (4) and (5) 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 [9][10][11][12][13][14][15][16][17][18][19][20][21]26]. Now we list the DFH statements necessary to move forward to CPD.
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 at D  , when it is a contact surface with some immobile solid It is possible to derive the differential equation for the total energy (that includes both Helmholtz and kinetic energy) with energy flux a J and dissipation rate function  : In accordance with Eqs. (11), (12), and (22) the inequality 0   holds. Also, in accordance with Eqs.
in consistency with the second law of thermodynamics.

Upscaling to macroscopic transport in porous medium: chemical potential drive
Now we consider transition from microscopic DFH to macroscopic CPD description of fluid transport in porous medium. We use variables related to microscopic pore-scale description in parallel to similar variables related to macroscopic description. If a is a microscopic variable, then similar macroscopic variable is denoted by symbol a . The latter variable is usually obtained from the former by using certain averaging procedure. For example, we define component molar density by is now sufficiently large to contain representative piece of saturated porous medium. The respective spatial region D encompasses pores belonging to a porous medium region  containing both solid skeleton and pores. This porous medium region  is supposed to be representative of pore structure in statistical terms but not necessarily identical to the neighbor regions of porous medium. This is different from homogenization approach with periodic geometry of pore system [27]. The volumes of two regions D and  are related by is porosity. The processes in a porous medium are assumed to be sufficiently slow in respect to the changing macroscopic component molar densities, so that the local distribution of chemical components and phases in pores are close to the equilibrium. Therefore, it is possible to define the macroscopic Helmholtz energy density with mixed no-slip and free boundary conditions. In terms of differential operations Eq. (20) can be represented as follows where the left side is positive elliptic self-adjoint operator in 2 ( ) L D acting on the velocity field. It admits the solution with a Green's function, which can also be represented in operator form Here the integral operator in Averaging Eq. (23) over different cross-sections of the region  produces the macroscopic flux of components generated by gradients of chemical potentials (Chemical Potential Drive or CPD) where the transport matrix Direct calculations lead to the macroscopic equivalent of

Boundary conditions between DFH and CPD domains
In highly heterogeneous porous material there can be situations when there are adjacent spatial regions where one can use alternatively DFH or CPD. It is necessary to discuss what boundary conditions should be set at the common boundary where variable X is determined by the additional boundary equation

Examples of CPD models
For the first example we assume that the porous medium is homogeneous and isotropic, and there is single-phase single-component flow. In this case the CPD transport law is equivalent to the Darcy law (see the last equation in (6)

Numerical examples
To demonstrate application of a workflow based on the joint DFH modeling within resolved porosity and CPD modeling within matrix containing unresolved smaller pores (DFH+CPD workflow) we chose a heterogeneous, but relatively well characterized chalk sample. The sample is grainstone composed of pelloids, skeletal grains and ooids. The grains are cemented by equant calcite spar cement. Intergranular porosity is the dominant pore system, additionally some moldic and micropores are present in partially leached grains. Pore throats are ranging between approximately 1 and 10 microns. Pore throat distribution is skewed to smaller pore throat sizes. Measured gas porosity of the sample is 26.2%, and gas permeability is 19.4 mD. An 8 mm diameter mini plug was scanned with 2.46 um per voxel resolution, and a portion of the mini plug was scanned with higher resolution of 0.82 m per voxel (Fig. 1).
Two 3D digital rock models (DRM) with these two different resolutions were constructed (Fig. 2). Both models cover the same spatial region with the size 0.82 mm x 0.82 mm x 0.82 mm, but they have different number of cells (voxels). The high-resolution DRM has 1000 3 (1000x1000x1000) cells and the low-resolution DRM has 333 3 (333x333x333) cells so that the cell sizes are 0.82 um and 2.46 um, respectively. In either model some portion of pores fell below resolution. However, it is evident from Fig. 2 that significantly higher proportion of pores remains below resolution in the low-resolution DRM than it does in the high-resolution DRM. The data on the resolved and unresolved porosity is shown in Table 1.  In addition to X-ray micro-CT, we have done scanning electron microscopy (SEM) of parts of the sample to reveal the pores below micro-CT resolution down to 50 nm. Using both micro-CT high-resolution (0.82 um) and SEM data we performed digital modeling of MICP experiment and reconstructed pore throat size distribution that is compared to the standard experimental MICP measurements in Fig. 3. The two curves in Fig. 3 match well meaning that most of the controlling pore throats have been adequately resolved on the corresponding scales.
In modeling on the high-resolution DRM, we assumed that the subresolution matrix portion of the model is effectively impermeable for oil at relevant oilwater capillary pressure levels so that only the DFH equations need to be solved within the resolved pore space. On the other hand, in modeling on the lowresolution DRM, we assumed that the subresolution matrix is permeable and therefore the full DFH+CPD simulation needs to be done. To set up a single-phase DFH+CPD simulation it is necessary to provide absolute permeability in the subresolution matrix cells so that the single-phase CPD transport law given in Eq. (33) could be resolved. In order to do this, we first simulated singlephase flow on the high-resolution DRM to obtain velocity distribution (Fig. 4a). Then, by interpreting the simulation results in terms of Eq. (33) and linking them with the gray scale distribution from the high-resolution images such as those shown in Fig. 1c, we established a correlation between the gray scale and absolute permeability values. We then applied the same correlation to the low-resolution gray scale images to inform the low-resolution DRM about the subresolution matrix permeability. The velocity distribution obtained on this 333 3 low-resolution DRM is presented in Fig. 4b. The absolute permeabilities calculated on the high-and low-resolution DRMs are 17.3 mD and 18.2 mD, correspondingly.
The next step was simulating two-phase flow and the steady-state relative permeability experiment. The two immiscible phases were assumed to be water and oil with the properties as follows:  is the oil-water interfacial tension. Since the sample is carbonate, it was assumed that it has mixed wettability. To model this, we first injected oil into initially water-wet 100% water saturated high-resolution DRM using our standard DFH workflow [11,[14][15][16][17][19][20][21]. (In all cases here and below injection was arranged in direction from bottom to top as shown in the 3D views.) Then at those places where oil touched the pore walls we switched wettability to moderately oil-wet while leaving the untouched pore walls water-wet. In this way most of the smaller pores as well as the pores isolated by thin pore throats remained water-wet, while the larger pores attained wettability towards oil. The distribution of water and oil at the residual water saturation state is shown in Fig. 5a. Using the obtained mixed-wet wettability distribution we simulated the steady-state relative permeability experiment by injecting oil and water mixture in several steps different in water content in the influx; at the final step, with only water injected, the residual oil saturation was obtained (Fig. 5c). To run DHD+CPD two-phase simulation, it is necessary to provide CPD transport matrix of Eq. (24) in subresolution matrix cells. Alternatively, when dealing with immiscible phases like oil and water that are each described by just one pseudo-component, it is possible to use simplified CPD transport model provided in Eq. (37) and based on the concepts of absolute and relative permeabilities and capillary pressure. We already discussed how absolute permeability was assigned to the subresolution matrix cells of the low-resolution DRM. The relative permeability was assigned in a similar way since all the velocity distribution, pressure, and water saturation for the subresolution matrix cells are known from the simulation on the high-resolution model. The capillary pressure function was constructed using the Young-Laplace equation and then scaled to the pore size distribution correlated by the gray scale values.
The distributions of oil and water obtained in simulation of the two-phase steady-state relative permeability experiment modeled using the DFH+CPD workflow on the low-resolution DRM are presented in Fig.5b,d. The chart in Fig. 6 compares the relative permeability curves obtained on both high-and lowresolution DRMs.
The final example we wish to present here is the simulation of both absolute and steady-state relative permeabilities on another low-resolution DRM with the size of 4.92 mm x 4.92 mm x 2.46 mm approximated with 500 x 500 x 250 cells with the cell size of 9.84 um. This model covers maximum size parallelepiped that can be extracted from the micro-CT scans of the 8 mm mini plug as shown by the dashed line in Fig. 1a. All of the parameters necessary to set up the DFH+CPD model were taken from the previous calculations on the highresolution DRM as described previously.  6. Comparison of the simulated steady-state relative permeability curves obtained on the high-resolution (orange and blue curves) and low-resolution (yellow and gray curves) DRMs of the same spatial region within the sample.
Simulated absolute permeability of this model is 19.8 mD. The distributions of water and oil corresponding to the residual water saturation and residual oil saturation as obtained during the steady-state relative permeability simulation are shown in Fig. 7. The simulated relative permeability curves are presented in Fig. 8.
In this work, the simulation results were interpreted in consistency with the conditions in Eq. (37) that is with the traditional Darcy transport. However, in complex fluid problems deviations from relative permeability concept are expected, which then must be treated in the frame of the wider CPD model.  Fig. 7. 3D view of the distribution of water and oil in the upscaled low-resolution DRM at (a) initial conditions corresponding to the residual water saturation, and (b) final state corresponding to the residual oil saturation. Water is shown in semitransparent blue, oil is shown in red, and skeleton is transparent. Fig. 8. Simulated steady-state relative permeability curves obtained on the upscaled low-resolution model.

Conclusion
We developed a new digital rock workflow applicable for simulation of multiphase flow on multiscale models representing heterogeneous core samples. The workflow called DFH+CPD is based on solving the standard DFH equations within resolved porosity and solving the new CPD equations within unresolved porosity or subresolution matrix. The CPD equations are derived in mathematically rigorous way from the DFH equations averaged over statistically representative ensemble of pores. The new workflow amounts to an effective upscaling procedure allowing multiphase simulation on coarse resolution digital rock models that include information that can be obtained on high-resolution models. The DFH+CPD workflow was implemented within the Schlumberger DHD simulator.
To validate the new workflow, we compared absolute and steady-state relative permeabilities simulation results on the high-resolution model and the low-resolution model of the same region within the heterogeneous carbonate core. Then we demonstrated simulation of the steady-state relative permeabilities on the upscaled digital rock model that covers almost entire 8 mm mini plug.
We thank Schlumberger for permission to publish this work.