SCA2013-014 1/12

DIRECT HYDRODYNAMIC SIMULATION OF

MULTIPHASE FLOW IN POROUS ROCK
D. Koroteev1, O. Dinariev1, N. Evseev1, D. Klemin1, A.Nadeev1, S. Safonov1,
O._Gurpinar2, S. Berg3, C. van Kruijsdijk3, R. Armstrong3, M. T. Myers4, L. Hathon4,
H. de Jong4
1
Schlumberger Moscow Research, 13, Pudovkina Str., Moscow 119285 Russia
2
Shlumberger Reservoir Characterization Group, 1675 Broadway, Suite 900, Denver, CO
80202 USA
3
Shell Global Solutions International, Kesslerpark 1, 2288 GS Rijswijk, Netherlands
4
Shell Global Solutions International, 3333 Highway 6 South, Houston, Texas 77251-
7171, USA
This paper was prepared for presentation at the International Symposium of the Society of Core
Analysts held in Napa Valley, California, USA, 16-19 September, 2013

ABSTRACT
We present various numerical studies conducted with a novel pore-scale simulation
technology called Direct Hydrodynamic (DHD) Simulation that can be used to study
multiphase flow at various scales ranging from individual pore-scale events to complex
scenarios like capillary de-saturation and relative permeability of digitized rock samples.
DHD uses a diffuse interface description for fluid-fluid interfaces that is implemented via
the density functional approach applied to the hydrodynamics of complex systems. In
addition to mass and momentum balance, a full thermodynamic energy balance is
considered. Hence the simulator inherently takes into consideration multi-phase and
multi-component behavior and is suited for non-isothermal cases which allows the
handling of many physical phenomena including multiphase compositional flows with
phase transitions, different types of fluid-rock and fluid-fluid interactions (e.g. wettability
and adsorption), and various types of fluid rheology.
The DHD simulator is a research prototype optimized for high performance computing
(HPC) and applied to porous media systems. We demonstrate the utility of DHD to
simulate two-phase flow displacement ranging from the classical "Lenormand" pore-scale
displacement events and Roof's snap-off criteria to more complex cooperative
phenomena like capillary de-saturation and relative permeability. The presented
simulation results are benchmarked against experimental data in core flooding, 2D
micromodel, and synchrotron-based x-ray microtomography experiments and provide
good agreement.

1
 SCA2013-014 2/12

INTRODUCTION
The interest in enhanced oil recovery (EOR) methods is growing reflecting the gradual
maturation of conventional reserves. To date there is a wide range of EOR methods based
on gas injection, chemical injection, thermal treatment and even microbial EOR. EOR
requires an integrated workflow. It involves various direct and indirect enablers
interacting in well-executed synchronization from laboratory analysis screening to full
field development and monitoring. All EOR applications have in common a close
interconnection between multiphase hydrodynamics and physical chemistry in porous
media. In order to select the best suitable EOR method for a particular situation the
understanding and prediction of multiphase flow in various displacement and recovery
systems is required [1]. Traditionally core flooding studies are part of the EOR de-risking
process where comparative studies are used for process optimization. However, usually
there is an enormous space of design parameters to be tested for a particular EOR method
[1]. SCAL tests can be expensive and time consuming. Moreover laboratory capacities
are typically limited and usually not sufficient to assess the full range of uncertainties of
EOR processes. In addition these tests are typically destructive (i.e. there is always an
alteration of the core during the test, and already during coring and handling) which
means that tests cannot be repeated on the same rock sample under the same conditions.
Ultimately, for many EOR processes the performance is determined by the efficiency of
the associated pore-scale displacement processes. Direct pore scale simulation of these
processes, often termed digital rock (DR), can be an early-stage de-risking tool for
enhanced oil recovery workflows. The Digital Rock approach is a method for prediction
of transport properties of rock samples over a wide range of (reservoir) conditions, flow
regimes, and fluid composition in a 3D digital model of the rock of interest. It does not
suffer from some of the shortcomings of traditional core floods. Numerical modelling can
be repeated as many times as necessary on the same or various rock models, for different
scenarios of oil production or various EOR techniques etc. at a reduced time and cost.
However digital rock does not replace but rather complements laboratory experiments by
assessing sensitivities and reducing uncertainty in the design of, for instance, chemical
EOR processes [1]. Over the past three years, Shell and Schlumberger have jointly
developed a digital rock approach using a novel direct hydrodynamic simulator (DHD)
[2]. It is based on classical continuum hydrodynamics employing mass, momentum and
energy balance concepts, where a diffuse interface approach for the fluid-fluid interfaces
is implemented via the density functional approach.
A key requirement for any digital rock application is the correct modelling of two-phase
flow (without EOR agents) ranging from elementary pore scale displacement processes to
cooperative and macroscopic (Darcy scale) phenomena. In this paper we present an
overview and examples of an extensive validation study demonstrating that DHD does
correctly model elementary pore-scale processes like Haines jumps and snap-off,
cooperative pore-filling processes, up to complex macro-scale processes like capillary de-
saturation and the rate dependency of relative permeability. These are validated against

