Determination of inorganic and organic permeabilities of shale

International Journal of Coal Geology 215 (2019) 103296

Contents lists available at ScienceDirect

International Journal of Coal Geology journal homepage: www.elsevier.com/locate/coal

Determination of inorganic and organic permeabilities of shale a

b

a,⁎

c

Sheng Li , Qian Sang , Mingzhe Dong , Peng Luo a b c

T

Department of Chemical and Petroleum Engineering, University of Calgary, Calgary, AB T2N 1N4, Canada College of Petroleum Engineering, China University of Petroleum (Huadong), Qingdao, Shandong, China Mining and Energy Division, Saskatchewan Research Council, Regina, SK S4S 7J7, Canada

ARTICLE INFO

ABSTRACT

Keywords: Shale oil Vacuum-imbibition test Inorganic permeability Organic permeability Diffusion Bundle-of-tube model

Shale oil resources have become important energy supplies worldwide, but their storage features and flow mechanisms remain unclear. It is now well-accepted that shale formations are complex systems with both inorganic and organic contents. This study aims to characterize the permeabilities for inorganic and organic pore systems separately in shale samples using an innovative modeling approach validated with experimental data. Oil and water vacuum–imbibition tests were first conducted along with helium porosity measurements for four shale rock samples. A mathematical model considering imbibition and diffusion was then proposed and used to match the lab-measured oil and water imbibition curves. The critical parameters for these shale samples, including inorganic/organic porosities, fluid saturations, and permeabilities were successfully determined using this approach. Among all the samples in the study, the results showed that 33.93 to 40.54% of total oil in place was stored inside kerogen as either the free oil within the organic pores or the dissolved oil. Moreover, the flow of this part of the oil was controlled by the organic permeability as well as by the diffusion coefficient, which was different than the flow character of oil in the inorganic pore system. The organic permeabilities of the samples used in this study were between 2.24 × 10−6 and 1.59 × 10−5 md, which were 243 to 2741 times less than each sample's corresponding inorganic permeability. The sensitivity analysis indicated that organic permeability significantly affected the oil imbibition rate. The proposed methodology should be adopted in future reservoir characterizations, and the determined permeabilities for both inorganic and organic pore material should accordingly be considered in future reservoir simulations in order to accurately describe fluid flow in shale reservoirs.

1. Introduction Hydrocarbon production from organic-rich shale has gained worldwide interest in the industry as shale reservoirs are proving to be globally abundant and showing the potential to be commercially exploited (McGlade et al., 2013). Unlike other conventional or tight sandstone, and carbonate reservoirs, shale reservoirs are characterized by their rich organic matter and clay content. After lithification, the organic matter is converted to kerogen from which fossil fuel is derived during geological time. Shale can be both source rock and reservoir rock for hydrocarbons. The exploitation of shale oil reservoirs is still at the preliminary stage due to the limited understanding of the storage features and flow mechanisms in the shale rock. Considering there are both inorganic and organic media in shale, the conventional petro–physical parameters which ignore the variations in organic pore systems are not capable of fully describing a shale reservoir. Some researchers realized the

⁎

importance of characterizing organic material using radiation technologies, which primarily consist of scanning electron microscopy (SEM), nuclear magnetic resonance (NMR), and elemental dispersive spectroscopy (EDS) measurements (Barnett et al., 2012; Wang et al., 2018; Yuan et al., 2014). However, some nanopores are difficult to detect using SEM, and interpretation of NMR spectra requires parameters which are also difficult to estimate, such as total pore surface area (Curtis et al., 2011; Sigal, 2015). In addition to radiation methods, fluid-invasion methods are also used for determining the porosity and fluid saturation. For example, a mercury injection capillary pressure (MICP) test is commonly used for tight rock samples (Clarkson et al., 2013). However, an MICP test fails to touch nanopores in unconventional tight rocks due to extremely high pressure which could fracture the rock (Labani et al., 2013). Another drawback of MICP is that it cannot distinguish between the organic and inorganic pore networks because mercury is a non-wetting phase for all pores (Shi et al., 2019). Spontaneous imbibition is another fluid-invasion technique for

Corresponding author. E-mail address: [email protected] (M. Dong).

https://doi.org/10.1016/j.coal.2019.103296 Received 18 May 2019; Received in revised form 16 September 2019; Accepted 17 September 2019 Available online 17 October 2019 0166-5162/ © 2019 Elsevier B.V. All rights reserved.

International Journal of Coal Geology 215 (2019) 103296

S. Li, et al.

characterizing shale rocks. Lan et al. (2015) concluded the strong oilwetting affinity of the shale samples to the presence of the water-repellant pores with or coated by solid bitumen and pyrobitumen. Therefore, the equilibrated imbibed volume of oil was significantly greater than that of brine. However, the research did not recognize that part of the oil was not stored in the organic pores; rather, it dissolved into the kerogen instead, as both kerogen and molecular oil were nonpolar. The dissolved hydrocarbon has been studied primarily for shale gas reservoirs because the adsorption and diffusion of gas into kerogen happens much faster than that of oil. Lopez and Aguilera (2018) employed a petro-physical model to determine the inorganic porosity, adsorbed porosity, organic porosity, as well as the fractional volumes of solid kerogen available for gas diffusion. Javadpour et al. (2007) stated that the gas existing in shale gas reservoirs assumes three distinct forms: compressed gas in pores, adsorbed gas on the surface of pore walls, and dissolved gas in kerogen. Yang et al. (2016) investigated the dissolved gas stored in shales and stated that the hydrocarbon dissolution into kerogen is a diffusion-controlled, time-dependent process and that equilibrium was not achieved instantaneously. Sang et al. (2018) determined dissolved oil by conducting the vacuum–imbibition tests for a series of shale rock samples using an improved approach, by also measuring the helium porosity, so that the inorganic pore volume, organic pore volume, and dissolution pore volume could be distinguished. However, these researchers only considered the equilibrium status in a fluid-invasion process when characterized the pore systems. As a result, the effects of the pore systems on the dynamic process in shale are ignored and the inorganic and organic permeabilities in shale remain indistinguishable. From a modeling perspective, generating a multi-system pore network model is desirable in order to fully reflect the actual flow path naturally occurring in shale. However, it is not practical based on current computing power. Alternatively, abstract models were generated to mimic flow, among which, a tube-bundle model with a set of parallel-placed capillary tubes is one of the most important ones that represent the pore space in porous media (Bartley and Ruth, 1999). Fluid imbibition behavior in porous media has been extensively investigated using the tube-bundle model. Cai et al. (2014) modified the Hagen-Poiseuille and Laplace-Young equations for modeling tortuous, irregularly-shaped capillaries. Some researchers examined the pressure equilibrium and fluid transfer between tubes and their effects on the imbibition profile (Dong et al., 1998; Dullien, 2002; Li et al., 2017). However, none of these researches can be directly applied to shale rocks because local wettability heterogeneity exists in shale. Lan et al. (2015) and Yassin et al. (2016) made some preliminary attempts to applying the tube-bundle model to shale rocks by separating inorganic and organic tubes, in their studies, imbibition profiles were matched to obtain the wettability index and relative permeability, but the dissolved oil in kerogen was ignored in their models. In a recent study from Shi et al. (2019), the imbibition profile was incorporated with fractal model to analyze the pore size distribution. However, the author considered all the imbibed oil as free phase and dissolved oil was not considered either. If considering the oil dissolution process, the movement of oil in organic tubes has two directions: axial imbibition and radial diffusion. To the authors' best knowledge, this coupled imbibition and diffusion process has not been modeled and calculated; therefore, the organic permeability of shale has never been accurately characterized. Based on numerous experimental and simulation studies aiming to characterize shale reservoirs, it is now well-accepted that shale rock is a very complex system with local heterogeneous wettability. Because wettability affects the relative permeability and the capillary pressure, the complex fluid saturation and distribution of formation fluids have a significant impact on fluid flow in shale reservoirs. In order to understand the fluid storage features and transport mechanisms of formation fluids in shale rocks, the objectives of this study are to determine the permeability of the inorganic and organic pore systems, and to develop

