The effect of a microscale fracture on dynamic capillary pressure of two-phase flow in porous media

Accepted Manuscript

The effect of a microscale fracture on dynamic capillary pressure of two-phase flow in porous media Mingming Tang , Huifang Ma , Shuangfang Lu , Hongbin Zhan , Wenyue Guo PII: DOI: Reference:

S0309-1708(17)30931-4 10.1016/j.advwatres.2018.01.015 ADWR 3071

To appear in:

Advances in Water Resources

Received date: Revised date: Accepted date:

28 September 2017 15 January 2018 15 January 2018

Please cite this article as: Mingming Tang , Huifang Ma , Shuangfang Lu , Hongbin Zhan , Wenyue Guo , The effect of a microscale fracture on dynamic capillary pressure of two-phase flow in porous media, Advances in Water Resources (2018), doi: 10.1016/j.advwatres.2018.01.015

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights  Lattice Boltzmann method is used to invest dynamic capillary pressure (DCP) effects Significant effects of micro scale fractures on DCP are observed



The orientation of micro-fracture will influence DCP coefficient

AC

CE

PT

ED

M

AN US

CR IP T



ACCEPTED MANUSCRIPT

The effect of a microscale fracture on dynamic capillary pressure of twophase flow in porous media

1

CR IP T

Mingming Tang1*, Huifang Ma2# , Shuangfang Lu1* , Hongbin Zhan3* and Wenyue Guo2

Research Institute of Unconventional Oil & Gas and Renewable Energy, China University of Petroleum, Qingdao, China, 266580; [email protected]; [email protected]

College of Science, China University of Petroleum, Qingdao, China, 266580;

AN US

2

[email protected]; [email protected] 3

Department of Geology & Geophysics, Texas A&M University, College Station, TX

PT

ED

M

77843-3115; [email protected]

* Corresponding authors:

CE

(Mingming Tang) [email protected]; Tel.:+86-13156891161; (Hongbin Zhan) [email protected]; Tel.: +1-979-862-7961;

AC

(Shuangfang Lu) [email protected]; Tel.:+86-0532-86983191 #

Huifang Ma is co-first author

2

ACCEPTED MANUSCRIPT

Abstract: Dynamic capillary pressure (DCP) effects, which is vital for predicting multiphase flow behavior in porous media, refers to the injection rate dependence capillary pressure observed during non-equilibrium displacement experiments. However, a clear picture of the effects of microscale fractures on DCP remains elusive. This study quantified the effects of microscale fractures on DCP and simulated pore-scale force and saturation

CR IP T

change in fractured porous media using the multiphase lattice Boltzmann method (LBM). Eighteen simulation cases were carried out to calculate DCP as a function of wetting phase saturation. The effects of viscosity ratio and fracture orientation, aperture and length on

AN US

DCP and DCP coefficient were investigated, where refers to the ratio of the difference of DCP and static capillary pressure (SCP) over the rate of wetting-phase saturation change versus time. Significant differences in values were observed between unfractured and fractured porous media. The values of fractured porous media were , which were one or two orders of magnitude lower than those of

M

5.68

to

. A horizontal fracture had greater

ED

unfractured porous media with a value of 4

effects on DCP and than a vertical fracture, given the same fracture aperture and length.

PT

This study suggested that a microscale fracture might result in large magnitude changes in

CE

DCP for two-phase flow.

Keywords: fracture; porous media; dynamic capillary pressure; multiphase flow; lattice

AC

Boltzmann method

3

ACCEPTED MANUSCRIPT

1. Introduction Understanding and predicting multiphase flow behavior in fractured porous media is a challenging subject in a number of disciplines and applications, such as CO2 storage in brine aquifers (Berg and Ott, 2012, Hall et al., 2016, Juanes et al., 2006, Kim and Santamarina,

CR IP T

2016, Pruess and Garcia, 2002, Qi et al., 2009), oil recovery in fractured reservoirs (Tang et al., 2017) and water flow back in shale gas development (Cao et al., 2017). Recently, several field and laboratory experiments and numerical analysis have been conducted to investigate the nonlinear transient and spatial characteristics of multiphase flow in fractured porous

AN US

media (Tang et al., 2017, Faybishenko, 2017). Fractures not only dramatically change the permeability of the porous media but also change the effective/average static capillary pressure. In addition, fractures usually show unique multiscale characteristics, whose main dimensions could span several orders of magnitude in space with numerous microscale

M

fractures and large faults (Zhang et al., 2017c).

ED

Since it is impractical to represent every microscale fracture in a numerical model, due to the enormous computational cost associated with a complete discretization of these fractures

PT

(Jerbi et al., 2017), matrix and microscale fractures are usually divided into small

CE

homogenized continuum blocks in field-scale simulations (Tang et al., 2017), such as equivalent continuum models (Wu, 1999, Tang et al., 2017), dual-domain or dual-porosity

AC

continuum models (Jerbi et al., 2017, Lichtner, 2000), and their various extensions (Pruess and Garcia, 2002). In these models, the mass balance and Darcy‟s law are assumed to be valid for each block (Tang et al., 2017, Landry et al., 2011). For two-phase flow, which is the primary concern of this study, the governing equation can be written as (Dahle et al., 2005, Yu et al., 2015, Tang et al., 2017, Sussman et al., 1998): (1)

4

ACCEPTED MANUSCRIPT

is the viscosity of phase , which can be a wetting phase (w) or a non-wetting

phase (n),

is the porosity,

for phase ,

is the absolute permeability,

is the saturation of phase

pressure of phase , t is time, and equation among

,

, and

bounded by a constraint of

,

is the

is the differential operator in space. The constitutive

, is the capillary pressure equation:

, where

is the relative permeability

(2)

CR IP T

where

is usually called the dynamic capillary pressure (DCP) (Hou et al., 2012) or

AN US

non-equilibrium fluid pressure difference (Bottero et al., 2011), which is dependent on actual flow dynamics, the subscript “c” means capillary pressure hereinafter, and the superscript “dyn” means dynamic hereinafter.

Different from DCP, another widely-used capillary pressure called the static capillary

M

pressure (SCP) or equilibrium fluid pressure difference,

, which is often measured under

ED

a static or quasi-static state in laboratory experiments (Bottero et al., 2011, Abidoye and Das, 2014, Das et al., 2005), is frequently used as a substitute of

in field-scale

PT

simulations (Tang et al., 2017, Yu et al., 2015). It is notable that the terms of non-equilibrium and dynamic have been used exchangeable in many previous studies, and the terms of

CE

equilibrium and static have also been used exchangeable by various investigators.

AC

Hassanizadeh et al. (2005), Bottero et al. (2011) gave a detailed discussion about the terms of non-equilibrium and dynamic and the terms of equilibrium and static in respect to the capillary pressure. In this study, the non-equilibrium state refers to a flow field that is still dynamically change with time, whereas the equilibrium state refers to a flow field that is not change with time anymore. However, many studies (Abidoye and Das, 2014, Das and Mirzaei, 2012, Niasar et al., 2010, Goel et al., 2016) have observed unique features of DCP that has not been seen in SCP 5

ACCEPTED MANUSCRIPT during non-equilibrium two-phase flow in porous media, suggesting that the conventional SCP concept may not accurately describe two-phase flow behavior under non-equilibrium flow conditions. Several quantitative equations have been proposed to describe DCP, considering the change in saturation rate based on experimental results (Bottero et al., 2011, Abidoye and

CR IP T

Das, 2014, Goel et al., 2016). A common practice of those studies was to relate DCP and SCP with a constant coefficient, and a widely used equation was proposed by Hassanizadeh et al. (2005):

where

AN US

(3)

denotes the difference of DCP and SCP, is called the DCP coefficient hereinafter ,

and Sw is the wetting-phase saturation in a two-phase flow problem

The value of is crucial for studying two-phase flow in porous media (Bottero et al.,

M

2011, Dahle et al., 2005, Hassanizadeh et al., 2005). However, up to date, only a number of values of have been published from experimental and modeling studies using unfractured

ED

porous media (Hou et al., 2012, Bottero et al., 2011, Hassanizadeh et al., 2005, Abidoye and

PT

Das, 2014, Das and Mirzaei, 2012, Goel et al., 2016). A comprehensive investigation of DCP in fractured porous media remain elusive.

CE

The objective of this study is to fill this knowledge gap. Specifically, we will investigate the effects of fracture characteristics, such as fracture orientation, aperture and length on DCP

AC

and , and will report innovative methods to determine the value. In addition, we will investigate the effects of viscosity ratio on DCP and

A realistic porous media model is