2
 SCA2013-014 3/12

reference cases in the classical literature, experimental data involving time-resolved 2D

micromodel and synchrotron-based micro-tomography [3].

DIGITAL ROCK WORKFLOW

The basic workflow of the digital rock approach has three basic steps.
Characterization. First the basic fluid, rock and fluid-rock properties are measured. This
involves PVT data, bulk fluid rheology, fluid-fluid properties (IFT), rock-fluid
interactions (surface energies, composition), and rock morphology.
Simulation. Based on the data acquired in the first step, digital models are built to
explore the dynamics of the rock fluid system in the porous media representative of the
reservoir under study.
Validation. Results of numerical simulations are validated against experimental results of
laboratory tests, which may include 2D micromodels with very well defined geometry
and conventional and special core analysis.

IMAGING AND IMAGE PROCESSING

X-ray microtomography (micro-CT) has been available for more than three decades and
has been a powerful tool for studying a wide array of processes in porous media systems
[4-7]. This imaging technique provides nondestructive and noninvasive 3dimensional
imaging of the interior of objects by mapping xray attenuation through a sample. Local
attenuation depends on the
physical composition of the
material and its mass density
which allows for the
examination of the internal
properties of a rock,
including microporosity and
pore morphology. An
example for a poorly
consolidated sandstone is
displayed in Figure 1 (left).
Typical micro-CT spatial
Figure 1. Greyscale and binarized microCT images of a poorly image resolution is between
consolidated sandstone. Voxel size is 2.15 m. 0.5 and 2 m (at 1-10 mm
field of view, respectively).
Typical acquisition time for benchtop machines is 2-24 hours, while faster scanning also
at dynamic flow conditions is possible at synchrotron beamlines [3, 8]. Within these
limitations, micro-CT scanning allows for resolving the pore space of sandstones with
features and length scales relevant to the confined fluid phases (i.e. oil and water). Clays
with micro-porosity and carbonates with multiple length scales still pose a challenge for
micro-CT scanning but the technology is rapidly developing: resolution of benchtop

3
 SCA2013-014 4/12

instruments has improved from 2-4 micrometers in 2006, to 50 nanometers in 2013 (for
field of view less than 1 mm) [9, 10]. This progress together with involvement of
different imaging techniques like high-resolution electron microscopy and confocal
microscopy [11] enables targeting more difficult rock types.
The digital rock model which is required for simulation of pore scale flow processes is
obtained by segmentation from grey-scale data obtained by micro-CT [12]. State-of-the
art image processing and segmentation workflows involve the application of edge
preserving filters to improve the signal-to-noise ratio (SNR), and complex methods like
active contour or watershed-based segmentation approaches [13] which are more robust
to reduce image artifacts and allow use of lower SNR data than global thresholding.
Moreover, such workflows are less sensitive to the subjectivity of the operator. An
alternative option is the application of geostatistics-based segmentation methods like
indicator Kriging [14] which make use of the noise structure of the data. The
segmentation result can be independently validated by comparison with measured
parameters, such as, porosity, pore size distribution and, more importantly, absolute
permeability (which can be directly computed on the pore space of the digital model by a
single-phase flow simulation with DHD). We observed a typical match of porosity and
permeability for sintered glass bead packs and sandstone rock within 10 and 15 %,
respectively. The permeability match might be more relevant for the simulation of flow-
related properties, because discrepancies in porosity might also be caused by micro-
porosity below imaging resolution which is also not contributing to flow.