a more accurate model for fluid flow in shale rocks. In this paper, we present the experimental procedures and results in Section 2. Section 3 introduces the detailed theory of the model for calculating the inorganic and organic permeabilities. Section 4 discusses the results of the model, petro–physical properties of the shale rock samples (i.e., inorganic/organic porosities, fluid saturations, and permeabilities), and the sensitivities of the organic permeability, total organic content (TOC), and diffusion coefficient, which are followed by conclusions in Section 5. 2. Experimental procedures and results The shale rock samples used in this study were collected from wells in the Shahejie Formation, Dongying Depression, located in Dongying, Shandong Province in eastern China. This formation was developed on the North China Platform basement in the Mesozoic–Cenozoic era. The lithology of the samples is identified as kerogen-rich limestone, argillaceous limestone, and sandstone. The mineral compositions of the samples include calcite, clay, and quartz, with small amounts of feldspar, dolomite, and pyrite. A total of 20 shale rock samples were tested in Sang et al. (2018). Among these samples, four representative ones which cover different wells, depths, and kerogen types were selected. Table 1 below summarizes the properties of the samples used in this study. In order to calculate the apparent density of the samples, the volume and weight of the core samples were measured and then the samples were crushed to small pieces with sizes of around 0.02 m. The volume of the crushed sample can be calculated from the weight and the apparent density. Helium was injected to the holder before and after the crushed, cleaned, and dried samples were put in, so that the volume of the holder, as well as the volume of free space (void volume + sample pore volume) could be measured, respectively. Then the void volume between samples and holder can be calculated from the volume of the holder and crushed sample bulk volume. After deducting the void volume between samples and holder, the pore volumes of samples were obtained and then converted to porosities which were defined as the ratio of pore volume to the bulk sample volume. The calculated porosities for each sample are also in the last row of Table 1. The pore structure of the shale sample was determined using SEM. In the previous study, the scanned images for samples in the same formation clearly indicated pores which were enclosed by organic material. Combined with the pore size distribution measurement, the sizes of organic pores within the kerogen ranged from 30 nm to 260 nm (Sang et al., 2018). The results were similar to Yassin‘s SEM/EDS analysis results of the samples from the Montney formation in Canada (Yassin et al., 2016). After measuring the porosity and characterizing the pore size distribution, the vacuumed sample was subjected to spontaneous imbibition tests using synthetic brine and n-dodecane, separately. The schematic diagram of the apparatus is shown in Fig. 1 below. The crushed samples were stored in two sealed test cells, and a high precision burette with graduation of 1 × 10−8 m3 (0.01 cm3) was connected to each test cell to measure the imbibed liquid volume. The top of the burettes Table 1 Geochemical characteristics, mass, and porosity of samples (Sang et al., 2018).

2

Number

#1

#2

#3

#4

Well name Depth (m) Mass (g) TOC (wt, %) Kerogen type Vitrinite reflectance (% Ro) Bulk volume (cm3) Porosity (%)

L67 3307 67.64 2.84 II 0.79 27.61 7.11

L67 3276 63.02 2.84 II 0.75 26.15 7.96

LY1 3624 55.03 2.39 III 1.02 23.72 13.23

NY1 3493 60.79 1.56 III N/A 24.32 13.00

International Journal of Coal Geology 215 (2019) 103296

S. Li, et al.

Fig. 1. Experimental apparatus of oil / water vacuum imbibition tests.

was sealed to avoid any liquid vaporization. The volume changes in the burettes were recorded at different times and converted to sample volume fractions. The tests were terminated when equilibrium had been attained; i.e., no more fluid could be imbibed into the rock sample. The detailed experimental procedure and results can be found in Sang et al. (2018). 3. Model development 3.1. Background As mentioned in the above section, previous studies have clearly shown the existence of superposed organic and inorganic pore systems. Due to the different types of pores, local wettability can vary in shale rocks: water resides primarily within the intraparticle and intergranular pores of the inorganic content, such as clay minerals; the organic material has a strong affinity (i.e., significant molecular interaction) for the hydrocarbon fluids, so part of the oil is stored within the organic compound as dissolved oil, or kerogen-bounded oil (Passey et al., 2010). Therefore, it is reasonable to consider the partitioning of total oil within the matrix to be: 1) free oil within the organic pores; 2) dissolved oil within the organic pores; and 3) free oil within the inorganic pores. In this paper, “oil within the organic pore space” refers to free oil within the organic pore channels, and “oil within the organic material” refers to both free oil within the organic pore channels and dissolved oil in kerogen. In this section, we modified this idealized conceptual model by considering capillary bundles with different wettabilities, thereby making the model more suitable for organic-rich shale rocks. As Fig. 2a and b show, the dual-wet pore network consists of inorganic capillaries with larger radii and organic capillaries with smaller radii. In the model, each phase follows its own path and obeys its own transport mechanisms. From a calculation in Sang et al.'s paper, water had to overcome a capillary pressure of up to 2.8 MPa in order to enter an organic pore with a radius of 1 × 10−7 m (Sang et al., 2018). As a result, water only flows in inorganic pores, as Fig. 2a shows. However, when the sample is pre-vacuumed, oil can flow in both inorganic and organic pores. During the course of oil flowing in the organic pores, it dissolves into the kerogen around the organic pores. The fuzzy boundaries of oil columns in the organic matrix in Fig. 2b illustrate the dissolved oil in kerogen. Similar flowing schemes were proposed in previous studies except that the dissolved oil portion was neglected (Lan et al., 2015; Yassin et al., 2016). In another words, we considered total pore volume as a place for free oil storage, and dissolved oil storage was associated with the volume of kerogen and the volume of oil which can dissolve into a unit volume of kerogen at equilibrium. In this paper, the latter parameter is defined as the term “unit oil intake”. Moreover, we assumed that only free oil can flow in the channels. The

Fig. 2. Schematic illustration of inorganic / organic pore channels. (a) Water imbibition, and (b) Oil imbibition and diffusion.

mobility of dissolved oil inside the kerogen was considered limited as oil molecules were bounded together with kerogen. Besides the dedicated flow paths for oil and water, we have also proposed these following assumptions when formulating the models: 1. Since the spontaneous imbibition process for oil and water was simulated for pre-vacuumed samples, we assumed that part of the invaded fluid evaporated to the gas phase upon entering the vacuumed pores. As a result, there was some diluted gas in pores during the imbibition process. The liquid/vapor interfacial tension (IFT) for water and n-dodecane at room temperature, which were calculated to be 0.036 and 0.020 kg/s2 (36 and 20 dyne/cm) by CMG WinProp module (Version 2016.10, Computer Modeling Group Ltd., Calgary), respectively, were used when calculating the corresponding fluid profile. 2. In the organic pore network, the radii of capillaries were significantly reduced, which led to an increased Knudsen number, Kn. Consequently, the flow regime was expected to shift from continuous flow to Knudsen or transient flow regime (Riewchotisakul and Akkutlu, 2016; Wasaki and Akkutlu, 2015). However, a different opinion came from Fathi: since shale might be over-pressured in the reservoir condition, transient flow was not likely to develop in organic capillaries even though the pore radii were small (Fathi et al., 2012). Hence, the inorganic and organic permeabilities which we introduced in the study were apparent permeabilities: for the inorganic pore network, its apparent permeability equaled its intrinsic permeability as no boundary slip occurred; for the organic pore network, its apparent permeability might deviate from its intrinsic value due to the slippage flow and due to molecular diffusion, and their effects to the deviation were beyond the scope of this study. 3. We also assume that the pore structures remain unchanged during the pre-treatments, i.e. crushing and vacuum, as well as the 3

International Journal of Coal Geology 215 (2019) 103296

S. Li, et al.

imbibition test. As the pore size for the shale is less than one μm and crushed samples are around one to two cm, the crushing only changes the number of pores per sample. Crushing may destroy some of the micro fractures in the core sample while the matrix permeabilities of inorganic and organic pores, which are the focus of the study, are intact. During the imbibition test, kerogen swells as the dissolved oil concentration is increasing in it. But the effects of kerogen expansion was ignored as 1) the highest TOC of the samples used in the study was only 2.84%; 2) the expansion of the coal samples was reported to be less than 5% in the literature (Larsen, 2004); and the swelling length of kerogen for Barnett and Eagle Ford shale was determined to be less than 0.12 nm (Pang et al., 2019). The swelling is significantly small relative to the organic pore radius 8 to 55 nm, and kerogen thickness 54 nm. 4. In the model, we ignored the adsorption of oil and water on the mineral surface of the inorganic pores, because: 1) the clay content of the samples used in the study is not high; 2) the surface adsorption amount is much less than the volume of free oil due to the small surface-area-to-volume ratio of the inorganic pores.

characteristic pore size corresponding to the pressure at which mercury first forms a continuous connected pathway through a sample (Katz and Thompson, 1986). Nelson adapted the form of the equation but chose a different point on the pore-throat-size distribution (35% mercury saturation) as the characteristic pore size (Nelson, 2009):

k

kinorg = 4.48dinorg 2

ninorg =

ninorg rinorg 2L

A

inorg

(2)

rinorg 2

When there is no rock-fluid interaction and fluid flow in capillaries therefore can be considered laminar, all the fluid stays inside the pore channels. At this time, the penetration distance of a liquid into a capillary tube, x, is described using the following equation:

rinorg 2 x

qinorg = ninorg

t

kinorg A P µ x

=

Vkerogen =

+ Patm

x=

(

2 cos rinorg

µ

(3)

Vkerogen = dkerogen =

(4)

)