developed based on X-ray micro-tomography (CT) image data of a sandstone sample, which is obtained from a borehole drilled to the depth of 1900m, from well Long-26 in Fuyu Formation, Songliao Basin, Northeast of China. The effects of injection rates and fluid mobility ratio (which is defined as the ratio of the viscosities of wetting phase and non6

ACCEPTED MANUSCRIPT wetting phase) are studied. A simplified artificial microscale fracture is added into the porous model to analyze the influence of fracture on two-phase flow, in particularly, to investigate DCP and its associated values. The lattice Boltzmann method (LBM) is used to simulate two-phase flow in the porous model in this study.

CR IP T

2. Methods Many numerical methods have been proposed to simulate multiphase flow in porous media (Cahn and Hilliard, 1958, Badalassi et al., 2003, Hirt and Nichols, 1981, Pilliod and Puckett, 2004, Osher and Sethian, 1988, Sussman et al., 1998, Sukop and Thorne, 2006, Bao

AN US

and Schaefer, 2013, Tang et al., 2015, Chen and Doolen, 1998). Among these methods, LBM has been proven to be an efficient and effective simulation method for a variety of problems (Liu et al., 2012, Gunstensen et al., 1991, Shan and Chen, 1993, Swift et al., 1995, He et al., 2000, Nourgaliev et al., 2003, Bao and Schaefer, 2013, Sukop and Thorne, 2006, Succi et al.,

M

2014, Cahn and Hilliard, 1958). Eshraghi et al. (2012) reported that the LBM model could be

ED

used on simulation of dendritic solidification with 36 billion grid points, and the parallel efficient is almost 100% by using central processing unit (CPU) cores of 1000~4068. Li et al.

PT

(2016) reported that the three-dimensional multi-phase LBM simulation cases (with a model ) could achieve a parallel efficiency of about 60% with 400K CPU

CE

size of

cores using a hybrid and heterogeneous, multiple data programming model.

AC

Up to present, there are at least three main three-dimensional LBMs (Suga et al., 2015), namely the D3Q15 model, the D3Q19 model, and the D3Q27 model. After a careful check of these three LBMs, we select the D3Q27 model (Figure 1) to simulate two-phase flow in the CT image-based digital rock model in this study, with a specific focus on analyzing the effects of a horizontal or vertical fracture on DCP and . Such a choice is based on two considerations. Firstly, Kang and Hassan (2013) found that the D3Q19 model violated Galilean invariance while the D3Q27 model did not. Secondly, Kuwata and Suga (2015) and 7

ACCEPTED MANUSCRIPT Suga et al. (2015) reported that both the D3Q19 model and the D3Q15 model produced false forces (Reynolds stresses) in particular directions while such an issue was negligible in the D3Q27 model. Therefore, the D3Q27 model is a proper choice for the type of problem investigated here. Figure 1 shows the lattice cell of the D3Q27 model, in which the number notes are the

CR IP T

discrete speed vectors. There are 27 discrete speed vectors in total. The 0th discrete speed vector is placed at the center of the cell with a discrete velocity value of zero, and the other discrete speed vectors are pointed from the center of the cell to the center points of faces (16) with a discrete velocity value of one, the middle points of edges (7-18) with a discrete

AN US

velocity value of √ , and the corners (19-26) with a discrete velocity value of √ , all in

CE

PT

ED

M

dimensionless forms (Liu et al., 2016, Suga et al., 2015).

AC

Figure 1. The lattice cell of D3Q27 model. In the D3Q27 model, the distribution function satisfies the following lattice Boltzmann

equation (Sukop and Thorne, 2006, Bao and Schaefer, 2013, Tang et al., 2015): ⃗⃗⃗⃗

(4)

8

ACCEPTED MANUSCRIPT where k is the discrete speed vector serial number from 0-26 (see Figure 1),

is the

microscopic distribution function of molecular in the direction of discrete speed vector ⃗⃗⃗⃗ for is the time step, “eq” in the superscript of

phase

state hereinafter,

represents the equilibrium

in the superscript (or subscript) means the phase

hereinafter,

is the

relaxation time which expresses the fact that molecular collisions tend to relax the

CR IP T

distribution function toward the equilibrium state (Bhatnagar et al., 1954), and its determination will be provided later in this section.

Eq. (4) is derived from the BGK equation (Chen et al., 1991, Sukop and Thorne, 2006,

⃗⃗⃗⃗

AN US

Bao and Schaefer, 2013, Chen and Doolen, 1998):

(5)

The fluid viscosity

is recovered as

M

(6)

where cs  1 / 3 is the sound speed in lattice units. at direction of ⃗⃗⃗⃗ can be

ED

The microscopic equilibrium distribution function

PT

written as (Sukop and Thorne, 2006, Succi et al., 2014):

where

(⃗⃗⃗⃗ ⃗⃗⃗⃗⃗

)

|⃗⃗⃗⃗⃗

|

(7)

CE

⃗⃗⃗⃗ ⃗⃗⃗⃗⃗

are weight factors, . is the sign of magnitude value of a vector,

is the fluid

AC

density, which could be obtained from: ∑

(8)

and the macroscopic fluid velocity ⃗⃗⃗⃗

is given by (Sukop and Thorne, 2006, Succi et

al., 2014): ⃗⃗⃗⃗

∑

⃗⃗⃗⃗

(9)

9

ACCEPTED MANUSCRIPT It is notable that the microscopic equilibrium distribution function in Eq. (7) is determined using the macroscopic fluid properties of Eqs. (8)-(9) (Sukop and Thorne, 2006, Succi et al., 2014). Using macroscopic properties to determine the microscopic equilibrium distribution function has commonly done in many previous LBM studies (Mohamad, 2011, Succi et al., 2014, Suga et al., 2015, Liu et al., 2016).

CR IP T

, the discrete speed vectors ⃗⃗⃗⃗ , and sound speed used

Table 1 lists the weight factors in the D3Q27 model.

Table 1. Parameters used in D3Q27 model ⃗⃗⃗⃗

Model (0,0,0)

√ ( (

(

), (0, ),( , ,

AN US

D3Q27

0), (0,0, ) ,0 ) (0, , )

)

M

Based on the D3Q27 model, we employ the Shan–Chen pseudo-potential concept in the

ED

present work due to its simplicity and flexibility (Shan and Chen, 1993, Huang et al., 2007). The basic idea of this concept is to use a pseudo-potential that depends on the local density

PT

and interaction strength to represent the pairwise microscopic molecular interactions at the macroscopic scale (Huang et al., 2007). Once the interaction strength is properly chosen, the

CE

phases will spontaneously segregate (Chen and Doolen, 1998, Huang et al., 2007). This

AC

concept represents the microscopic physics closely, and is flexible in dealing with the surface tension force (Tang et al., 2015, Tang et al., 2017). Because of its remarkable simplicity, flexibility and clear representation of the underlying microscopic physics of phase segregation, the pseudo-potential concept has been successfully used in a broad range of studies about multiphase flow (Shan and Chen, 1993, Landry et al., 2014, Huang et al., 2007), especially for multiphase flow in porous media (Landry et al., 2014, Landry et al., 2011, Tang et al., 2015). 10

ACCEPTED MANUSCRIPT According to this model, the total force acting on the particles of fluid phase , denoted as

, can be written as (Succi et al., 2014): ∑̅

where ̅,

∑

̅

⃗⃗⃗⃗

̅

is the pseudo-potential,

̅

⃗⃗⃗⃗

(10)

is the interaction strength between phases

and

is the number of phases, and b is 26 for the two-phase flow of this study.

CR IP T

Based on the interaction force, there are several implementations of LBM for multiphase flow. Our implementation is achieved through a common velocity of ⃗ ∑

⃗

∑

⃗⃗⃗⃗

∑

∑

(Huang et al., 2007, Succi et al.,

2014). The transient velocity of each phase ⃗⃗⃗⃗

⃗⃗⃗⃗

modified by the interaction force

AN US

be written as (Huang et al., 2007): ⃗

, defined as

̅

(11) is a constant G in all

M

Since we are only considering two-phase flow in this study,

can

simulation cases. Additional fluid phases can be implemented by changing

̅

into an

ED

interaction strength matrix.

PT

Contact angle is also controlled by interaction forces between the solid and fluid phases,

CE

which can be easily implemented by replacing potential

AC

solid and fluid

⃗⃗⃗⃗

, and replacing

̅

̅

⃗⃗⃗⃗

in Eq. (10) with a solid

in Eq. (10) with the interaction strength between

(Bao and Schaefer, 2013, Tang et al., 2015), where

and

in the

superscript represent the fluid phase and solid phase, respectively, and “a” in the subscript is an abbreviation of adhesive force between the fluid and solid phases. The solid surface force between fluid phase

and solid wall ∑

could be calculated as: ⃗⃗⃗⃗

⃗⃗⃗⃗

11

(12)

ACCEPTED MANUSCRIPT

where

is the interaction strength between solid wall

⃗⃗⃗⃗

and fluid phase ,

is equal to 1 for solid wall and 0 for fluid. The total interaction force acting on the fluid node is:

. The dimension analog is applied in all sixteen simulation cases of this study to convert

CR IP T