DHD DESCRIPTION
The DHD simulator is based on the density functional (DF) method applied for
multiphase compositional hydrodynamics. The theoretical background of the DHD
simulator is explained in detail in [2]. It combines continuum fluid mechanics and
thermodynamic principles by considering mass, momentum and energy balance together
with a diffuse interface description. The diffuse interface approach is a physically
consistent and efficient way to model the evolution of the fluid-fluid interfaces in
multiphase flow. The DHD framework combines concepts from physical chemistry,
statistical physics and physics of solids with hydrodynamics and naturally takes into
account interfacial surface tension, interfacial tension at contact with solid surfaces
(wettability), moving contact lines and dynamic changes of topology of interfaces. DHD
is capable of modelling the following classes of multiphase hydrodynamic problems:
complex compositional fluids with phase transitions (gas-liquid, liquid-liquid, liquid-
solid); flow in complex geometries of boundary surfaces; wettability and adsorption;
surfactants, solvents, polymers; complex fluid rheology and presence of mobile solid
phase; and thermal effects. The fluid phase behavior, which is traditionally characterized
by an equation of state (EoS), has to be converted into a thermodynamic fluid model
specified by Helmholtz free energy functions used to model phase behavior in DHD.
High performance computing with a massively parallel GPU realization of the DHD code
together with enhanced algorithms of cross-machine and cross-GPU communications

4
 SCA2013-014 5/12

interleaved with computations, allows modelling several tens of billions (~1010) of cells
on a modern medium-sized GPU cluster. The HPC version of DHD code utilizes domain
decomposition parallelization method with asynchronous cross-machine
communications, which gives excellent scalability and nearly linear parallel speedup.
This guarantees practical simulation times on a todays GPU cluster systems. Currently,
characteristic computational times for complex multiphase flows in representative sub-
volumes of digitized rock samples are around 1 - 3 days, while simple geometries can be
tackled within several minutes.
The resolution of the current DHD-based DR approach goes down to about 500 nm and is
sufficient to discretize the typical pore scale (tens to hundreds of micrometers) to
accurately capture the thermodynamic and transport processes mentioned before.
However it is significantly above the thickness of interfaces and wetting layers that range
from 0.1 to 10s of nanometers [15-17]. The macroscopically relevant internal dynamics
of the interfaces are given as input to DHD and can be obtained from experiments or
higher resolution models (such as molecular dynamics [18-21], dissipative particle
dynamics [22-24]). Given the functionality and resolution, DHD is designed to provide
the right amount of rigor at the scale most appropriate for the EOR processes.

(a) (b)
(c)

Figure 2. (a) Piston-type flow from three ducts in imbibition scenario. (b): Piston-type flow from two
ducts in imbibition scenario: meniscus collapse. (c): Piston-type flow from two ducts in drainage scenario.
Wetting phase shown in blue and non-wetting phase is red semitransparent

In [2] a set of fundamental analytical benchmarks for DHD are documented which
include static and dynamic capillary phenomena like sessile drops, capillary bridges, and
processes which involve topological changes of interfaces like breakup and coalescence.

5
 SCA2013-014 6/12

In the following, we present a set of more specific benchmarks for multi-phase flow in
the visco-capillary regime in porous media.

EXAMPLES OF MULTIPHASE FLOW SIMULATIONS

Piston-type motion and meniscus collapse in two crossing ducts has been identified by
[25, 26] as a set of elementary processes in pore-scale displacement. When non-wetting
fluid displaces wetting fluid the flow occurs only when the pressure at the entrance of a
capillary equals or exceeds a threshold pressure determined by the capillarys geometry
and interfacial tension between wetting and non-wetting fluids. When a wetting fluid
displaces non-wetting fluid the flow starts spontaneously once the wetting fluid is in
contact with the capillarys entrance.

DHD reproduces the classical "Lenormand