(5)

The time for fluid breaking through the model, t, is also calculated by transforming the above equation:

µ

t = L2 2kinorg

(

inorg

2 cos rinorg

+ Patm

)

(9)

[(rorg + dkerogen )2 Vkerogen norg L

+ rorg 2

rorg 2] norg L rorg

(10)

(11)

In which norg is the number of organic tubes in a sample, and rorg is the representative radius of orgnaic tubes. We define the outer boundary of the kerogen layer to be rorg + dkerogen as rout. Assume that C(x, r, t) represents the oil concentration at a certain point in the kerogen layer at a certain time, where x indicates the position along the imbibition direction (axial direction), r has the range of rorg < r < rout and indicates the position at a certain thickness of the kerogen (radial direction). and t indicates time. Because the kerogen thickness is far less than the organic tube length, and the diffusion process is much slower than the fluid advancement, it is rational to assume that the primary direction of diffusion is radial and that diffusion in the kerogen along the axial direction is negligible. The diffusion process is described using Fick's law in cylindrical coordinates:

+ Patm t

inorg

rock TOC

where ρrock and ρkerogen are the densities of the rock and of the kerogen, respectively. Assume that kerogen is only wrapped around the organic pores with a homogeneous thickness. The kerogen wall thickness, dkerogen, is then calculated by transforming the equation for the volume of kerogen

where Patm is the atmospheric pressure. After substituting Eq. (4) into Eq. (3), then integrating on both sides, we have the expression of distance x for the inorganic pore system:

2kinorg

AL

kerogen

where qinorg is the imbibition rate in the inorganic tubes, ∆P is the difference of the fluid pressure between the front of the column and the inlet, μ is the viscosity of the fluid, and kinorg is the apparent permeability of the inorganic pore system. It is defined as the permeability of sample with the same size as if there were only inorganic pore networks existing. For vacuum imbibition, we have:

2 cos P = (Pc + Patm) = rinorg

(8)

2

3.2.2. Oil imbibition in organic tubes In the organic pores, as the oil interface advances, oil contacts fresh kerogen and dissolution occurs; further, we assume that the dissolution of the surface layer of oil attains equilibrium instantaneously. Once the concentration gradient is formed along the radial direction inside the kerogen, the dissolution of oil in the kerogen is controlled by a diffusion process. As a result, the oil flow in organic pores can be divided into two parts: interface advancement in the axial direction and diffusion in the radial direction. The arrows in the above Fig. 2b have shown the oil movement in these two directions. The kerogen volume, Vkerogen, is calculated from TOC measurement using the following equation

(1)

AL

inorg

where dinorg is the representative pore diameter for inorganic pore system and ∅inorg is the inorganic porosity. The former is a tuning parameter and the later can be calculated using the final water imbibition volume divided by the sample bulk volume.

3.2.1. Fluid imbibition in inorganic tubes As Fig. 2 shows, in a capillary bundle with cross sectional area of A, straight length of L, in which the inorganic tube radius is rinorg, and the tortuosity of these tubes is τ, we have the following relationships for inorganic porosity ∅inorg and inorganic tube number ninorg:

=

(7)

2

where k is the permeability in millidarcies and d is the characteristic pore diameter in micrometers. As the flow paths along inorganic and organic pore networks are independent and the inorganic and organic pores were superimposed in shale matrix, the Eq. (7) needs to be applied to both inorganic and organic pore systems separately. Specifically for the inorganic pore system, we have:

3.2. Model formulation

inorg

4.48d 2

C D C = r . at any position x t r r r

(6)

(12)

which satisfies the following boundary conditions and the initial conditions:

The flow rate can be calculated by substituting Eq. (5) back into Eq. (3). From percolation theory, Katz and Thompson derived the following relationship between permeability and porosity by introducing the

C (x , r , t ) = Co for 0 < x 4

xt and r = rorg

(13)

International Journal of Coal Geology 215 (2019) 103296

S. Li, et al.

C (x , r , t ) = 0 at r = rout r C (x , r , 0) = 0 for 0 < x

L and rorg

r

rout

(14)

Jo ( rout ) =

(15)

Yo ( rorg )

where xt is the penetration distance and Co is the unit oil intake expressed as cm3 of dissolved oil per cm3 of kerogen. The volume of the dissolved oil was the difference between final oil imbibition volume and the sample pore volume. And the volume of kerogen was calculated from the mass, TOC of the sample, as well as the density of kerogen, which was assigned as 1400 kg/m3. The radial diffusion equation can be solved by the method of separation of variables. First we subtract the equilibrium solution to make the boundary conditions homogeneous:

C =C

=

C D C = r t r r r

=

(19)

After substitution, the initial condition becomes:

Co for 0 < x

L and rorg

r

rout

( rorg ) 4

22

2242

2

rout + 2

ln

rout rorg

an =

(20)

Jo

1/

(21)

2DT (x )

1/

(25)

qt =

where Jo(λr) and Yo(λr) are the first and the second kinds of Bessel functions with order of 0, respectively. After combining the integral constants, the general solution of C is:

Qt =

where a, b, and λ need to be determined based on the boundary and initial conditions:

Jo ( rout )

=

Jo ( rorg ) = 1

22

2

+

2242

4

rorg 6 224262

6

+…

(33)

n r ) rdr

( n r )2 [J0 2

Yo

n

2

{

r )2

( n 2

2(

r )|rrout org

n r ) J1 (

( r )]2 +

( Co ) Yo ( 2(

n

( n r )2 [J0 ( 2

r )]2

}

rout rorg

(34)

n r ) rdr

n r ) rdr n

2(

n r ) Y1 (

[Y0 ( r )]2 +

r )2

( n 2

r )|rrout org [Y0 ( r )]2

}

rout rorg

(35)

exp ( D

n

2t )[a J ( r ) n 0 n

+ bn Y0 (

n r )]

(36)

Qt = t

xt 0

xt

t

0

0

2 rorg

2 rorg

D

D

C (x , rorg , t ) r

C (x , rorg , t ) r

dx

dxdt

(37) (38)

rorg 2 x

norg =

t

A

=

korg A P µ x

In this equation, the left-hand side of the equation represents the incremental oil which remains in the pore channel, the first term on the right-hand side of the equation is the total imbibed oil per unit of time. korg is the organic permeability, which is defined as the permeability of sample with the same size as if there were only organic pore networks existing. Similar to Eq. (2), we have:

(29)

rorg 4

…

(28)

The expressions of J0, Y0, J0 , and Y0 can be found below in Eqs. (30) to (33).

rorg 2

2 2 J0 ( rout ) + rout

norg

Yo ( rorg ) Yo ( rout )

2

(27)

In order to find A and B which can satisfy both Eqs. (27) and (28) simultaneously, we need to have:

Jo ( rorg )

(32)

where Qt is the cumulative absorption at time t. The dissolved oil does not stay in the pore space, so the fluid which pushes the moving interface forward is less than the total imbibed fluid at the inlet:

(26)

aJo ( r ) + bYo ( r ) = 0 at r = rout

…

Once the oil concentration at each location of kerogen is known, the instantaneous rate of oil dissolving into kerogen, qt, is written as:

Also note that Eq. (23) is the Bessel differential equation with order of 0; the general solution is:

aJo ( r ) + bYo ( r ) = 0 at r = rorg

( rorg )6 1 1 1 + 2 2 2 1+ + 2 246 2 3

1+

n=1

(24)

C = Co + exp ( D 2t )[aJo ( r ) + bYo ( r )]

{

C = Co +

Note that Eq. (22) has the general solution of:

R (r ) = AJo ( r ) + BYo ( r )

(31)

After combiningan, bn, and λn, the final solution is:

(23)

T (t ) = Cexp ( D 2t )

+…

n r ) rdr

( Co )1/

=

and

r 2R (r ) + rR (r ) + r 2 2R (r ) = 0

2

rinner router rinner

bn =

(22)

=0

n

router

This equation only makes physical sense when both terms are equal to the indicated negative constant −λ2. Thus, the original problem has now been transformed into two boundary value problems, namely:

2(

( Co )1/

=

1

T (t ) +

8

2

J0 ( rout ) +

( Co ) Jo ( rout

rorg

Then we assume that C = R(r)T(t), which, after substituting this expression into Eq. (17), yields: 2

rout 7 2242628

6

( rout )3 1 ( r )5 1 1 1+ + 2out2 1+ + 224 2 2 46 2 3

∗

R (r ) + r R (r ) T (t ) = = DT (x ) R (r )

rout 5 22426

Where γ is Euler's constant. Once we substitute the expressions for Jo, Yo, Jo , and Yo back into Eq. (29), it becomes a polynomial equation of λ, and λ is the only unknown. We define λn is the nth positive roots of λ. Determination of the coefficients of an and bn, which correspond to each λn by applying the Sturm-Liouville theory to the radial regime with kerogen, we have:

(18)

C (x , r , t ) = 0 at r = rout r

+

J0 ( rorg ) +

( rorg )2

rout 2

which is then bounded by:

xt and r = rorg

+

2

4

Yo ( rout )

(17)

C (x , r , t ) = 0 for 0 < x

rout 3 224

2

rorg

ln

(16)

Co

C (x , r , 0) =

2

rout 2

norg qt

(39)

org

rorg 2

(40)

where ∅org is the organic porosity of the sample which can be calculated from the pore volume and final oil imbibition volume. Similar to

(30) 5

International Journal of Coal Geology 215 (2019) 103296

S. Li, et al.

process for both water and oil. The detailed procedure for conducting the simulation is listed: (1) We determined these following parameters using lab measured results prior to simulation: a. Inorganic porosity was determined using the final water imbibition volume and sample bulk volume; b. Organic porosity was determined using the final oil imbibition volume, sample bulk volume, and sample pore volume. c. Dissolved oil volume was determined using the final oil imbibition volume and sample pore volume; d. Kerogen volume was determined using TOC, sample mass, and kerogen density; e. Equilibrium dissolved oil concentration was determined using dissolved oil volume and kerogen volume. (2) Water imbibition profile was matched using Eq. (5) by assigning different inorganic pore radii. With each inorganic pore radius, inorganic tube numbers ninorg, was calculated according to Eq. (2). Inorganic permeability was then calculated using Eq. (8) and inorganic pore radius which reached the best fit. (3) The oil imbibition in inorganic pores was calculated using inorganic pore radius from Step (2). After that oil imbibition profile in organic pores was calculated by deducting oil intake into the inorganic pores from the total oil imbibition curve measured in the lab. (4) Oil imbibition profile in organic pores was then matched. It was much more complicated than simulating water profile as it incorporated both flow and diffusion. The flow chart for calculating the imbibition profile in the organic pores can be found in Fig. 4 below.

Fig. 3. Illustration of time-dependent dissolution calculation with changing equilibrium dissolution.

Eq. (8), we have the same relationship between porosity and permeability in the organic pore system:

korg = 4.48dorg 2

org

(41)

2

where dorg is the representative pore diameter for organic pore system and it is a tuning parameter. As the interface moves forward, the oil contacts more kerogen and the total equilibrium/maximum dissolved oil, Qmax, t, increases correspondingly. The below Eq. (42) expresses the equilibrium dissolved oil corresponding to a given time t and interface position xt:

Qmax, t =

(rout 2

rorg 2)

xt 0

xCo dx

At first, we start from a time step t with known equilibrium dissolved oil volume Qmax, t, instantaneous dissolved oil volume Qt, and interface distance xt. In order to find the updated interface distance xt+∆t, we use a trial and error method by adding a small incremental distance ∆x:

(42)

With this set of conditions, we have developed a diffusion problem with changing boundary conditions. In order to solve it, a new solution approach that we designate the “equivalent time method” is proposed in this study. Fig. 3 below illustrates the calculation process. The blue and orange curves are the volumes of dissolved oil versus time, with the maximum dissolved oil volumes defined as Qmax, t and Qmax, t+∆t, respectively. Qt is the instantaneous dissolved oil at time t. Below Fig. 3 is the step–by–step procedure to calculate instantaneous dissolved oil at time t + ∆t, with known Qt, Qmax, t, and Qmax, t+∆t.

x t(+i)

(46)

= xt + i x

The superscript (i) designates the ith iteration within the time step. The incremental distance ∆x is sufficiently small at every iteration to make sure that stability is satisfied for the solution process. For the known interface position xt+∆t(i), these following terms were calculated explicitly:

(1) At time step t, with known Qmax, t, Qt is a function of Qmax, t, D, and t, and the analytical solution is given by Eqs. (36) to (38).

a. The total imbibed volume from time t to t + ∆t using Darcy's law:

(43)

Qt = f (Qmax , t , D, t )

t

(2) When the maximum absorption changes to Qmax, t+∆t, the equivalent time t∗ for the corresponding instantaneous Qt is found. In Fig. 3, a horizontal line at V = Qt was drawn and it crossed the orange line at the point (t∗, Qt). The explicit formula of t∗ is difficult to write and therefore expressed as the inverse function of Qt. In the calculation, t∗ was obtained by trial and error.

t =f

(44)

1 (Q max , t + t , D, Qt )

(3) The incremental absorption is then calculated by proceeding to the next time step along the orange line. Correspondingly, in the below equation, the updated equilibrium volume of oil dissolution Qmax, ∗ t+∆t and new time step t + ∆t are used.

Qt +

t

= f (Qmax , t + t , D, t +

t)

(45)

After introducing the “equivalent time method” to calculate the instantaneous adsorption change within one time step, we need to incorporate it into the entire workflow to simulate the full imbibition

Fig. 4. Flow chart of model computations. 6

International Journal of Coal Geology 215 (2019) 103296

S. Li, et al.

Vt(+i)

=

t

korg A P t µ xt(+i ) t

number, norg, was updated according to Eq. (40), and the thickness of kerogen layer was calculated from Eq. (11).

(47)

4. Results and discussion

b. The updated equilibrium volume of dissolved oil Qmax, t+∆t(i) corresponding to xt+∆t(i) using Eq. (42). c. Instantaneous dissolved volume Qt+∆t(i) corresponding to Qmax, (i) using Eq. (45) in the “equivalent time method.” t+∆t d. Then we verify whether the below equation is satisfied:

Vt(+i)

= norg [(Qt(+i )

t

t

Qt ) + (xt(+i )

(48)

x t ) rorg 2]

t

The proposed model in Section 3 has been used to calculate the oil and water interface movement and imbibition volumes for shale rock samples. The calculation started from the time when oil and water first contacted the inlet of the pre-vacuumed sample, and ended when the fluid had completely saturated the pore spaces and that the diffusion had attained equilibrium, if applicable. In Section 4.1, the model is shown to have obtained the inorganic and organic permeabilities of four shale rock samples by matching the vacuum imbibition profiles. Moreover, there is a summary of the inorganic/organic permeabilities, porosities, and fluid saturations for the four shale rock samples. In Section 4.2, model calculation results for rock sample #1 are discussed in detail, and different states of the imbibed oil are differentiated. In Section 4.3, the sensitivities of the organic permeability, TOC, and diffusion coefficient in the calculation for Sample #1 are discussed.

xt+∆t(i)

If yes, we can then treat the as the correct interface at the end of the time step t + ∆t. Otherwise, additional iterations are required by keeping on adding other incremental distances, ∆x, as Eq. (47) shows, to the assumed interface position at the end of the time step until Eq. (48) is satisfied. Once the converged correct position is determined, the equilibrium and instantaneous volumes of oil dissolution are updated before moving on to the next time step.

Qmax, t + Qt +

t

t

(i ) = Qmax , t+

= Qt(+i)

(49)

t

4.1. Model validation

(50)

t

Fig. 5 below represents the experimental and simulation results for oil/water spontaneous imbibition tests for four shale rock samples. In these figures, the blue dots show the imbibed water volume percentage converted from measured water volume, the red dots show the imbibed

The Step (4) was repeated with different organic pore radii until the best fit for oil imbibition into organic material was found. Then the organic permeability can be calculated based on Eq. (41). Every time when different organic pore radii were assigned, the organic tube 20

20

1

2 16

Volume percentage

Volume percentage

16

12

8

Simulated water Measured water Simulated oil Measured Oil

4

12

8

Simulated water Measured water Simulated oil Measured Oil

4

0

0 0

300

600

900

1200

1500

0

300

Time, hour

600

900

1200

20

20

3

4 16

Volume percentage

16

Volume percentage

1500

Time, hour

12

8

Simulated water Measured water Simulated oil Measured Oil

4

12

8

Simulated water Measured water Simulated oil Measured Oil

4

0

0 0

400

800

1200

1600

2000

0

Time, hour

300

600

900

Time, hour

Fig. 5. Experimental and simulation results for water and oil vacuum–imbibition test using Samples 1 to 4. 7

1200

1500

International Journal of Coal Geology 215 (2019) 103296

S. Li, et al.

corresponding organic permeabilities from 2.24 × 10−6 md to 1.59 × 10−5 md. (3) One sees that the inorganic permeability is approximately 200 to 3000 times greater than the organic permeability for each sample, consistent with the saturation time and radii comparisons of inorganic and organic pores. It must be stated that not all the pores in the inorganic or organic systems share the same pore size and pores with different sizes have different imbibition rates. The effects of the pore size to the imbibition profile will be discussed later in 4.3.1. The ideal and most accurate way to simulate the overall imbibition profile is to run the model as many times as possible for different pore sizes and then add them together based on the pore size distribution. However, the pore sizes distribution for shale samples were hard to obtain as mercury injection is not suitable for pores smaller than 2 nm whereas nitrogen adsorption cannot characterize pores larger than 50 nm (Barnett et al., 2012; Nelson, 2009); and either these two methods can separate the organic pores from the inorganic pores. In order to simplify the calculation, one representative pore radius for each inorganic pores and organic pores were selected to simulate the fluid imbibition profile. The reason to call them the representative inorganic/organic pore radius is that by selecting appropriate pore radius which achieves the best-fit of the imbibition profiles, one can directly calculate the inorganic/organic permeability based on Eqs. (8) and (41) and the yielded permeability values can represent the flow ability of the inorganic and organic pore networks in the samples. The organic porosity, inorganic porosity, organic fluid saturation, and inorganic fluid saturation are defined and calculated using the below Eqs. (51) to (54), respectively.

10 (a) Oil in inorganic pores

Volume percentage, %

8

Oil inorganic pores Dissolved oil

6

4

2

0 0

5

10

15

20

25

30

35

40

Time, Hour 10 (b)

Volume percentage, %

8

6

4 org

Oil in inorganic pores

=

PVorg

2

Dissolved oil

inorg

0 0

200

400

600

800

1000

1200

=

Sorg =

1400

PVinorg (52)

Vbulk PVorg + Vdis

PVorg + PVinorg + Vdis

Time, Hour Sinorg =

Fig. 6. Model calculated imbibed volume percentage for oil versus time for one shale sample, for (a) first 40 h, (b) the whole process.

(51)

Vbulk

Oil inorganic pores

(53)

PVinorg PVorg + PVinorg + Vdis

(54)

where dissolved oil volume is represented by the difference between total oil volume and helium pore volume, and organic pore volume is calculated by the difference between helium pore volume and water pore volume. Table 2 below summarizes the inorganic/organic porosities and fluid saturations for the four shale rock samples, along with their inorganic/organic permeabilities which were obtained by matching the experimental imbibition curves for oil and water saturation tests. It can be seen that the maximum dissolution of oil accounts for 0.70 to 1.96% of the total bulk volume for the four samples and that the dissolved oil also causes the discrepancy between porosity and fluid saturation. From Samples #1 to #4, the organic pores only account for 24.2%, 24.2%, 33.3%, and 28.8% of the total pore volume, while 40.6%, 34.2%, 36.7%, and 33.9% of oil is stored in the organic pore systems, respectively. This restates again the importance of accounting for the dissolved oil and the possible significant underestimation of original oil in place if dissolved oil is ignored during production.

oil volume percentage, and the black solid and dash lines are the simulated water and oil volume intake curves, respectively, using the model described in Section 3.2. We see that a satisfactory match has been achieved for both water and oil imbibition for all four cases, corroborating the validity of the model and its parameters. The vacuum–imbibition results of all shale oil rock samples show a similar pattern: (1) The equilibrium volume of water is always less than that of oil. The scales of all the components in Fig. 6 have been set to the same to show the trend clearly. The reason was stated above when developing the model: water can only occupy the inorganic pore spaces in shale rocks, while there are three forms of oil in the shale pore systems: free oil in inorganic pores, free oil in organic pores, and dissolved oil in kerogen. When comparing the above results with the helium measured porosity data in Table 1, one sees that the helium saturation volume sits between the water volume and the oil volume, consistent with the fact that the helium volume represents the total effective pore volume of both inorganic and organic media, but the dissolved oil volume is not included. (2) The saturation of oil requires a much longer time to equilibrate than that of water. This is due to the different sizes of the inorganic and organic pores. For the four samples, the representative inorganic pore radii range from 1.9 × 10−7 m to 4.0 × 10−7 m, with the corresponding inorganic permeabilities from 2.35 × 10−3 md to 1.37 × 10−2 md; the representative organic pore radii range from 8.0 × 10−9 m to 5.5 × 10−8 m, with the

4.2. Model results Fig. 6 below shows the simulation results for an oil spontaneous imbibition test for shale rock sample #1 in Table 1. The determined inorganic and organic permeabilities were 8.33 × 10−3 md and 1.59 × 10−5 md, with corresponding inorganic and organic pore radii of 4 × 10−7 m and 5.5 × 10−8 m, respectively. The diffusion coefficient was determined to be 5 × 10−22 m2/s. It can be found that the 8

International Journal of Coal Geology 215 (2019) 103296

S. Li, et al.

Table 2 Summarized petro-physical properties of samples. Number

#1

#2

#3

#4

Sample name Water volume (ml/g) Helium volume (ml/g) Oil volume (ml/g) Inorganic porosity (% sample V) Organic porosity (% sample V) Helium porosity (% sample V) Max dissolution (% sample V) Total imbibed oil (% sample V) Inorg. fluid saturation (% pore V) Org. fluid saturation (% pore V) Representative inorg. radius (m) Representative org. radius (m) Inorganic permeability (md) Organic permeability (md) K ratio (Kinoganic/Korganic)

L67-3307 0.022 0.029 0.037 5.39 1.72 7.11 1.96 9.07 59.46 40.54 4.0 × 10−7 5.5 × 10−8 8.33 × 10−3 1.59 × 10−5 524

L67-3276 0.025 0.033 0.038 6.03 1.93 7.96 1.21 9.17 65.79 34.21 1.9 × 10−7 4.0 × 10−8 2.35 × 10−3 1.06 × 10−5 243

LY1-3624 0.038 0.057 0.060 8.82 4.41 13.23 0.70 13.93 63.33 36.67 2.1 × 10−7 8.0 × 10−9 6.14 × 10−3 2.24 × 10−6 2741

NY1-3493 0.037 0.052 0.056 9.25 3.75 13.00 1.00 14.00 66.07 33.93 3.0 × 10−7 2.0 × 10−8 1.37 × 10−2 1.00 × 10−5 1370

concentration gradient still exists in the kerogen along the radial direction. Under these circumstances, part of the oil is transported from the organic pore spaces into the deeper kerogen by diffusion. Due to the very small diffusion coefficient and low organic permeability, oil dissolution into kerogen and oil moving in organic pores occur in a very slow manner. As a result, the overall oil imbibition rate in the second stage is much less than in the first stage. The last stage starts after oil fully saturates the organic tubes. In this stage, both inorganic and organic pores have already been fully saturated with oil so that the only oil intake is that of oil diffusion into kerogen; i.e., we still observed slightly more oil intake at this stage. The value for the diffusion coefficient used in this study is 5 × 10−22 m2/s, which may seem to be very low, so we have compared it to a literature value. Yang et al. simulated the gas diffusion process into a spherical kerogen particle (Yang et al., 2016). The parameters used in that study include: an apparent diffusion coefficient of 6.2 × 10−12 m2/s, a particle radius of 200 μm, and a characteristic time for the entire diffusion process to attain equilibrium of 1 h. Here, we proposed a dimensionless time parameter, tD = Dt/r2, to describe the time required for a diffusion process to attain equilibrium. In the process simulated in Yang's study, tD is calculated to be 0.56. In this study, specifically for sample #1, if we take the characteristic time to attain equilibrium as 1000 h and a kerogen thickness of 5 × 10−8 m as calculated from the TOC, the dimensionless time from this study is calculated to be 0.72. Therefore, close dimensionless times were found between this study and previous research despite the diffusion coefficient in these two studies being significantly different. From these two studies, it should also be noted that the diffusion coefficient for liquid into solid is 10 orders of magnitude less than that for gas into solid.

clear advantages of the proposed model over previous experimental study done by Sang et al. include: (1) different states of imbibed oil can be easily differentiated instead of only the total imbibed oil volume as recorded in the experiments; (2) instead of using the final imbibition values for oil and water to determine the inorganic/organic porosities and saturations, this study utilized the time dependent imbibed volumes and incorporated them with the proposed diffusion-imbibition model so that the inorganic/organic permeability can be obtained by matching the imbibition profile (Sang et al., 2018). In Fig. 6, the total imbibed volume has been divided into three parts and colored separately. The blue area represents the volume of oil stored in the inorganic pores in terms of percentage of the total sample volume; the orange area represents the oil volume percentage entering into the organic pores; and the green area represents the dissolved oil volume converted to the percentage of the total sample volume. The final bulk volume percentages of the three states of oil are 5.39%, 1.72%, and 1.96%, respectively. It can be seen that, although the majority of the pore space is inorganic, due to the existence of kerogen around the organic tubes, rock-fluid interaction helps to absorb a significant amount of oil within the organic pore system. Specific to this sample, oil in the organic material accounts for 40.6% of the total imbibed oil, yet the organic pores account for only 24.2% of total pore volume. As a result, if the dissolved oil is neglected, significant discrepancy is expected during reservoir simulation and the original oil in place and production can be significantly underestimated. It is also worth mentioning that the distribution among the three states of oil can vary and depends on other factors, such as kerogen volume. One of the findings from Sang et al.'s research is that a greater amount of dissolved oil was expected for shale rock samples with higher TOCs, and the 20 samples had dissolved oil ranging from 4.3 to 47.1%, with an average of 16.2% of the total oil volume (Sang et al., 2018). From the temporal domain, the imbibition process is also divided into three stages, which are bounded using vertical black dashed lines in Fig. 6. It can be seen that the first stage lasts approximately 20 h. At this stage, most of the imbibed oil fills the inorganic pores. This is logical because inorganic media have larger pore sizes and greater permeabilities. The first stage ends when oil fully saturates the inorganic tubes. The second stage starts at 20 h and extends to 1150 h. Within this stage, as inorganic pores cannot imbibe more oil, the blue area in Fig. 6 does not increase. Oil instead fills the organic pores and part of the oil is also dissolved into the kerogen as dissolved oil. The driving force for oil dissolving into kerogen comes from two mechanisms: (1) as oil moves forward in the organic pores, it contacts more “fresh” kerogen which has not been contacted previously with oil; as a result, surface adsorption occurs, followed by oil dissolution into the deeper kerogen; (2) although the kerogen behind the fluid interface has already been contacted with oil, equilibrium has not been established and an oil

4.3. Sensitivity analysis After examining the typical results for the model and comparing the model results with experimental data, different mechanisms have been shown for the oil flowing in the inorganic and organic pore systems: the former is controlled by the inorganic permeability, and the latter is controlled by the organic permeability and the diffusion coefficient, as well as the kerogen content. The conventional method to calculate fluid flow using only one permeability and one porosity clearly does not reasonably represent the flow mechanisms and therefore is not applicable. The sensitivities of these critical petro-physical parameters in this study have to be investigated. We used the results from Sample #1 as the example to examine how a change of these variables affects the interface advancement and the fluid saturation curve in our model. Because the solution for inorganic pore systems was obtained analytically, the impact of the variables can be directly observed from the equations; therefore, only the imbibition volume and interface 9

International Journal of Coal Geology 215 (2019) 103296

S. Li, et al.

Table 3 List of Parameters Used in Sensitivity Analysis 4.3.1. Organic perm, md −3

1 2 3 4 5

1.59 × 10 1.59 × 10−4 1.59 × 10−5 1.59 × 10−6 1.59 × 10−7

Representative pore size, m

Kerogen thickness, m

−7

−7

5.5 × 10 1.7 × 10−7 5.5 × 10−8 1.7 × 10−8 5.5 × 10−9

5.4 × 10 1.7 × 10−7 5.4 × 10−8 1.7 × 10−8 5.4 × 10−9

Diffusion coefficient, m2/s

0.015

−20

5 × 10 5 × 10−21 5 × 10−22 5 × 10−23 5 × 10−24

Interface position, m

Case No.

0.02

movement profile in the organic tubes are discussed in this part of study. The variables selected for analysis include organic permeability, TOC, and diffusion coefficient. 4.3.1. Effects of organic permeability In this section, four different organic permeabilities, 1.59 × 10−7, 1.59 × 10−6, 1.59 × 10−4, and 1.59 × 10−3 md, along with the original value, 1.59 × 10−5 md, were used in the model. The change of the organic permeability was achieved by adjusting the representative pore size. To increase the organic permeability from 1.59 × 10−7 md to 1.59 × 10−3 md, the pore size was increased from 5.5 × 10−9 m to 5.5 × 10−7 m. In order to keep the effects of the diffusion the same for different pore sizes and kerogen thicknesses, the diffusion coefficients was changed from 5 × 10−24 md to 5 × 10−20 md accordingly. Below the Table 3 listed corresponding pore sizes, kerogen thickness, as well as the diffusion coefficients for each cases. The calculated imbibed volume percentage and interface movement profiles are plotted in Figs. 7 and 8, respectively. From these two figures, it can be seen that the greater the organic permeability i.e. the larger pore size, the greater the volume of oil that is imbibed within a certain time, and the faster the interface moves. When the organic permeability was set to 1.59 × 10−3 md and 1.59 × 10−4 md, the time when the oil volume attains equilibrium significantly fell behind the time when the oil interface reaches the end of the model. It can be found from Eq. (47), when organic permeability was increased, more oil could enter into the organic pore system per unit of time. The extra oil tended to stay in the pore channels and push the interface forward, as the effects of diffusion were the same. This phenomenon was also observed by previous study in which the LucasWashburn equation was used to calculate the interface movement versus time for channels with different pore sizes in a porous media (Shi

Korg = 1.59*10-4 md

Korg = 1.59*10-5 md

0

Korg = 1.59*10-6 md

Organic saturation, %

0 600

800

1000

1200

1400

800

1000

1200

1400

1600

4.3.2. Effects of TOC, the same unit oil intake In this study, TOC, which is a measurement of kerogen volume, is represented by the thickness of kerogen around the organic tubes in the model as the organic porosity was fixed constant. The unit oil intake (cm3 of dissolved oil per cm3 of kerogen) is assumed to be constant so that the total equilibrium dissolved oil is proportional to TOC; i.e., samples with higher TOC contain more dissolved oil. This is a valid assumption because the properties of kerogen are similar among samples with similar lithologies. It has also been found by Sang et al. that greater TOC leads larger dissolved oil content by analyzing imbibition results for 20 shale samples (Sang et al., 2018). Figs. 9 and 10 exhibit the profiles of imbibition oil volume and interface movement, respectively. From these two figures, it can be seen that the greatest TOC corresponds to the greatest dissolved oil volume and the slowest interface movement. As the organic pore radius, i.e. the organic permeability, was the same among these cases, the total imbibition rate for oil entering the organic pore system was similar. However, for a greater kerogen volume which can hold more dissolved oil, there is a lower dissolved oil concentration and it, in turn, creates higher flux for oil dissolved into kerogen from organic pore channels. The oil stays in the pore channel, which is the difference between total imbibition rates and dissolution rate, is less when TOC is higher and it corresponds to slower interface movement. For the cases with TOCs equal to 0.57 wt%, 1.42 wt%, and 2.84 wt%, the equilibrium volume percentages are 2.11%, 2.69%, and 3.65% of bulk sample volume, respectively. After deducting the volume of organic pore space, 1.72%,

1

400

600

et al., 2018). As a result, the interface reached the end of the model at an earlier time and it took a much longer time for diffusion to attain equilibrium. This is also the reason for the similar patterns in Figs. 7 and 8 when the organic permeabilities are 1.59 × 10−6 md and 1.59 × 10−7 md. When the organic permeability was small, diffusion was able to consume all of the imbibed oil in the organic pore channels, causing a very slow interface movement. The kerogen behind the interface almost attained equilibrium, and we can sum the interface movement and diffusion together, treat the case with instantaneous dissolution as a case without diffusion but with large pore space. The relationship between pore size and imbibition rate can be extended to the inorganic pore system as well. When there is no diffusion involved, Eq. (5), which is used to calculate the interface movement for the inorganic pores, has exactly the same form as it in Shi's paper, although the flow directions in these two studies are different (Shi et al., 2018).

2

200

400

Fig. 8. Simulation results showing interface position versus time under different organic permeabilities.

3

0

200

Times, hour

Korg = 1.59*10-7 md

4

Korg = 1.59*10-3 md Korg = 1.59*10-4 md Korg = 1.59*10-5 md Korg = 1.59*10-6 md Korg = 1.59*10-7 md

0.005

0

5

Korg = 1.59*10-3 md

0.01

1600

Times, hour

Fig. 7. Simulation results showing imbibed volume under different organic permeabilities. 10

International Journal of Coal Geology 215 (2019) 103296

S. Li, et al.

6

5

TOC = 14.20% TOC = 5.68% TOC = 2.84% TOC = 1.42% TOC = 0.57%

D = 5*10-18 m2/s D = 5*10-20 m2/s D = 5*10-22 m2/s

4

D = 5*10-24 m2/s D = 5*10-26 m2/s Organic saturation, %

Organic saturation, %

4

2

3

2

1

0

0 0

200

400

600

800

1000

1200

1400

0

1600

200

400

600

Times, hour

1000

1200

1400

1600

Fig. 11. Simulation results showing imbibed volume under different diffusion coefficients.

Fig. 9. Simulation results showing imbibed volume under different TOCs but constant unit oil intake of kerogen.

0.02

0.02

0.015

Interface position, m

0.015

Interface position, m

800

Times, hour

0.01

TOC = 14.20% TOC = 5.68% TOC =2.84% TOC = 1.42% TOC = 0.57%

0.005

0.01

D = 5*10-18 m2/s D = 5*10-20 m2/s D = 5*10-22 m2/s D = 5*10-24 m2/s D = 5*10-26 m2/s D=0

0.005

0

0 0

200

400

600

800

1000

1200

1400

0

1600

200

400

600

800

1000

1200

1400

1600

Times, hour

Times, hour

Fig. 12. Simulation results showing interface position versus time under different diffusion coefficients.

Fig. 10. Simulation results showing interface position versus time under different TOCs but constant unit oil intake of kerogen.

this value is due to dissolution. The figure indicates that, when the diffusion coefficient equals 5 × 10−26 or 5 × 10−24 m2/s, the imbibition volume is significantly less than the cases with larger diffusion coefficients. However, the original case with a diffusion coefficient of 5 × 10−22 m2/s shows a similar profile to the cases with diffusion coefficients which are 100 and 10,000 times larger than the original value. The corresponding green, red, and purple dashed lines show that these cases attained the equilibrium total organic oil saturation of 3.64% of sample bulk volume, while the blue and black dashed lines representing the cases with lower diffusion coefficients did not attain equilibrium at the end of the simulation. Moreover, it can be seen that, when the diffusion coefficient equals 5 × 10−22 or 5 × 10−24 m2/s, the imbibition profiles are essentially identical from the start. It is therefore logical to conclude that, once the diffusion coefficient exceeded 5 × 10−22 m2/s, the diffusion process attained equilibrium essentially instantaneously and that higher diffusion coefficients did not further facilitate the oil dissolution process. In Fig. 12, an orange dashed line, representing the interface position versus time when there is no rock-fluid interaction — i.e., no diffusion— was added to compare with cases with different diffusion

from the volume of oil in the organic pore system, the dissolved oil volumes for the three cases are 0.39%, 0.97%, 1.93%, respectively, which are proportional to the respective TOCs. We do not observe these relationship for cases with greater TOCs as equilibrium has not been attained for these cases at the ends of simulation. In below Fig. 9, organic saturation still increases after 1500 h when TOC is 14.20% and 5.68%. It can be seen that, for cases with greater TOCs, some dissolved oil is stored in the kerogen far away from the organic pores. As a result, greater diffusion coefficients would be required to guarantee sufficient mass transfer over the durations of these simulations. 4.3.3. Effects of diffusion coefficient In this section, four different diffusion coefficients, 5 × 10−18 m2/s, 5 × 10−20 m2/s, 5 × 10−24 m2/s, and 5 × 10−26 m2/s, were simulated, along with another case with original diffusion coefficient of 5 × 10−22 m2/s. The imbibed volume percentages are plotted against time in Fig. 11, and the interface movements are reflected in Fig. 12. In Fig. 11, a horizontal black solid line at y = 1.72 is plotted to represent the organic porosity. The organic imbibed volume beyond 11

International Journal of Coal Geology 215 (2019) 103296

S. Li, et al.

coefficients. It can be seen that, when a diffusion coefficient is extremely small, say 5 × 10−26 m2/s, the movement of the interface is almost identical from the case without diffusion. This means that, when the diffusion coefficient is small, the imbibed fluid is more likely to stay within the pores because it takes time for it to transfer into the kerogen, but the newly imbibed fluid keeps filling in. As a result, less time was required for oil to fully fill the organic pore spaces with a smaller diffusion coefficient. Another phenomenon which can be observed in Fig. 12 is that, in the early stage, the interface with the diffusion coefficient of 5 × 10−22 m2/s, represented by the green dashed line, travels faster than the interfaces of the cases with greater diffusion coefficients; however, in the later stage, it became slower. The reason is that, at the early stage, the oil volume dissolved into the kerogen was less for the run with a smaller diffusion coefficient. In another words, imbibed oil was more likely to stay within the pore space and the interface accordingly moves faster. However, the faster interface movement helped the oil contact more ‘fresh” kerogen, which in turn accelerated the oil dissolution. Another factor that changes the profile in the early and later stages is that the imbibition rate reduced as the imbibed length increased. The case with faster interface movement in the early stage experienced more significant rate drop. Taken together, these reasons caused the reverse of the interface profiles from early to late stages during the entire imbibition process.

both the organic permeability and the diffusion coefficient, as well as the kerogen content. The effects of these parameters have been investigated in this study as well. This study emphasizes the importance of considering shale rock as a coupled system with both inorganic and organic compartments when characterizing the fluid flow in a shale reservoir. The organic permeability is hundreds to thousands of times smaller than the inorganic permeability. Although it has been a widely used past practice, using an average permeability or effective permeability cannot adequately describe the flow of shale oil properly because different mechanisms apply and must be considered when oil flows in the dual organic and inorganic pore systems. Our future study will focus on modeling the recovery process in shales under this premise as well as how to adapt and apply this concept in a large scale reservoir simulation. From the production perspective, understanding under which conditions the dissolved oil can be activated and how to facilitate oil flow in the organic pore networks with very low permeability is of practical significance. Nomenclature a an A b bn C Co C(x, r, t) C∗

Integral constant for first kind Bessel function Coefficient of A corresponding to λn Cross-sectional area of the rock sample Integral constant for second kind Bessel function Coefficient of B corresponding to λn Dissolved oil concentration Equilibrium/maximum oil concentration Dissolved oil concentration at x, r, t Concentration of dissolved oil which deviates from its equilibrium value C∗(x, r, t) Concentration of dissolved oil which deviates from its equilibrium value at x, r, t d Diameter of representative pores dinorg Diameter of representative inorganic pores dkerogen Thickness of the kerogen layer around organic pores dorg Diameter of representative organic pores D Diffusion coefficient of oil in kerogen i Iteration steps Jo First kinds Bessel function with order of 0 Jo First order derivative of Jo kinorg Inorganic permeability korg Organic permeability L Rock sample length ninorg Numbers of the inorganic tubes in a rock sample norg Numbers of the organic tubes in a rock sample Atmospheric pressure Patm Pc Capillary pressure PVinorg Volume of inorganic pores PVorg Volume of organic pores Imbibition rate in the inorganic pores qinorg qt Instantaneous rate of oil dissolving into kerogen Equilibrium/maximum absorption volume at t Qmax, t Qmax, t+∆t Equilibrium/maximum absorption volume at t + ∆t Qt Cumulative absorption volume at t Qt+∆t Cumulative absorption volume at t + ∆t r Radius indicator in the diffusion equation rinorg Radius of inorganic tubes Radius of organic tubes rorg rout Outer radius of kerogen layer R(r) Component of radius function in the solution of C∗ R′(r) First order derivative of R(r) R′′(r) Second order derivative of R(r) Sinorg Inorganic fluid saturation Sorg Organic fluid saturation t Time indicator

5. Conclusions In this study, a comprehensive methodology was proposed to determine the inorganic and organic permeabilities of shale and it revealed the microscopic storage features and flow mechanisms for shale rocks. The methodology incorporated the water and oil vacuum–imbibition tests in shale rock samples with a bundle–of–tube model in which the diffusion process and the flow process were coupled together. The algorithm to solve the model's constitutive equations was developed, and the model was then used to successfully match the results of the vacuum–imbibition tests using four shale rock samples. After that, a series of sensitivity analyses were conducted. The main conclusions are listed below: (1) The lab tests exhibited the following phenomena for all four shale rock samples: imbibed oil volume > gas pore volume > imbibed water volume. The reason is that the oil can enter the pore space within kerogen and as well it can dissolve into the kerogen while water cannot. Therefore, the imbibed water volume, the difference between the gas pore volume and the imbibed water volume, and the difference between the imbibed oil volume and the gas pore volume indicate the inorganic pore space, organic pore space, and dissolved oil volume, respectively. (2) The four shale rock samples used in this study have the following parametric ranges: inorganic porosities of 5.39 to 9.25% of bulk sample volume; organic porosities of 1.72 to 4.41% of bulk sample volume; inorganic saturations of 59.46 to 66.07%; organic saturations of 33.93 to 40.54%; inorganic permeabilities of 2.35 × 10−3 to 1.37 × 10−2 md; and organic permeabilities of 2.24 × 10−6 md to 1.59 × 10−5 md. Each sample's organic permeability is 243 to 2741 times less than its corresponding inorganic permeability. From the above values, it can be seen that a significant amount of oil is stored in the organic pore systems, but the low organic permeabilities make this part of the oil difficult to access. This was also observed in the vacuum–imbibition test: it took more than 1000 h for oil to completely saturate the organic pore system. (3) The first of its kind proposed model considered both the inorganic and organic permeabilities when simulating shale oil flow. Imbibition in the inorganic pore system dominates the early stage. Oil dissolves into the kerogen once it enters the organic pore systems, which slows the filling process in the organic pores. As a result, the oil flowing in the organic pore system is controlled by 12

International Journal of Coal Geology 215 (2019) 103296

S. Li, et al.

t∗ TOC T(t) T′(t) Vbulk Vdis Vkerogen x xt Yo Yo

Equivalent time when instantaneous absorption is Qt and equilibrium/maximum absorption is Qmax, t+∆t Total organic content Component of time function in the solution of C∗ First order derivative of T(t) Bulk volume of a shale sample Dissolved oil volume Volume of kerogen Position indicator Penetration distance of a liquid into a capillary tube at time t Second kind Bessel function with order of 0 First order derivative of Yo

Petroleum Engineers, pp. 14–16. https://doi.org/10.2118/144391-MS. Dong, M., Dullien, F.A.L., Zhou, J., 1998. Characterization of waterflood saturation profile histories by the ‘Complete’ Capillary number. Transp. Porous Media 31, 213–237. https://doi.org/10.1023/A:1006565621860. Dullien, F.A.L., 2002. Porous Media, AcademicPress, SanDiego. Springer Berlin Heidelberg, Berlin, Heidelberg. https://doi.org/10.1007/978-3-662-04999-0. Fathi, E., Tinni, A., Akkutlu, I.Y., 2012. Correction to Klinkenberg slip theory for gas flow in nano-capillaries. Int. J. Coal Geol. 103, 51–59. https://doi.org/10.1016/J.COAL. 2012.06.008. Javadpour, F., Fisher, D., Unsworth, M., 2007. Nanoscale gas flow in shale gas sediments. J. Can. Pet. Technol. 46, 55. https://doi.org/10.2118/07-10-06. Katz, A.J., Thompson, A.H., 1986. Quantitative prediction of permeability in porous rock. Phys. Rev. B 34, 8179–8181. https://doi.org/10.1103/PhysRevB.34.8179. Labani, M.M., Rezaee, R., Saeedi, A., Al Hinai, A., 2013. Evaluation of pore size spectrum of gas shale reservoirs using low pressure nitrogen adsorption, gas expansion and mercury porosimetry: a case study from the Perth and Canning Basins, Western Australia. J. Pet. Sci. Eng. 112, 7–16. https://doi.org/10.1016/j.petrol.2013.11.022. Lan, Q., Dehghanpour, H., Wood, J., Sanei, H., 2015. Wettability of the montney tight gas formation. SPE Reserv. Eval. Eng. 18, 417–431. https://doi.org/10.2118/171620-PA. Larsen, J.W., 2004. The effects of dissolved CO2 on coal structure and properties. Int. J. Coal Geol. 57, 63–70. https://doi.org/10.1016/j.coal.2003.08.001. Li, S., Dong, M., Luo, P., 2017. A crossflow model for an interacting capillary bundle: development and application for waterflooding in tight oil reservoirs. Chem. Eng. Sci. 164, 133–147. https://doi.org/10.1016/J.CES.2017.01.059. Lopez, B., Aguilera, R., 2018. Petrophysical quantification of multiple porosities in shalepetroleum reservoirs with the use of modified pickett plots. SPE Reserv. Eval. Eng. 21, 187–201. https://doi.org/10.2118/171638-PA. McGlade, C., Speirs, J., Sorrell, S., 2013. Unconventional gas – a review of regional and global resource estimates. Energy 55, 571–584. https://doi.org/10.1016/j.energy. 2013.01.048. Nelson, P.H., 2009. Pore-throat sizes in sandstones, tight sandstones, and shales. Am. Assoc. Pet. Geol. Bull. 93, 329–340. https://doi.org/10.1306/10240808059. Pang, Y., He, Y., Chen, S., 2019. An innovative method to characterize sorption-induced kerogen swelling in organic-rich shales. Fuel 254, 115629. https://doi.org/10.1016/ j.fuel.2019.115629. Passey, Q.R., Bohacs, K., Esch, W.L., Klimentidis, R., Sinha, S., 2010. From oil-prone source rock to gas-producing shale reservoir - geologic and petrophysical characterization of unconventional shale gas reservoirs, in: international oil and gas conference and exhibition in China. Soc. Petrol. Eng. 740–741. https://doi.org/10. 2118/131350-MS. Riewchotisakul, S., Akkutlu, I.Y., 2016. Adsorption-enhanced transport of hydrocarbons in organic nanopores. SPE J. 21, 1960–1969. https://doi.org/10.2118/175107-PA. Sang, Q., Zhang, S., Li, Y., Dong, M., Bryant, S., 2018. Determination of organic and inorganic hydrocarbon saturations and effective porosities in shale using vacuumimbibition method. Int. J. Coal Geol. 200, 123–134. https://doi.org/10.1016/j.coal. 2018.10.010. Shi, Y., Yassin, M.R., Dehghanpour, H., 2018. A modified model for spontaneous imbibition of wetting phase into fractal porous media. Colloids Surf. A Physicochem. Eng. Asp. 543, 64–75. https://doi.org/10.1016/j.colsurfa.2017.12.052. Shi, Y., Yassin, M.R., Yuan, L., Dehghanpour, H., 2019. Modelling imbibition data for determining size distribution of organic and inorganic pores in unconventional rocks. Int. J. Coal Geol. 201, 26–43. https://doi.org/10.1016/j.coal.2018.11.010. Sigal, R.F., 2015. Pore-Size Distributions for Organic-Shale-Reservoir Rocks from NuclearMagnetic-Resonance Spectra Combined With Adsorption Measurements. https://doi. org/10.2118/174546-PA. Wang, Yu, Wang, L., Wang, J., Jiang, Z., Jin, C., Wang, Yanfei, 2018. Characterization of organic matter pores in typical marine and terrestrial shales, China. J. Nat. Gas Sci. Eng. 49, 56–65. https://doi.org/10.1016/j.jngse.2017.11.002. Wasaki, A., Akkutlu, I.Y., 2015. Permeability of organic-rich shale. SPE J. 20, 1384–1396. https://doi.org/10.2118/170830-PA. Yang, Z., Wang, W., Dong, M., Wang, J., Li, Y., Gong, H., Sang, Q., 2016. A model of dynamic adsorption–diffusion for modeling gas transport and storage in shale. Fuel 173, 115–128. https://doi.org/10.1016/j.fuel.2016.01.037. Yassin, M.R., Dehghanpour, H., Wood, J., Lan, Q., 2016. A Theory for relative permeability of unconventional rocks with dual-wettability pore network. SPE J. 21, 1970–1980. https://doi.org/10.2118/178549-PA. Yuan, W., Pan, Z., Li, X., Yang, Y., Zhao, C., Connell, L.D., Li, S., He, J., 2014. Experimental study and modelling of methane adsorption and diffusion in shale. Fuel 117, 509–519. https://doi.org/10.1016/J.FUEL.2013.09.046.

Greek symbols ∆P ∆t ∆V ∆x θ λ λn μ ρkerogen ρrock σ τ ∅inorg ∅org

Pressure difference Incremental time Total imbibed volume in organic pores within unit time step Incremental distance Contact angle Constant introduced when solving C∗ The nth positive roots of λ Viscosity of the fluid Density of kerogen Density of rock Interfacial tension Tortuosity of a rock sample Inorganic porosity Organic porosity

Acknowledgements The authors acknowledge the financial support from Discovery Grants of Natural Sciences and Engineering Research Council (NSERC) of Canada and the National Natural Science Foundation of China (No. 51774310). References Barnett, M., Units, D., Chalmers, G.R., Bustin, R.M., Power, I.M., Sciences, O., 2012. Characterization of gas shale pore systems by porosimetry, pycnometry, surface area, and field emission scanning electron microscopy/transmission electron microscopy image analyses: examples from the Barnett, Woodford, Haynesville, Marcellus, and Doig units. Am. Assoc. Pet. Geol. Bull. 96, 1099. https://doi.org/10.1306/ 10171111052. Bartley, J.T., Ruth, D.W., 1999. Relative permeability analysis of tube bundle models. Transp. Porous Media 36, 161–188. https://doi.org/10.1023/A:1006575231732. Cai, J., Perfect, E., Cheng, C.-L., Hu, X., 2014. Generalized modeling of spontaneous imbibition based on hagen–poiseuille flow in tortuous capillaries with variably shaped apertures. Langmuir 30, 5142–5151. https://doi.org/10.1021/la5007204. Clarkson, C.R., Solano, N., Bustin, R.M., Bustin, A.M.M., Chalmers, G.R.L., He, L., Melnichenko, Y.B., Radliński, A.P., Blach, T.P., 2013. Pore structure characterization of North American shale gas reservoirs using USANS/SANS, gas adsorption, and mercury intrusion. Fuel 103, 606–616. https://doi.org/10.1016/j.fuel.2012.06.119. Curtis, M.E., Ambrose, R.J., Sondergeld, C.H., Rai, C.S., 2011. Transmission and scanning electron microscopy investigation of pore connectivity of gas shales on the nanoscale. In: North American Unconventional Gas Conference and Exhibition. Society of

13