the simulation data from lattice units into physical units. For velocity, kinematic viscosity, relaxation time and surface tension coefficient, the following equations are satisfied (Frish et al., 1987, Gunstensen et al., 1991, Shan and Chen, 1993, Liu et al., 2016)

̅

(14) (15)

(16)

M

̅

AN US

(13)

where the superscripts LB and p represent the units in lattice space and physical space, is the surface tension coefficient in lattice units, and is determined by the

ED

respectively,

and

and

PT

interaction strength G, the space step

is determined by the resolution of a porous model,

are the characteristic velocities in lattice space and physical space,

CE

respectively. The time step

is determined by Eq. (13), the relaxation time

determined by Eq. (15) based on viscosity

is

, ̅ in Eq. (16) is a reference density

(1000kg/m3 in this study). The space step

AC

and time step

in the lattice units are all

kept as unity, similar to what has been done in many previous LBM-based multiphase flow studies (Hao and Cheng, 2010, Frish et al., 1987, Liu et al., 2012, Bao and Schaefer, 2013).

12

CR IP T

ACCEPTED MANUSCRIPT

AN US

Figure 2. Schematic of pressure distribution in a single capillary tube containing two fluids with equal viscosity separated by an interface at x(t). A and F are respectively the inlet and outlet of the capillary tube, C and D are respectively the points of the non-wetting and wetting phases at the interface of x(t). B is the middle point of segment AC, and E is the middle point of segment DF. and are the pressure drop caused by viscous force in non-wetting phase and wetting phase, respectively.

Another important subject related to capillary pressure is upscaling. To the best of our knowledge, there are at least three main types of upscaling methods for capillary pressure (Miri et al., 2014), volume averaging

M

(Figure 2), including boundary averaging type type

(Dahle et al., 2005), and energy conservation law-based type

(Raeini et al.,

ED

2014b, Raeini et al., 2014a). Figure 2 is a schematic diagram to interpret the differences

experiment,

PT

among these three capillary pressure upscaling methods. In general, during a displacement .

CE

The boundary averaging type

is widely used in experiment-based capillary pressure

AC

studies, since it is convenient to measure the average pressure pressure

at the outlet in laboratory, thus

. The volume averaging type

at the inlet and the average

can be easily calculated by

is more popular in analytical and numerical studies

because the pressure at every point can be determined precisely. However, usually larger than the physically-based capillary pressure experiment since of

and

and

and

are

during a displacement

(Figure 2) are also used to calculate the pressure difference

. 13

ACCEPTED MANUSCRIPT The energy conservation law-based type

is also calculated based on numerical

simulations. Applying the energy conservation law to porous media (Raeini et al., 2014b, Raeini et al., 2014a), one can obtain the following balance equation for energy rate (17)

energy rate of

is the energy rate of capillary pressure,

∑

,

phase and non-wetting phase,

is the total energy rate of viscous force in wetting is the energy rate of viscous force in phase

AN US

calculated by: ( ∑

(

,

and

at

are respectively the pressure and volume of grid cell (i, j, k) at grid cell (i, j, k), and is treated as a constant

PT

is the viscosity of phase

CE

in this study.

)

are the simulated x, y, and z components of velocities of phase

grid cell (i, j, k), with phase ,

)

18)

M

,

)

ED

where

, which is

)