displacements" very well as displayed in Figure 2
using the same flow parameters (geometry,
wettability and interfacial tension (IFT)) as in [25,
26].
In addition to the processes in Figure 2, two
different types of snap-off processes are part of the
elementary pore scale displacements:
Figure 3. DHD modeling of geometry
related snap-off Geometry related snap-off was identified by Roof
[27] as one of the key processes leading to
disconnection of the non-wetting phase. According
to Roof's model, snap-off happens when the radius
of the detaching drop is approximately twice the
radius of the non-wetting phase jet in the
constriction. That corresponds to a situation where
the capillary pressure in the constriction and in the
detaching droplet are equal. The DHD simulation
in circular capillaries with constrictions of
different aspect ratio (see an example in Figure 3)
shows good correspondence with Roofs criteria .
In the second type of snap-off the non-wetting
phase is disconnected through an alternative flow
path for the wetting phase with higher capillary
pressure like in pore doublets. Figure 4 illustrates
Figure 4. DHD modeling of imbibition in an example of imbibition phenomena in pore
pore doublets doublets of different radii. DHD results are
coherent with experimental data from [28]
Microfluidic devices: modeling vs. experiment

6
 SCA2013-014 7/12

While in the previous examples the focus was on spatial distribution of fluid in
displacement processes, in the next step DHD's ability to correctly model the dynamics is
probed. In addition to the shapes of interfaces, in particular the time scale of filling events
is tested. For that purpose, drainage experiments in 2D micromodels were imaged with a
high-speed camera at 2000 frames per second and then compared to DHD simulations in
a similar geometry (Figure 5). We find the correct time scale for drainage events (~12
ms) is also captured with DHD demonstrating that DHD also correctly resolves the fast
dynamics of interfaces during the pore filling events. Note that the time scale and
dynamics of pore drainage events ultimately affect front propagation and fluid phase
topology during fluid displacement [29-31]. Thus, accurate simulation of the time scales
and dynamics of pore filing events is critical for any digital rock simulation platform.
Details of the experiment will be published separately [31].

Figure 5. Drainage dynamics in micromodels.

The next level of complexity in the two-phase flow validation are cooperative processes
with direct macroscopic relevance.
Capillary Desaturation is the process where the non-wetting phase is mobilized by
increasing the balance of (mobilizing) viscous forces over the (trapping) capillary forces,
i.e. the capillary number. Capillary de-saturation is the basis of many EOR processes like
surfactant flooding. For the validation of DHD to literature data [32] simulations are
performed on a digitized water-wet sandstone sample with a porosity equal to 0.23 and an
absolute permeability of 1150 mD. Following the experimental workflow in [32] the
model was initially filled with 100% of the non-wetting oil. Subsequently it was flooded
by the wetting phase (water) as shown in Figure 6. Water injection was modeled at
different capillary numbers (Ca) and continued until the oil fraction in the produced fluid
dropped below 2%. The capillary desaturation curve (Figure 7) calculated by DHD is in
very good correspondence with literature data on sandstone rock [32]. Note that these
were not identical samples which explains the small differences. A more detailed
discussion on capillary desaturation will be published separately [33].
Relative Permeability

7
 SCA2013-014 8/12

In addition to the residual oil saturation also macroscopic flow parameters like relative
permeability show a dependence on capillary number [34]. Results of DHD simulations
are displayed in Figure 8 where a steady-state scenario was modelled. Initial conditions
were obtained by distributing the phases in accordance with given saturations and local
wettability of the digital rock model. The wettability of the solid surfaces is specified by

Figure 6. DHD modeling of oil (red) recovery by water (blue) flooding

Figure 7. Capillary desaturation curve for the Figure 8. Relative permeability curves for the digitized
digital rock model of the Berea sandstone core sandstone sample at different flow regimes.
sample. Dashed and dotted curves for
comparison are from [32]
means of the specialized boundary conditions used in our method [2]. According to the
DHD framework each cell of the model (each voxel of the digital rock model) can be
assigned with individual wettability respecting the local mineralogical content. In the
particular simulation we have used a stochastic wettability model representing the
existing knowledge about the nature of the sample. Thus, the wettability distribution was
created in such a way that the model was made generally water-wet with some regions
wet better than the others. Before running the relative permeability simulations, the fluids
were distributed in accordance with specific saturations and then a preliminary
simulations were carried out on a closed model (without external fluxes) to let fluids

8
 SCA2013-014 9/12

reach thermodynamic equilibrium in respect to the local wetting properties. The flow was
imposed by applying a body force in one spatial direction. The dependence of relative
permeability curves on Ca shows very similar trends as typical core flooding data [34].
Pore filling events: comparison with synchrotron data

Figure 9. Pore filling events. Red color represents oil phase. Water phase and grains are set to be
transparent. Left column is experimental data showing oil distribution within two different ROIs before
pore filling event. Middle column is experimental data showing oil distribution after pore filling event.
Right column show results of DHD modeling: oil distribution after pore filing event. Numerical modelling
of the pore filling event used distributions from the left as initial condition.

The last example presented here is the direct modelling of selected regions of interest of a
core flooding experiment where the pore scale displacements have been imaged at real-
time, i.e. under dynamic flow conditions, by using fast synchrotron based micro-
tomography. Drainage and imbibition experiments performed on a 4 mm water wet Berea
sandstone sample were imaged at a time resolution of a few seconds to capture the
capillary equilibrium states between pore-filling events that last on the order of a
millisecond [3, 8]. Here we present the modeling of pore-filling events from
characteristic regions of interest (ROIs) during drainage (see Figure 9). For the
simulations we used exactly the same fluid properties, pore geometries and flow regimes
as in the experiment. The experimental fluid distribution before the pore filling event was
used as initial condition for simulating the dynamics of the events. One can see that
experimental and numerical geometries of the oil phase after the pore filling events are

9
 SCA2013-014 10/12

very similar. Also the displaced volumes are in very good agreement: for the
displacement event in the top panel of Figure 9 we obtain 5.9 nL for the experiment and
5.92 nL for the DHD simulation. Differences between simulation and experiment are
attributed to the fact that only a relatively small field of view of a larger experiment was
modeled for which it is difficult to define the correct flow/pressure inflow boundary
conditions [35].

SUMMARY
The DHD-based digital rock approach has been extensively benchmarked versus
experimental data for two-phase flow displacement processes with respect to microscopic
details, dynamics and macroscopic cooperative effects. In summary, excellent agreement
has been found which ultimately establishes the reliability of DHD for immiscible
displacement processes. This opens new perspectives for applications to advanced SCAL
characterization, such as full cycle drainage-imbibition relative permeability simulation
and chemical EOR de-risking. The foreseeable growth of computational power will allow
more practical and affordable simulation times and access larger sample sizes, which will
enable transition of the digital rock techniques from an R&D tool to practical commercial
implementation in the near future.

ACKNOWLEDGMENT
We thank Shell and Schlumberger for permission to publish this work. We acknowledge
TomCat beam line at PSI for active participation in conducting the synchrotron core
flooding experiments. Mark Andersen, Dmitry Korobkov and Bart Suijkerbuijk are
acknowledged for helpful discussions and reviewing the manuscript.

REFERENCES
1. Koroteev D. et al., Application of Digital Rock Technology for Chemical EOR Screening (2013)
SPE-165258 to be available Jul 2013
2. Demianov A., Dinariev O. and Evseev N., Density functional modelling in multiphase compositional
hydrodynamics, Can. J. Chem. Eng. (2011), 89, 206 226
3. Berg S. et al., Multiphase Flow in Porous Rock imaged under dynamic flow conditions with fast X-
ray computed micro-tomography, (2013) SCA 2013-011.
4. Wildenschild D., et al. Quantitative analysis of flow processes in a sand using synchrotronbased
xray microtomography, Vadose Zone Journal (2005), 4, 1, 122126.
5. Wildenschild D., et al., Using xray computed tomography in hydrology: systems, resolutions and
limitations, Journal of Hydrology (2002), 267, 285297.
6. AlRaoush, R. I. Impact of wettability on porescale characteristics of residual nonaqueous phase
liquids, Environmental Science and Technology (2009), 43, 13, 47964801.

10
 SCA2013-014 11/12