( (

is the

CR IP T

where Q is discharge at the inlet,

3. Data preparation

AC

Figure 3-a shows the realistic digital porous model developed based on CT image data

from ZEISS Xradia CT scanning with a voxel resolution of 0.1

. The total size of the

digital porous model is 1700 × 800 × 400, and the porosity is 22.4%. The mean pore radius is 1.6

and the permeability is

. To investigate the effects of fracture on DCP

and its associated DCP coefficient of , an idealized single fracture consisting of two parallel planes is added into above porous models with different orientations, apertures and lengths as 14

ACCEPTED MANUSCRIPT shown in Figure 3 and Table 2. The grid cells in original porous media (Figure 3-a) are divided into two types as done in the work of Sukop and Thorne (2006). One type of grid cell is called solid grid cell which represent the solid, the left grid cells are called fluid grid cells which represent the space where the fluids could flow through. In our simulations, the type of each grid cell is represented by a Boolean number B. The solid grid cells are marked as B=1,

CR IP T

and the fluid grid cells are marked as B=0. To ensuring the connecting of the adjacent pore space and fracture space, we implement the idealized single fracture by assigning all the grid cell type number B to 0 between the two parallel fracture planes, meaning that the space

AC

CE

PT

ED

M

AN US

between the two fracture planes can only be occupied by fluids, not by solids (see Figure 3).

Figure 3. Three dimensional models: (a) unfractured porous model; (b) porous model with a horizontal fracture whose aperture is 2 and length is 170 ; (c) porous model with a vertical fracture whose aperture is and length is 170 ; (d) porous model with a horizontal fracture whose aperture is 4 and length is 170 ; (e) porous model with a horizontal fracture whose aperture is 1 and length is 80 ; (f) porous model with a horizontal fracture whose aperture is 1 and length is 150 ; (g) porous model with a horizontal fracture whose aperture is 2 and length is 80 ; (h) unfractured porous model

15

ACCEPTED MANUSCRIPT with a larger porosity of 34% comparing with (a) with a porosity of 23.40% after eroding process (The scales of the model are: Lx=170 , Ly=80 , Lz=40 ). To investigate the effects of fracture length on DCP, we built two models shown in Figures 3-e and 3-f, with different fracture lengths of 150

and 80

, respectively. We

also built a model, shown in Figure 3-g, with the same fracture length and aperture as the

CR IP T

model in Figure 3-c. One additional porous model with a much larger porosity of 34% is also built using an eroding process carried out using the IMERODE function in MATLAB with an eroding radius of 3.

For the two-phase displacement simulations using above eight porous models, the

AN US

densities of displacing and displaced fluids are set as unity because the gravitational force is usually negligible compared to the capillary and viscous forces in our case, however it must be noted that gravitational forces can be very significant for chlorinated solvents. The contact

M

angle is set to be 140o, the same as the contact angle used in the conventional mercury injection capillary pressure experiments (Zhang et al., 2017a), which is a standard method to

ED

measure SCP. The viscosities are all set to be 2mPa.s for the sake of simplicity. Different

PT

viscosities can be used if necessary. Ultra-low flux rates are used in the simulation cases, and the max Mach number is around Ma=0.0029, which is suitable for studying incompressible

CE

transient flow problems based on LBM, according to the work of Hazi (2003) and Reider and Sterling (1995).

AC

To minimize the boundary effects, one inlet zone of 400 layers and one outlet zone of

30 layers are added to the model (Figure 3). We use more layers in the inlet zone than the outlet zone to implement the constant injection velocity boundary before the non-wetting phase reaching the porous zone. For obtaining the high-resolution simulation results, the number of elements for the porous model is as large as 2130 × 800 × 400. The Ada highperformance computer cluster at the Texas A&M University (TAMU) at College Station

16

ACCEPTED MANUSCRIPT campus of USA is used to simulate these models. Each simulation case requires approximately 3000 CPUs for about 10.5 days. 4. Results and discussions Eighteen simulation cases associated with above eight porous models (Tables 2 and 3)

CR IP T

are considered to investigate the effects of a horizontal or vertical fracture on DCP and DCP coefficient . In each simulation case, the porous media is first saturated with a wetting phase, and then a non-wetting phase is injected from the inlet side to displace the wetting phase. The contact angle is set as 140o for all simulation cases to eliminate

AN US

the influence of wettability. The value of G which is the interaction strength for twophase flow, is set as 2.0 to implement surface tension coefficient

=25mN/m

according standard processes used in LBM simulations (Shan and Chen, 1993, Huang et

M

al., 2007). Two injection fluxes, 1mm/s and 1.5mm/s, are applied at the inlet boundaries to observe flow rate-dependent effects on DCP. Continuous boundary

ED

condition is applied to the outlet boundaries. Non-slip wall boundaries are applied to other boundaries. The space step and time step are respectively

PT

=1× 10-9s, and the characteristic velocities are

=1× 10-7m,

=1× 10-5 and

=1× 10-3

CE

Table 2. Parameters of simulation cases (groups A-C)

AC

Cases A1 A2 A3 A4 A5 Velocity(mm/s) 1 1.5 1 1.5 1 Wettability(o) 140 140 140 140 140 Fracture orientation Null Null Null Null Null Null Null Null Null Null Fracture aperture( ) Null Null Null Null Null Fracture length( ) Surface tension coefficient(mN/m) 25 25 25 25 25 Viscosity of wetting phase (cP) 1.0 1.0 1.0 1.0 1.0 Viscosity of non-wetting phase(cP) 1.0 1.0 1.5 1.5 1.0 Viscosity ratio 1.0 1.0 1.5 1.5 1.0 Porosity(%) 22.4 22.4 24.7 24.7 34 (1H is the abbreviation of horizontal, 2V is the abbreviation of vertical)

17

A6 1.5 140 Null Null Null 25 1.0 1.0 1.0 34

B1 1 140 H1 2 170 25 1.0 1.0 1.0 24.7

B2 1.5 140 H 2 170 25 1.0 1.0 1.0 24.7

C1 1 140 V2 2 80 25 1.0 1.0 1.0 23.1

ACCEPTED MANUSCRIPT

CR IP T

Table 3. Parameters of simulation cases (groups C-G) Cases C2 D1 D2 E1 E2 F1 F2 Velocity(mm/s) 1.5 1 1.5 1 1.5 1 1.5 Wettability(o) 140 140 140 140 140 140 140 Fracture orientation V H H H H H H 2 4 4 1 1 1 1 Fracture aperture( ) 80 170 170 80 80 150 150 Fracture length( ) Surface tension coefficient(mN/m) 25 25 25 25 25 25 25 Viscosity of wetting phase (cP) 1.0 1.0 1.0 1.0 1.0 1.0 1.0 Viscosity of non-wetting phase (cP) 1.0 1.0 1.0 1.0 1.0 1.0 1.0 Viscosity ratio 1.0 1.0 1.0 1.0 1.0 1.0 1.0 Porosity(%) 23.1 27.2 27.2 23.1 29.9 29.9 23.2 These eighteen cases are classified into seven groups: A, B, C, D, E, F and G. A1, A2,

G1 1 140 H 2 80 25 1.0 1.0 1.0 23.2

A3, A4, A5, and A6 are in group A, with no fracture; B1 and B2 are in group B with one horizontal fracture with an aperture of 2

and a fracture length of 170

and a fracture length of 80

AN US

group C with one vertical fracture with an aperture of 2

and D2 are in group D with one horizontal fracture with an aperture of 4 length of 170

; D1

and a fracture

; E1 and E2 are in group E with one horizontal fracture with an aperture of

and a fracture length of 80

with an aperture of 1

; F1 and F2 are in group F with one horizontal fracture

M

1

; C1 and C2 are in

and a fracture length of 150

ED

horizontal fracture with an aperture of 2

; G1 and G2 are in group G with one

and a fracture length of 80

. A1, A2, A5 and and to

PT

A6 in group A are designated as the reference cases to calculate the DCP coefficient

investigate the effects of a fracture in groups B–F. Porous models in cases A5 and A6 are

CE

generated from the porous models in Figure 3-a. A3 and A4 in group A are designated to

AC

investigate the effect of viscosity ratio on capillary pressure. Groups B, C, and G are designated to investigate the effects of fracture orientation on DCP. Groups E and F are designated to investigate the fracture length on DCP. Groups B, D, E, and F are designated to investigate the effects of fracture aperture on DCP. 4.1. DCP of unfractured porous media

18

G2 1.5 140 H 2 80 25 1.0 1.0 1.0 23.1

ACCEPTED MANUSCRIPT Figure 4 shows the simulated phase evolution of A1 and A2 (Table 2) at different injection rates. There are three main features in these snapshots. First, convex non-wetting phase fronts form at the beginning (Figures 4-a1 and 4-b1). Next, at the dynamical state, the non-wetting phase does not have enough time to flow through all the flow paths (Figures 4-a2

AN US

CR IP T

and 4-b2), and the channeling of the non-wetting phase is evident.

ED

M

Figure 4. Snapshots of the phase distribution of simulation cases A1 and A2: (a1) the snapshot of A1 at 40ms; (a2) the snapshot of A1 at 60ms; (a3) the snapshot of A1 at 80ms; (b1) the snapshot of A2 at 30ms; (b2) the snapshot of A2 at 40ms; (b3) the snapshot of A2 at 50ms (The red color represents the displacing non-wetting phase; the light gray color represents the pore space filled with the wetting phase).

PT

Figure 4 also shows that the channeling pattern becomes more obvious when the injection flux is higher, since more pore space filled with the wetting phase is left behind for

CE

the high injection flux case. The flow pattern in Figure 4 can be partially interpreted by capillary number Ca=U

/( σcos(θ)) (Friedman, 1999) , where U is the inlet velocity,

AC

the wetting phase viscosity,

is

is the porosity, σ is the surface tension coefficient, θ is the

contact angle, and viscosity ratio of the non-wetting and wetting phases M= (Lenormand et al., 1988), where

is the viscosity of the non-wetting phase. Since

,

, , σ, M are the same in cases A1 and A2, so the flow pattern is determined by Ca, which is controlled by the injecting flux, and the value of Ca is larger for case A2 than for case A1.

19

ACCEPTED MANUSCRIPT Figure 5-a shows the evolution of DCP and saturation of A1, A2, A3 and A4, where DCP increases almost linearly before reaching a peak value and then decreases slowly to an asymptote. Similar results are obtained in the experiments of Bottero et al. (2011) and numerical simulations of Hao and Cheng (2010). Figure 5-b shows that the saturation of the

PT

ED

M

AN US

CR IP T

non-wetting phase first increases rapidly and then slows down to approach an asymptote.

AC

CE

Figure 5. Evolution of dynamic capillary pressure and saturation of non-wetting phase: (a) dynamic capillary pressure versus time of A1 and A2; (b)Saturation of nonwetting phase versus time of A1 and A2; (c) dynamic capillary pressure versus time of A3 and A4; (b)Saturation of nonwetting phase versus time of A3 and A4.

Figure 6-a shows DCP versus saturation of the wetting phase of A1, A2, A3 and A4.

There are three main features in this figure as well. Firstly, a higher injection flux will lead to a higher DCP. DCP of A2 is higher than that of A1, which agrees with Eq. (3). Secondly, the viscosity ratio has notable influence on DCP curves. DCP of A3 is higher than that of A1. Thirdly, DCP curves show non-monotonic trends at lower saturation of the wetting phase,

20

ACCEPTED MANUSCRIPT which is sometimes called "overshoot" in DCP experiments by some previous investigators (Hassanizadeh et al., 2005). As Eq. (3) indicates that SCP is needed to calculate the DCP coefficient , simulations with ultra-low flux rates are needed to get SCP. However, simulations with ultra-low flux rates will require much more computational resources, so a standard constant-rate mercury

CR IP T

injection experiment (CRMI) is carried out with Automated System for Pore Examination 730 with the same sample used in CT scanning. The injecting rate is set at an ultra-low state of 0.0001cc/min. The CRMI results are plotted in Figure 6-a. The curve is converted using a multiplier factor of

, where

and

are the

= 140o and

are

and

AN US

contact angle and interfacial tension used in simulation cases , respectively, and their values =25mN/m, while

and

are the corresponding values of

used in the experiment, and their values are

= 140o and

=485mN/m.

M

According to Eq. (3), to calculate the DCP coefficient , three terms are needed: DCP,

ED

SCP and the rate of change of saturation. To facilitate the calculation, DCP and SCP at specific saturation are firstly calculated by a linear interpolation between two adjacent data

PT

points. Then, values of the rate of change of saturation are calculated using a centraldifference scheme based on the data of Figure 5-b. Next, at a given saturation, the DCP

CE

difference

and

are given in pairs, and plotted in Figure 6-b for the

AC

range of 0.45
,

, which provides an additional data

) in Figure 6-b. All these data are used to fit a regression line of

in Figure 6-b, and the slope of the fitting line is the estimated value of .

21

ACCEPTED MANUSCRIPT In previous experimental works of Hassanizadeh et al. (2005) , Bottero et al. (2011) and Hou et al. (2012), all data under different injection fluxes are used to fit the DCP coefficient In this study, we calculate the DCP coefficient for each simulation case based on Figure 6-a, based on a regression procedure using the least-squares fitting method (Bottero et al., 2011) (see Table 4). The slope of the fitting line of cases A1&A2 in Figure 6-b is

CR IP T

.

Table 4 gives the DCP coefficients , coefficient of determination (R2) and 95% confidence

PT

ED

M

AN US

bounds of the fitting lines in Figure 6-b.

AC

CE

Figure 6. DCP and DCP coefficient: (a) dynamic capillary pressure and quasi-static capillary pressure;(b)Fitting curve of versus according Equation (2-3) (The saturation of wetting phase is calculated based on saturation of nonwetting phase according to ). Figure 6-b also shows that the viscosity ratio has notable influence on DCP coefficient

which is

for case A3&A4 with viscosity ratio of 1.5. This indicates that

increasing viscosity ratio may lead to increasing dynamic capillary pressure coefficient. Such a phenomenon is also observed in the work of Goel and O'Carroll (2011). We find that the ratio of between A1&A2 and A3&A4 is 1.61, which is close to ratio of the non-wetting phase viscosities between A1&A2 and A3&A4, which is 1.5 (see Table 2). In the work of 22

ACCEPTED MANUSCRIPT Niasar et al(2010) a similar conclusion is also drawn at the wetting phase saturation range of about 0.5
Cases

A2

A1&A2

A3

A4

0.94

0.89

0.90

0.97

3.74-4.11

4.01-4.23

3.96 -4.23

6.28-7.08

A3&A4 0.97

6.45-6.80

6.47-6.80

CR IP T

R-square 95% confidence bounds

A1 3.97 0.81

4.2. Effects of a horizontal or vertical fracture on two-phase flow

Figure 7 shows the results of evolution of phase spatial distribution with time for cases

AN US

A1, A5, B1, C1, D1, E1, F1 and G1. It is found that a fracture has a significant influence on two-phase flow in porous media, which is controlled by fracture orientation, aperture and length. The non-wetting fluid phase moves forward much faster in the fracture than in the adjacent porous matrix in B1 and D1, indicating the high conductivity of these horizontal

M

fractures. Despite the fact that the fracture orientations in B1 and D1 are the same (horizontal), the saturation distributions of the non-wetting fluid in these two cases are clearly

ED

different. Specifically, a greater fracture aperture leads to a smaller saturation of the non-

PT

wetting phase in the adjacent porous matrix. Quantitative analysis of such phenomenon is important for studying oil accumulation paths in fractured reservoirs or block faulted

CE

reservoirs (Lyu et al., 2017), since the crude oil in these systems needs to be extracted from water in the fractures and porous media during the accumulation process.

AC

The special parameter settings in our simulations indicate that when the fracture aperture

is larger than

, none of the non-wetting phases will flow into the adjacent porous

matrix owing to the much lower capillary pressure of fracture compared to those of the adjacent porous matrix. When the aperture of fracture is less than

, the non-

wetting fluid may partially saturate the adjacent porous matrix, since decreasing the fracture

23

ACCEPTED MANUSCRIPT aperture also decreases the difference in permeability and capillary pressure between the

M

AN US

CR IP T

fracture and the adjacent porous matrix (Wang, 2017, Zhang et al., 2016).

ED

Figure 7. Spatial distribution of fluid phases with an injection flux of 1 mm/s at 50ms, 60ms, 70ms, and 80ms (The red color represents the displacing non-wetting phase; the light gray color represents the pore space filled with the wetting phase).

PT

Figure 7 also indicates that fracture orientation is a crucial parameter in analyzing the flow paths in fractured porous media. Compared to a horizontal fracture (G1), a vertical

CE

fracture has less influence on flow during drainage (C1), provided that their apertures and

AC

lengths are the same as their counterparts in a horizontal fracture (see Tables 2 and 3). However, a non-equilibrium delayed effect on two-phase flow can be observed for a vertical fracture by comparing the simulation results of C1 and A1 at 60ms. The spatial distribution of the non-wetting phase shows no difference before reaching the vertical fracture, but the nonwetting phase displacement of the wetting phase becomes slower on the right of the vertical fracture in C1, as compared with A1 at 60ms. This is because the non-wetting phase will preferentially fill the vertical fracture first. Furthermore, Figure 7 shows that the fracture 24

ACCEPTED MANUSCRIPT length is another controlling factor of the flow pattern during the non-wetting phase displacement of the wetting phase in fractured porous media, as the longer fracture (B1) shows more influence on the flow pattern than the shorter fracture (G1) when the fracture

ED

M

AN US

CR IP T

apertures are the same.

CE

PT

Figure 8. Evolution of saturation of the non-wetting phase and DCP: (a) DCP versus time of cases A1, A5, B1, C1, D1, E1, F1 and G1; (b) DCP versus time of cases A2, A6, B2, C2, D2, E2, F2 and G2; (c) Saturation versus time of cases A1, A5, B1, C1, D1, E1, F1 and G1; (d) Saturation versus time of cases A2, A6, B2, C2, D2, E2, F2 and G2.

AC

Figure 8-a and Figure 8-b show the effects of microfractures on the DCP curves. It is observed that DCPs of the fractured porous media (B1, B2, C1, C2, D1, D2, E1, E2, F1, F2 and G1, G2) are lower than that of the original unfractured porous medium (A1 and A2) at the same injection flux, and a smaller fracture aperture will result in a higher DCP (A1, A2, B1, B2 and D1, D2) at the same injection flux; a shorter fracture length will lead to a higher DCP as well (A1, A2, B1, B2 and G1, G2; A1, A2, E1, E2 and F1, F2) at the same injection flux. Figure 8-c and Figure 8-d show the corresponding effects of microfractures on the 25

ACCEPTED MANUSCRIPT saturation of the non-wetting phase. The residual saturation values of the wetting phase in the fractured porous media (B1, B2, C1, C2, D1, D2, E1, E2, F1, F2, G1 and G2) are higher than those in the unfractured porous media (A1, A2) at the same injection flux. Figures 7-9 also show that porosity has crucial influences on DCP and saturation evolution. However, the relationship between porosity and DCP is very complex. This is

CR IP T

primarily due to two reasons. Firstly, DCP is not only controlled by fracture aperture and fracture length, but also controlled by fracture orientation (C1 and G1), which has no notable influence on porosity. Secondly, cases A1, A5, B1 and D1 show that the way the porosity is increased will also influence the effects of porosity on DCP. For instance, although the

AN US

increased porosity is higher in A5 comparing with D1, D1 actually has more influence on DCP than A5. The issue of porosity-DCP relationship deserves further investigation in the future.

M

Nevertheless, the change of DCP in D1 follows a trend similar to those in A1, B1, C1, , which is very

close the theoretical value of SCP in a fracture, which is

, where

ED

E1, F1 and G1. The DCP in D1 decreases to an asymptote of about

is

PT

the value of half of the fracture aperture (Zheng et al., 2017). This finding also agrees with the discussion on the spatial distribution of the non-wetting phase in Figure 7, which

CE

indicates that flow exchange between the fracture and the adjacent porous matrix ceases for a fracture with an aperture higher than

.

AC

Figure 8-a also shows that the vertical fracture has a transient delayed effect on DCP.

Comparing with A1, one can see that before reaching the peak value, DCP of C1 goes through a delayed period, marked as 2 in Figure 8-a at approximately 50ms. Hence, DCP of C1 first increases almost linearly (marked as 1 in Figure 8-a) and then slows down at a delayed period (marked as 2). After the delayed period, DCP of C1 increases again to the peak value (marked as 3) which is followed by a rapidly decreasing period (marked as 4) In 26

ACCEPTED MANUSCRIPT contrast, DCP of A1 increases almost linearly (marked as 1)to a peak value higher than that in B1 in the delayed period, and DCP of A1 decreases slowly after reaching the peak value (marked as 2 and 3). Comparing the snapshots of C1 in Figure 7 with DCP of C1 in Figure 8-a, it is found that

CR IP T

the non-wetting phase also invades the vertical fracture around the delayed period of 5060ms. C2 also goes through a similar delayed period as C1 at about 40ms (Figure 8-b) when the non-wetting phase invades the vertical fracture at 40ms (see Figure 9). Therefore, one

AC

CE

PT

ED

M

AN US

may conclude that the delayed effect in Figure 8-a is dominated by the vertical fracture.

Figure 9. Spatial distribution of fluid phases with an injection flux of 1.5mm/s at 34ms, 40 ms, 46 ms, and 52ms (The red color represents the displacing non-wetting phase; the light gray color represents the pore space filled with wetting phase).

Figure 9 shows that the evolution rate of the spatial distribution of the non-wetting phase with an injection flux of 1.5mm/s is much faster than that with an injection flux of 1mm/s. The DCP curves are also higher for the cases with higher saturation rate changes of the non27

ACCEPTED MANUSCRIPT wetting phase (see Figure 8-b). This agrees with Eq. (3), which is also validated in DCP experiments reported in a number of studies (Bottero et al., 2011). Previous numerical analysis and experiments mostly used constant-pressure conditions (Das and Mirzaei, 2012, Hou et al., 2012, Bottero et al., 2011), while the simulations of this study show that such a

CR IP T

phenomenon is also observed under constant-flux conditions. 4.3. Effects of a horizontal or vertical fracture on the DCP coefficient

Figure 10 shows that the fracture has a significant influence on the DCP versus saturation curves for the non-wetting phase. The absolute value of DCP not only varies greatly with

AN US

changes in the characteristics of fracture, but the maximum saturation is also dominated by the fracture. Figure 10 also shows a non-monotonic trend of DCP with increasing nonwetting phase saturation, which is different from the monotonic trend observed in static (Zhang et al., 2017b, Tang et al., 2017) or quasi-static capillary pressure experiments (Wang

M

and Alvarado, 2017). Such a non-monotonic trend or overshoot has also been observed in

AC

CE

PT

ED

high-pressure non-equilibrium displacement experiments before (Hassanizadeh et al., 2005).

Figure 10. DCP versus saturation of the non-wetting phase : (a) simulation cases with an injection flux of 1mm/s; (b) simulation cases with an injection flux of 1.5mm/s. The saturation of non-wetting phase and the saturation of wetting phase is bounded by a constraint of .

28

ACCEPTED MANUSCRIPT Although overshoot has been observed in experimental investigations before (Bottero et al., 2011, Hassanizadeh et al., 2005), the essential mechanism is still not fully understood to the best of our knowledge. From the simulation results, we find that with a higher injection rate, a portion of the non-wetting phase will enter nearby smaller pore spaces rather than the farther larger pore spaces at early stage, such as shown in the snapshots in Figure 7. This in

CR IP T

turn could cause DCP to be higher than SCP at early stage. As time processes, the whole system will eventually reestablish its equilibrium state, and the non-wetting phase will reenter the larger pore spaces. As a result, DCP declines with time. It should be noted that above explanation is based on the observation of the non-wetting phase in our simulation results,

AN US

not from a rigorous theoretical perspective.

According to Eq. (3), a straightforward way to calculate the DCP coefficient is to use one DCP curve, one SCP curve, and the change rate of saturation. The imbedded

M

microfracture will change the SCP curve. However, comparing with DCP, measuring or simulating SCP is often time consuming (Bottero et al., 2011, Hassanizadeh et al., 2005).

,

generated from the same porous model with different injection

PT

curves

ED

Here, we propose an alternative way to calculate the DCP coefficient based on two DCP

conditions listed in Tables 2 and 3, and two corresponding saturation rate change curves .

CE

and

,

,

and

AC

The procedure is briefly explained as follows. Firstly, curves

at the same saturation position are calculated by a linear interpolation between two

data points (Bottero et al., 2011). Then, one can obtain the following equations: (19) (20) By subtracting Eq. (19) from Eq. (20), one can eliminate the SCP term

29

, and gets:

ACCEPTED MANUSCRIPT

(21) According to Equation (21), the DCP coefficient could be calculated at specific saturation of the wetting phase. The resulting DCP coefficient

values are plotted versus the

wetting phase saturation in Figure 11. For very high viscosity ratios of , Goel and O'Carroll (2011), Goel et al. (2016) found that the DCP coefficient will

CR IP T

increase with the decrease of the wetting phase saturation. Our simulation results also show an increase trend of the DCP coefficient when the wetting phase saturation drops from 0.65 to 0.45 (case A1&A2). However, such an increase is very minor as compared with the

AN US

experimental results of Goel and O'Carroll (2011), Goel et al. (2016), but it is similar to the experimental results of Bottero et al. (2011). It must be noted that in our study the DCP coefficient is not calculated over the entire range of saturation and the viscosity ratio is smaller (1 and 1.5) than that used in the experiments of Goel and O'Carroll (2011), Goel et al.

M

(2016). This probably explains the discrepancy of the simulation results of this study and the

ED

experimental results of Goel and O'Carroll (2011), Goel et al. (2016). However, further

AC

CE

PT

research is necessary to decipher such a discrepancy.

Figure 11.The DCP coefficient versus local saturation of wetting phase Sw for fractured and unfractured porous media: (a) versus Sw of A1&A2, A5&A6, C1&C2, E1&E2, F1&F2, G1&G2 ; (b) versus Sw of B1&B2, D1&D2. Since the DCP coefficient varies a large scale, semi logarithmic plot is used in (a) to avoid compressing the curves.

30

ACCEPTED MANUSCRIPT

In this study, an alternative procedure inspired by Bottero et al. (2011) and Hou et al. (2012) is used to calculate the DCP coefficient . In this procedure, the plots of difference of DCP versus difference of change rate of saturation in Figure 12, are used to determine the values of DCP coefficient (see Equation (21)). According to the work of Bottero et al. , then

to

, leading

CR IP T

(2011), when and

This provides one additional

fitting data point. Fitting is optimized by applying a least-squares algorithm, with the difference in the saturation rate change term

as the dependent variable. The fitting

AN US

and the difference in DCP

as the independent variable,

lines of the fractured and unfractured porous models are shown in Figure 11. Table 5 lists the corresponding and coefficient of determination (R2) of the fitting lines in Figure 11. Table 5

M

indicates that fracture, regardless horizontal or vertical, has a significant effect on the DCP coefficient of fractured porous media, and the values of fractured porous media are one to

ED

two orders of magnitude lower than those of unfractured porous media. For unfractured porous media (A1 and A2), is

, which is close to the value calculated

PT

using another method in section 4.1.

CE

It should be noted that the previously reported experimental DCP coefficient varies over a few orders of magnitude, despite the similarities in the porous media and fluids used

AC

(Hassanizadeh et al., 2002). Hou et al. (2012) argued that the changing rate of the wetting phase saturation in those experiments should be corrected before calculating the DCP coefficient

After such a correction, the experimental DCP coefficient was about

at the wetting phase saturation of

in Hou et al. (2012), which is slightly

smaller than the DCP coefficient of cases A1&A2 of this study (see Table 5).

31

ACCEPTED MANUSCRIPT Figure 12-a indicates that both fractures (D1 and D2) and eroding process (A5 and A6) have significant influence on the DCP coefficient . Figure 11-b shows that fracture orientation has a strong influence on . The DCP coefficient of the porous model with a vertical fracture (C1 and C2) is

; while the DCP coefficient

fracture length of 80

, with the same

.

CR IP T

model with a horizontal fracture (G1 and G2) reduces to

of the porous

It is found that a greater fracture aperture leads to a smaller value of . The DCP coefficient of porous model with a fracture aperture of 4

(D1 and D2) is

(B1

AN US

which is 22 times smaller than that of the porous model with a fracture aperture of 2

,

and B2) , and it is about 40 times smaller than that of the porous model with a fracture aperture of 1

(E1 and E2, F1 and F2).

The DCP coefficient for a porous model with a horizontal fracture is one order of

M

magnitude smaller than that of the porous model with a vertical fracture (C1 and C2), and two orders of magnitude smaller than that of the unfractured porous model (A1 and A2). , doubling the fracture length will cause the reduction of

ED

When the fracture aperture is 2

0

(G1 and G2) to

PT

the value by 21.29%, from

B2). When the fracture aperture decreases by half to 1

, doubling the fracture length causes (F1 and F2) to

CE

the reduction of the value by 15.78%, from

(B1 and

AC

(E1 and E2). Above analysis reflects the complex interrelated effects of fracture

length and fracture aperture on the DCP coefficient. The effects of DCP and DCP coefficient of fracture on porous media are important for

many practical applications, such as petroleum development (Tang et al., 2017, Yu et al., 2015) and storage of CO2 in stimulated saline aquifers (Raziperchikolaee et al., 2013). In petroleum development, CO2 huff-n-puff and hydraulic fracturing are effective methods to increase the oil recovery (Tang et al., 2017). One key aspect to simulate and predict the 32

ACCEPTED MANUSCRIPT multi-phase flow in CO2 huff-n-puff is analyzing the properties of hydraulic fractures. From a field-scale reservoir simulation perspective, the actual hydraulic fracture is usually treated as a wide „pseudo‟ fracture to increase the numeric stability (Tang et al., 2017, Yu et al., 2015). The properties of the „pseudo‟ fracture will affect the oil-gas distribution and oil recovery, and one key property of the „pseudo‟ fracture is DCP. The simulation results from our studies

CR IP T

could be used to understand the performance of hydraulic fractures. It should be noted that the actual fracture geometry in oil field is more complex than those used in our model. So further research is still needed to expand this research to more complex fracture networks in order to better understand the performance of hydraulic fractures in reservoir scale.

AN US

Since the storage of CO2 in stimulated saline aquifers also involves the hydraulic

fracturing technique, thus the „pseudo‟ fracture concept is also applicable. Therefore, the effects of fracture on DCP and DCP coefficient will also affect the simulated CO2 plume,

M

whose shape will affect the long term trapping in stimulated saline aquifers (Raziperchikolaee

AC

CE

PT

ED

et al., 2013).

Figure 12. The DCP coefficient of fractured and unfractured porous media: (a) effects of fracture and eroding process; (b) effects of fracture length, aperture and orientation. Difference of is , difference of the change rate of is .

33

ACCEPTED MANUSCRIPT Table 5. The fitting results of the DCP coefficient A1&A2 A5&A6 B1&B2 C1&C2 D1&D2 E1&E2 F1&F2 G1&G2 39.8 3.40 2.44 5.68 0.11 4.11 0 2 R 0.61 0.87 0.95 0.96 0.83 0.96 0.98 0.97 (R2 is calculated based on CFTOOL in MATLAB ) Cases

5. Conclusions A series of two-phase displacement flow simulations are carried out to investigate the

CR IP T

effect of a microfracture on DCP and the DCP coefficient . The two-phase LBM and CT image-based three-dimensional porous modeling method are used to execute the simulations. Based on the simulated results, DCP and the DCP coefficient

are calculated for two

AN US

unfractured porous models and six fractured models with different fracture orientations,

apertures and lengths. In this study the focus has been to investigate the influence of fracture on the DCP and DCP coefficient . In addition, we also investigate the effects of viscosity ratio on the DCP and DCP coefficient . We have proposed an innovative way to calculate

M

the DCP coefficient ε based on two DCP curves from the same porous model with different

ED

injection conditions. This method yields the value of τ which agrees with that obtained using the conventional method based on the SCP and DCP curves. From the simulation results,

PT

non-monotonic trends in DCP versus saturation curves are observed. The DCP curves are higher for the cases with higher saturation rate changes of the non-wetting phase. The

CE

relationship between porosity and DCP is very complex and the way the porosity is increased will also influence the effects of porosity on DCP. The results also show that increasing

AC

viscosity ratio may lead to increasing DCP coefficient τ. This study shows that fracture have significant influence on two-phase flow in porous media, which is controlled by fracture orientation, aperture and length. With the same fracture aperture value, a horizontal fracture has a greater influence on DCP and the DCP coefficient than a vertical fracture. With the same fracture orientation, a smaller fracture aperture of the horizontal fracture leads to a higher DCP, and a longer length of the horizontal fracture will result in a smaller DCP. 34

ACCEPTED MANUSCRIPT A microfracture has

significant effects on the DCP coefficient . For unfractured porous

media, the DCP coefficient is about

. For fractured porous media, it is

for the porous model with a horizontal fracture with an aperture of 2

4

for the porous model with a horizontal fracture with an aperture of , and it is

aperture of 2

for the porous model with a vertical fracture with an . In addition,the horizontal fracture length and aperture have complex

CR IP T

it is

interrelated effects on DCP and the DCP coefficient . The vertical fracture has a non-

AN US

equilibrium delayed effect on DCP and the DCP coefficient .

Acknowledgments: This study was supported by the Natural Science Foundation of China (NO: 41330313, 41402108, 41406050) , the Foundation of Shandong Province (NO: ZR2014DL003) and Talent Introduction Funding (NO: 2013059), China University of

M

Petroleum (East China) Financing Innovation Funding (NO: 2014010336 and NO:

ED

27R1510046A). We would like to thank the China Scholarship Council for their financial support. We are also thankful to the Editors and anonymous reviewers for their constructive

PT

and insightful comments and suggestions.

CE

Author Contributions: Mingming Tang, Hongbin Zhan and Shuangfang Lu conceived the study; Mingming Tang conducted the research work under the supervision of Hongbin Zhan

AC

and Shuangfang Lu. All co-authors revised and commented the final draft of the paper. Conflicts of Interest: None. References ABIDOYE, L. K. & DAS, D. B. 2014. Scale dependent dynamic capillary pressure effect for two-phase flow in porous media. Advances in Water Resources, 74, 212-230.

35

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

BADALASSI, V. E., CENICEROS, H. D. & BANERJEE, S. 2003. Computation of multiphase systems with phase field models. Journal of Computational Physics, 190, 371-397. BAO, J. & SCHAEFER, L. 2013. Lattice Boltzmann equation model for multi-component multi-phase flow with high density ratios. Applied Mathematical Modelling, 37, 18601871. BERG, S. & OTT, H. 2012. Stability of CO2-brine immiscible displacement. International Journal of Greenhouse Gas Control, 11, 188-203. BOTTERO, S., HASSANIZADEH, S. M., KLEINGELD, P. J. & HEIMOVAARA, T. J. 2011. Nonequilibrium capillarity effects in two-phase flow through porous media at different scales. Water Resources Research, 47, W10505. CAHN, J. W. & HILLIARD, J. E. 1958. Free energy of a nonuniform system. I. Interfacial free energy. The Journal of Chemical Physics, 28, 258-267. CAO, P., LIU, J. S. & LEONG, Y. K. 2017. A multiscale-multiphase simulation model for the evaluation of shale gas recovery coupled the effect of water flowback. Fuel, 199, 191-205. CHEN, S., CHEN, H., MARTNEZ, D. & MATTHAEUS, W. 1991. Lattice Boltzmann model for simulation of magnetohydrodynamics. Physical Review Letters, 67, 3776-3779. CHEN, S. & DOOLEN, G. D. 1998. Lattice Boltzmann method for fluid flows. Annual Review of Fluid Mechanics, 30, 329-364. DAHLE, H. K., CELIA, M. A. & HASSANIZADEH, S. M. 2005. Bundle-of-tubes model for calculating dynamic effects in the capillary-pressure- saturation relationship. Transport in Porous Media, 58, 5-22. DAS, D. B., HANSPAL, N. S. & NASSEHI, V. 2005. Analysis of hydrodynamic conditions in adjacent free and heterogeneous porous flow domains. Hydrological Processes, 19, 2775-2799. DAS, D. B. & MIRZAEI, M. 2012. Dynamic effects in capillary pressure relationships for two-phase flow in porous media: Experiments and numerical analyses. AIChE Journal, 58, 3891-3903. ESHRAGHI, M., FELICELLI, S. D. & JELINEK, B. 2012. Three dimensional simulation of solutal dendrite growth using lattice Boltzmann and cellular automaton methods. Journal of Crystal Growth, 354, 129-134. FAYBISHENKO, B. 2017. Detecting dynamic causal inference in nonlinear two-phase fracture flow. Advances in Water Resources, 106, 111-120. FRIEDMAN, S. 1999. Dynamic contact angle explanation of flow rate-dependent saturation pressure relationships during transient liquid flows in unsaturated porous media. J Adhesion Sci Technol 13, 1495-1581. FRISH, U., D'HUMIERES, D., HASSLACHER, B., LALLEMAND, P., POMEAU, Y. & RIVET, J. P. 1987. Lattice gas hydrodynamics in two and three dimensions. Complex System, 1, 649-707.

36

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

GOEL, G., ABIDOYE, L. K., CHAHAR, B. R. & DAS, D. B. 2016. Scale dependency of dynamic relative permeability–satuartion curves in relation with fluid viscosity and dynamic capillary pressure effect. Environmental Fluid Mechanics, 16, 945-963. GOEL, G. & O'CARROLL, D. M. 2011. Experimental investigation of nonequilibrium capillarity effects: Fluid viscosity effects. Water Resources Research, 47, 15. GUNSTENSEN, A. K., ROTHMAN, D. H., ZALESKI, S. & ZANETTI, G. 1991. Lattice Boltzmann model of immiscible fluids. Physical Review A, 43, ID4320. HALL, M. R., RIGBY, S. P., DIM, P., BATEMAN, K., MACKINTOSH, S. J. & ROCHELLE, C. A. 2016. Post-CO2 injection alteration of the pore network and intrinsic permeability tensor for a Permo-Triassic sandstone. Geofluids, 16, 249-263. HAO, L. & CHENG, P. 2010. Lattice Boltzmann simulations of water transport in gas diffusion layer of a polymer electrolyte membrane fuel cell. Journal of Power Sources, 195, 3870-3881. HASSANIZADEH, S. M., CELIA, M. A. & DAHLE, H. K. 2002. Dynamic Effect in the Capillary Pressure-Saturation Relationship and its Impacts on Unsaturated Flow. Vadose Zone Journal, 1, 38-57. HASSANIZADEH, S. M., OUNG, O. & MANTHEY, S. 2005. Laboratory experiments and simulations on the significance of non-equilibrium effect in the capillary pressuresaturation relationship. In: SCHANZ, T. (ed.) Unsaturated Soils: Experimental Studies: Proceedings of the International Conference “From Experimental Evidence towards Numerical Modeling of Unsaturated Soils,” Weimar, Germany, September 18–19, 2003 Volume I. Berlin, Heidelberg: Springer Berlin Heidelberg. HAZI, G. 2003. Accuracy of the lattice Boltzmann method based on analytical solutions. Pysical Review E, 67, ID-056705. HE , W. S., YI , J. S. & VAN , N. T. 2000. Two-phase flow model of the cathode of PEM fuel cells using interdigitated flow fields. AIChE Journal, 46, 2053-2064. HIRT, C. W. & NICHOLS, B. D. 1981. Volume of fluid (VOF) method for the dynamics of free boundaries. Journal of Computational Physics, 39, 201-225. HOU, L., CHEN, L. & KIBBEY, T. C. G. 2012. Dynamic capillary effects in a small-volume unsaturated porous medium: Implications of sensor response and gas pressure gradients for understanding system dependencies. Water Resources Research, 48, W11522. HUANG, H., THORNE, D. T., JR., SCHAAP, M. G. & SUKOP, M. C. 2007. Proposed approximation for contact angles in Shan-and-Chen-type multicomponent multiphase lattice Boltzmann models. Physical Review E, 76, ID 066701. JERBI, C., FOURNO, A., NOETINGER, B. & DELAY, F. 2017. A new estimation of equivalent matrix block sizes in fractured media with two-phase flow applications in dual porosity models. Journal of Hydrology, 548, 508-523. JUANES, R., SPITERI, E. J., ORR, F. M., JR. & BLUNT, M. J. 2006. Impact of relative permeability hysteresis on geological CO2 storage. Water Resources Research, 42, W12418.

37

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

KANG, S. K. & HASSAN, Y. A. 2013. The effect of lattice models within the lattice Boltzmann method in the simulation of wall-bounded turbulent flows. Journal of Computational Physics, 232, 100-117. KIM, S. & SANTAMARINA, J. C. 2016. Geometry-coupled reactive fluid transport at the fracture scale: application to CO2 geologic storage. Geofluids, 16, 329-341. KUWATA, Y. & SUGA, K. 2015. Anomaly of the lattice Boltzmann methods in threedimensional cylindrical flows. Journal of Computational Physics, 280, 563-569. LANDRY, C. J., KARPYN, Z. T. & AYALA, O. 2014. Relative permeability of homogenous-wet and mixed-wet porous media as determined by pore-scale lattice Boltzmann modeling. Water Resources Research, 50, 3672-3689. LANDRY, C. J., KARPYN, Z. T. & PIRI, M. 2011. Pore-scale analysis of trapped immiscible fluid structures and fluid interfacial areas in oil-wet and water-wet bead packs. Geofluids, 11, 209-227. LENORMAND, R., TOUBOUL, E. & ZARCONE, C. 1988. Numerical models and experiments on immiscible displacements in porous media. Journal of Fluid Mechanics, 189, 165-187. LI, D. L., XU, C. F., WANG, Y. X., SONG, Z. F., XIONG, M., GAO, X. & DENG, X. G. 2016. Parallelizing and optimizing large-scale 3D multi-phase flow simulations on the Tianhe-2 supercomputer. Concurrency and Computation-Practice & Experience, 28, 1678-1692. LICHTNER, P. C. 2000. Critique of dual continuum formulations of multicomponent reactive transport in fractured porous media. Geophys. Monogr. Ser., 122, 281-298. LIU, H., KANG, Q., LEONARDI, C. R., SCHMIESCHEK, S., NARV EZ, A., JONES, B. D., WILLIAMS, J. R., VALOCCHI, A. J. & HARTING, J. 2016. Multiphase lattice Boltzmann simulations for porous media applications. Computational Geosciences, 20, 777-805. LIU, H., VALOCCHI, A. J. & KANG, Q. 2012. Three-dimensional lattice Boltzmann model for immiscible two-phase flow simulations. Physical Review E, 85. LYU, W. Y., ZENG, L. B., ZHANG, B. J., MIAO, F. B., LYU, P. & DONG, S. Q. 2017. Influence of natural fractures on gas accumulation in the Upper Triassic tight gas sandstones in the northwestern Sichuan Basin, China. Marine and Petroleum Geology, 83, 60-72. MIRI, R., SHADIZADEH, S. R. & KHARRAT, R. 2014. Fracture capillary pressure based on the liquid bridge dynamic stability study. Energy Sources, Part A: Recovery, Utilization, and Environmental Effects, 36, 2536-2545. MOHAMAD, A. A. 2011. Lattice Boltzmann method-fundamentals and engineering applications with computer codes, Springer-Verlag London. NIASAR, V. J., HASSANIZADEH, S. M. & DAHLE, H. K. 2010. Non-equilibrium effects in capillarity and interfacial area in two-phase flow: dynamic pore-network modelling. Journal of Fluid Mechanics, 665, 38-71.

38

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

NOURGALIEV, R. R., DINH, T. N., THEOFANOUS, T. G. & JOSEPH, D. 2003. The lattice Boltzmann equation method: Theoretical interpretation, numerics and implications. Int. Multiphase Flow, J, 29, 117-169. OSHER, S. & SETHIAN, J. A. 1988. Fronts propagating with curvature-dependent speed algorithms based on Hamilton-Jacobi formulations. Journal of Computational Physics, 79, 12-49. PILLIOD, J. E. & PUCKETT, E. G. 2004. Second-order accurate volume-of-fluid algorithms for tracking material interfaces. Journal of Computational Physics, 199, 465-502. PRUESS, K. & GARCIA, J. 2002. Multiphase flow dynamics during CO2 disposal into saline aquifers. Environmental Geology, 42, 282-295. QI, R., LAFORCE, T. C. & BLUNT, M. J. 2009. Design of carbon dioxide storage in aquifers. International Journal of Greenhouse Gas Control, 3, 195-205. RAEINI, A. Q., BIJELJIC, B. & BLUNT, M. J. 2014a. Numerical modelling of sub-pore scale events in two-phase flow through porous media. Transport in Porous Media, 101, 191-213. RAEINI, A. Q., BLUNT, M. J. & BIJELJIC, B. 2014b. Direct simulations of two-phase flow on micro-CT images of porous media and upscaling of pore-scale forces. Advances in Water Resources, 74, 116-126. RAZIPERCHIKOLAEE, S., ALVARADO, V. & YIN, S. 2013. Effect of hydraulic fracturing on long-term storage of CO2 in stimulated saline aquifers. Applied Energy, 102, 1091-1104. REIDER, M. B. & STERLING, J. D. 1995. Accuracy of discrete-velocity BGK models for the simulation of the incompressible Navier-Stokes equations. Computers & fluids, 24, 459-467. SHAN, X. & CHEN, H. 1993. Lattice Boltzmann model for simulating flows with multiple phases and components. Physical Review E, 47, 1815-1819. SUCCI, S., MORADI, N., GREINER, A. & MELCHIONNA, S. 2014. Lattice Boltzmann modeling of water-like fluids. Front. Physics, 2, 1-13. SUGA, K., KUWATA, Y., TAKASHIMA, K. & CHIKASUE, R. 2015. A D3Q27 multiplerelaxation-time lattice Boltzmann method for turbulent flows. Computers & Mathematics with Applications, 69, 518-529. SUKOP, M. C. & THORNE, D. T. 2006. Lattice Boltzmann modeling-an introduction for geoscientists and engineers, Springer-Verlag Berlin Heidelberg. SUSSMAN, M., FATEMI, E., SMEREKA, P. & OSHER, S. 1998. An improved level set method for incompressible two-phase flows. Computers & Fluids, 27, 663-680. SWIFT, M. R., OSBORN, W. R. & YEOMANS, J. M. 1995. Lattice Boltzmann simulation of nonideal fluids. Physical Review Letters, 75, 830-833. TANG, M., LU, S., LIANG, H., YAN, B., MA, H. & SHEN, S. 2015. Numerical modelling and simulation of seepage in tight sandstone based on parallel lattice Boltzmann method. Acta Geologica Sinica - English Edition, 89, 412-414. TANG, M., ZHAO, H., MA, H., LU, S. & CHEN, Y. 2017. Study on CO2 huff-n-puff of horizontal wells in continental tight oil reservoirs. Fuel, 188, 140-154. 39

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

WANG, X. & ALVARADO, V. 2017. Effects of low-salinity waterflooding on capillary pressure hysteresis. Fuel, 207, 336-343. WANG, Y. F. 2017. On subsurface fracture opening and closure. Journal of Petroleum Science and Engineering, 155, 46-53. WU, Y. S. 1999. On the effective continuum method for modeling multiphase flow, multicomponent transport and heat transfer in fracture rock, Washington, D. C., American Geophysical Union. YU, W., LASHGARI, H. R., WU, K. & SEPEHRNOORI, K. 2015. CO2 injection for enhanced oil recovery in Bakken tight oil reservoirs. Fuel, 159, 354-363. ZHANG, C., CHENG, Y. & ZHANG, C. M. 2017a. An improved method for predicting permeability by combining electrical measurements and mercury injection capillary pressure data. Journal of Geophysics and Engineering, 14, 132-142. ZHANG, Q., DONG, Y. H., LIU, S. M., ELSWORTH, D. & ZHAO, Y. X. 2017b. Shale pore characterization using NMR cryoporometry with octamethylcyclotetrasiloxane as the probe liquid. Energy & Fuels, 31, 6951-6959. ZHANG, Q. F., HUANG, Z. Q., YAO, J., WANG, Y. Y. & LI, Y. 2017c. A multiscale mixed finite element method with oversampling for modeling flow in fractured reservoirs using discrete fracture model. Journal of Computational and Applied Mathematics, 323, 95-110. ZHANG, Z. B., LI, X. & HE, J. M. 2016. Numerical study on the permeability of the hydraulic-stimulated fracture network in naturally-fractured shale gas reservoirs. Water, 8, 14. ZHENG, J., LIU, H., WANG, K. & YOU, Z. 2017. A new capillary pressure model for fractal porous media using percolation theory. Journal of Natural Gas Science and Engineering, 41, 7-16.

40