7. Porter, M. L et al. Measurement and prediction of the relationship between capillary pressure,
saturation, and interfacial area in a NAPL water glass bead system, Water Resource Research (2010),
46, W08512.
8. Berg, S., et al., Real-time 3D imaging of Haines jumps in porous media flow, Proceedings of the
National Academy of Sciences (PNAS) (2013), 110, 10, 37553759.
9. SkyScan 2011 X-ray nanotomograph with a stated resolution of 150-200 nm
10. XRadia UltraXRM-L200 with a stated resolution of 50 nm
11. Ibrahim et al., An Automated Petrographic Image Analysis System: Capillary Pressure Curves Using
Confocal Microscopy (2010), SPE-159180
12. Coles, M., et al., Developments in synchrotron x-ray microtomography with applications to flow in
porous media, SPE Res Eval & Eng (1998), 1, 4, 288296.
13. Sheppart, A. P., Sok, R. M. and Averdunk, H., Techniques for image enhancement and segmentation
of tomographic images of porous materials, Phyisca A (2004) 339, 145-151.
14. Oh W., Lindquist W., Image thresholding by indicator kriging IEEE Trans. Pattern Anal. Mach.
Intell. (1999) 21, 590.
15. Davis, H. T., Scriven, L. E., Stress and Structure in Fluid Interfaces, Advance in Chemical Physics
(1982) XLIX, 357-454.
16. Evans, R., The nature of the liquid-vapour interface and other topics in the statistical mechanics of
non-uniform, classical fluids, Advances in Physics (1979) 28, 2, 143-200.
17. Hirasaki, G. J., Wettability: Fundamentals and Surface Forces, SPE Formation Evaluation (1991),
217-226.
18. Ashurst, W. T., Hoover, W. G., Dense-fluid shear viscosity via nonequilibrium molecular dynamics,
Physical Review A (1975), 11, 2, 658-678.
19. Frenkel D., Smit B, Understanding Molecular Simulation: from algorithms to applications. San Diego,
California: Academic Press (2001).
20. Allen, M. P., Tildesley, D. J., Computer Simulation of Liquids, Clarendon Press, Oxford (1987).
21. Werder T., et.al, On the Water-Carbon Interaction for Use in Molecular Dynamics Simulations of
Graphite and Carbon Nanotubes, J. Phys. Chem. B (2003), 107, 1345-1352.
22. Hoogerbrugge, P. J. and Koelman, J. M. V. A., Simulating Microscopic Hydrodynamic Phenomena
with Dissipative Particle Dynamics, Europhys. Lett. (1992) 19, 155.
23. Buijse M, et al, Surfactant Optimization for EOR using Advanced Chemical Computational
Methods, (2012) SPE 154212
24. Liu M., Meakin P., Huan H., Dissipative particle dynamics simulation of pore-scale multiphase flow,
Water Resources Research (2007) 43, W04411.
25. Lenormand R., Zarcone, C. and Sarr, A., Mechanisms of the displacement of one fluid by another in
a network of capillary ducts J. Fluid Mech. (1983) 135, 337353.
26. Lenormand R., Liquids in porous media, J. Phys.: Condens. Matter (1990) 2, SA79SA88.
27. Roof, J.R. Snap-off of oil droplets in water-wet pores, SPE Journal (1970) 10, 1, 8590.
28. Chatzis I. and Dullien, F.A.L., Dynamic immiscible displacement mechanisms in pore doublets:
Theory versus experiment, J. Coll. Inter. Sci. (1983) 91, 1, 199222.
29. Moebius, F. and D. Or, Interfacial jumps and pressure bursts during fluid displacement in interacting
irregular capillaries, Journal of Colloids and Interface Science (2012), 377, 406-415.
30. Mohanty, K. K., Davis, H. T., and L. E. Scriven, Physics of oil entrapment in water-wet rock, SPE
Reservoir Engineering (1987), SPE 9406, 113-128.
31. Armstrong R. T. and S. Berg Interfacial velocities and capillary pressure gradients during Haines
jumps (2013 in progress).

11
 SCA2013-014 12/12

32. Lake L. W., Enhanced Oil Recovery, Prentice Hall (1989).

33. Evseev N. et al, Capillary desaturation modelling with density functional hydrodynamics (2013 in
progress)
34. Skauge, A., Thorsen, T., and Sylte A., Rate selection for waterflooding of intermediate wet cores,
(2001) SCA 2001-20.
35. Porter Mark L., et al., Lattice-Boltzmann simulations of the capillary pressuresaturationinterfacial
area relationship for porous media, Advances in Water Resources (2009), 32, 16321640.

12