Process-based and dynamic 2D modeling of shale samples: Considering the geology and pore-system evolution

International Journal of Coal Geology 218 (2020) 103368

Contents lists available at ScienceDirect

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

Process-based and dynamic 2D modeling of shale samples: Considering the geology and pore-system evolution ⁎

T

⁎

Yuqi Wua,b, Pejman Tahmasebib, , Chengyan Lina, , Chunmei Donga a b

School of Geosciences, China University of Petroleum (East China), Qingdao 266580, China Department of Petroleum Engineering, University of Wyoming, Laramie 82071, USA.

A R T I C LE I N FO

A B S T R A C T

Keywords: Shale Pore system Organic matter Dynamic modeling Shale permeability

Accurate characterization of shale samples is of great importance for predicting and evaluating their intrinsic properties as they undergo various geological processes. As such, constructing shale models not only entails accounting for the numerous components and pores, but also requires considering the effects of thermal maturation on the evolutions of the pore systems and components. In this study, a new dynamic and process-based modeling for considering the geological evolutions occurred in pore systems is developed, which makes a significant improvement in reproducing the complex geological features in shale samples. Three types of pores, including organic matter (OM), interparticle (interP), and intraparticle (intraP) pores, and several other components including kerogen, pyrite, clay minerals, globigerinids, calcite, feldspar, quartz, and petroleum products, are all considered in our process-based modeling. In order to demonstrate the performance of our process-based method, two-dimensional (2D) dynamic shale models of multiple cases with different geological processes are constructed. Then, to characterize the constructed samples, the occurrences and distributions of all components are studied using a two-point correlation function. Aside from statistical comparisons, the gas flow in these shale models is also simulated to exhibit the permeability variations produced in the process-based models. Moreover, the degree of cementation, the fractions of interP, intraP, and OM pores are varied to mimic new geological scenarios, which aims to further evaluate the performance of the proposed technique. The results indicate that our process-based model presents an excellent performance. Furthermore, the presented modeling technique can be extended to investigating the effects of various components and pores on the intrinsic properties of shales in the future.

1. Introduction Thanks to a substantial amount of natural gas under the earth, shale gas has emerged as one of the most promising current and future energy resources (Tahmasebi, 2018a). As such, tremendous attention has been put on the characterization of shale reservoirs (Daigle et al., 2017; Hu et al., 2017; Javadpour, 2009; Kelly et al., 2016; Li et al., 2015; Loucks et al., 2012, 2009; Ma et al., 2014; Yang et al., 2016; Yu et al., 2017). Advanced imaging techniques, for instance, high-resolution scanning electron microscope (SEM), back-scattered scanning electron microscope (BSEM), transmission electron microscopy (TEM), field-emission scanning electron microscopy (FE-SEM), focused ion beam scanning electron microscope (FIB-SEM), and X-ray computed tomography (CT) have been used extensively to better characterize the features of mineralogy and pore systems of unconventional reservoirs (Karacan, 2003; Loucks et al., 2012; Milliken et al., 2013; Munawar et al., 2018;

⁎

Tahmasebi et al., 2016; Wu et al., 2019c, 2019b, 2019a). The characteristics of pore systems of unconventional reservoirs are also analyzed based on the indirect experiments, such as the mercury injection capillary pressure (MICP)(Ghanbarian et al., 2019), nuclear magnetic resonance (NMR) (Daigle et al., 2014), gas (i.e. N2 and CO2) adsorption (Karacan and Okandan, 2000; Okolo et al., 2015), and small-angle neutron scattering (SANS) (Sun et al., 2018; Wu et al., 2019c). On the other hand, various permeability measurement methods such as transient pulse decay (Cui et al., 2009), Gas Research Institute (GRI) method (Heller et al., 2014), and canister desorption test are proposed to estimate the permeability of shales. Apart from the experimental techniques, the computational methods have been widely used to generate reproduce the geological features observed in the shale samples (i.e., digital shale models) (Chen et al., 2015; Gerke et al., 2015; Houben et al., 2014; Mehmani and Prodanović, 2014; Naraghi and Javadpour, 2015; Naraghi et al., 2018;

Corresponding authors. E-mail addresses: [email protected] (P. Tahmasebi), [email protected] (C. Lin).

https://doi.org/10.1016/j.coal.2019.103368 Received 10 June 2019; Received in revised form 9 December 2019; Accepted 10 December 2019 Available online 16 December 2019 0166-5162/ © 2019 Elsevier B.V. All rights reserved.

International Journal of Coal Geology 218 (2020) 103368

Y. Wu, et al.

Table 1 History of pore-scale stochastic modeling of shales. Shale modela

Considered components

Pore types

References

Static

Not divided iOM; OM;

InterP pores; OM pores

Static Static

Pyrite; Calcite; Clay minerals; OM iOM; OM

Static Static Static

Not divided iOM; OM; Not subdivided iOM; OM; iOM; OM

Static Static Static Static Dynamic Static Dynamic

Not divided iOM; OM; Not divided iOM; OM Pyrite; Slits; Other iOM; OM Clay minerals; Feldspar; Quartz; Calcite; OM Clay minerals; OM; Other iOM Pyrite; Calcite; Slits; Dolomite; Ankerite; OM Pyrite; Clay minerals; Chemically stable minerals (e.g., silica); Chemically unstable minerals (e.g., calcite; dolomite; feldspar); Crystals (e.g., calcite; dolomite); Globigerinids (i.e., fossil bodies); OM (kerogen and petroleum products)

InterP InterP pores; InterP InterP InterP pores; InterP InterP InterP InterP InterP InterP InterP pores;

Naraghi and Javadpour (2015) Chen et al. (2015) Tahmasebi et al. (2016, 2015a) Yang et al. (2015) Gerke et al. (2015) Tahmasebi (2018a)

a

pores; OM pores pores; IntraP OM pores pores; OM pores pores; OM pores pores; IntraP OM pores pores; OM pores pores; OM pores pores; OM pores pores; OM pores pores; OM pores pores; OM pores pores; IntraP OM pores

Song et al. (2018a, 2018b) Cao et al. (2018) Li et al. (2018) Naraghi et al. (2018) Yao et al. (2018) Ji et al. (2019) This study

Dynamic shale models represent the generated models that can evolve with the change of the situations, for instance, thermal maturation.

ensemble-based model and used it to estimate the apparent permeability in shale samples (Naraghi and Javadpour, 2015). In their models, the iOM was not divided into various minerals. In addition, they did not take the intraP pores into consideration and the shape of OM is limited to square objects. After that, Chen et al. (2015) proposed a method, elementary building block algorithm, to generate the shale models (Chen et al., 2015), in which the presence of pyrite, calcite, and clay minerals are considered. However, the intraP pores are not taken into account. Meanwhile, Tahmasebi and his co-authors generated the stochastic shale models based on the widely available high-resolution 2D SEM images using the multiple-point geostatistical algorithm coupled with cross-correlation functions. In their method, the SEM images are used directly and without any parameter estimation/extraction. Despite the fact that the generated models in their studies are extremely similar to the realistic shales, they did not take the dynamic variation of shale models into account (Tahmasebi et al., 2015a, 2016). Yang et al. (2015) put forward a superposition algorithm for combining the stochastic models with inorganic and organic pores together. In their study, multiple-point statistics and Markov chain Monte Carlo methods are used to construct the inorganic and organic pores, respectively. However, intraP pores are not added into their models and the iOM is not carefully separated into different components (Yang et al., 2015). Similar issues are also encountered in other studies (Cao et al., 2018; Gerke et al., 2015; Song et al., 2018b, 2018a). Therefore, Li et al. (2018), Naraghi et al. (2018), and Ji et al. (2019) divided iOM into different components and inserted them into the shale models (Ji et al., 2019; Li et al., 2018; Naraghi et al., 2018). Even so, intraP pores are not taken into consideration and their models are still static. Generally speaking, the stochastic models can be divided into two groups, static and dynamic models, according to whether the models can evolve by varying the thermal and pressure (Øren and Bakke, 2002; Tahmasebi, 2018b). The static models cannot evolve with the change of temperature and pressure. However, both the temperature and pressure (i.e., thermal maturation) make a significant effect on the physical and chemical properties of rocks (Curtis et al., 2012; Ko et al., 2018, 2016; Loucks and Reed, 2014; Pommer and Milliken, 2015). Such a factor is more important in the shale formations as they manifest a variety of minerals and organic matter, which can change during the time. Hence, accurate modeling of shales requires accounting for the effect of the thermal and pressure maturations on the components and pore systems. Yao et al. (2018) used process-based mothed to construct the dynamic pore networks considering clay aggregates, but they did not take intraP pores into account and their models are relatively simple compared to the existing SEM images of shale samples (Yao et al., 2018).

Song et al., 2018a; Tahmasebi, 2018a; Tahmasebi et al., 2015a; Yu et al., 2019; Zhao et al., 2018). Due to the intrinsic complexity manifested in these formations, stochastic methods are recognized as the most reliable way of producing the variability. Compared to the experimental techniques, the computational methods not only overcome some drawbacks of the experiments, but are also time- and moneysaving. For example, conducting a single experiment, e.g., SEM, MICP, gas adsorption, or SANS, is very difficult if a multiscale characterization (e.g. from nanometer to micrometer) is sought. As such, it is very laborious to investigate and study the effects of different components and pores on the shale properties using the experimental techniques. By contrast, these issued can be easily addressed using digital shale models and computational methods (Chen et al., 2015; Tahmasebi, 2018a; Wu et al., 2019b; Zhao et al., 2018). For example, the deep learning techniques have recently been developed to enhance the resolution of shale images (Kamrava et al., 2019a) and also predict the permeablity within the general scope of rock physics (Kamrava et al., 2019b). In terms of the digital stochastic shale models, numerous studies have been reported in the last five years, listed in Table 1. However, before reviewing the history of stochastic shale models, it is necessary to briefly introduce the components of shale since the pore systems are strongly controlled by various parameters such as the components (i.e. compositions), the temperature, and the pressure of the formations (Ko et al., 2018). Based on the SEM images, it has been shown (see Fig. 1) that the components of shale samples comprise inorganic and organic matters. For many shale factions, the inorganic matrix (iOM) mainly constitutes quartz, calcite, feldspar, clay, and pyrite (Curtis et al., 2012). Furthermore, the pore systems in shale are extremely complex as they are not only composed of multiscale and multiple-type pores (or fractures). The pores in shale systems are also associated with the matter components of shale, which usually consist of interparticle (interP) pores, intraparticle (intraP) pores (Fig. 1), and organic matter (OM) pores (Ko et al., 2018, 2017, 2016; Liu et al., 2017; Loucks et al., 2012; Slatt and O'Brien, 2011). In a similar classification, Loucks et al. (2012) classified the iOM into three groups based on the mechanical and chemical stability of minerals (Loucks et al., 2012). Among the three groups, for example, the relatively chemically stable minerals represent that they cannot be easily dissolved. These minerals include silica and pyrite. On the chemically unstable side, the minerals mostly contain carbonate, feldspar, phosphate, and others. The intraparticle pores are usually generated in these unstable minerals. To date, various stochastic methods for the characterization of shale samples are presented. Table 1 summarizes the available methods. Naraghi and Javadpour (2015) constructed a simple stochastics 2

International Journal of Coal Geology 218 (2020) 103368

Y. Wu, et al.

Fig. 1. Components and pores in realistic shales. (A) InterP pores between quartz (Qtz) and calcite grains with cement overgrowths. IntraP and OM pores are also present; 8427 ft. (2569 m), vitrinite reflectance (Ro) = 1.5%, Lower Cretaceous (lower Bexar Shale Member) Pearsall Formation, Maverick County, Texas. (B) IntraP pores within a pyrite framboid. 6900 ft. (2103 m), 0.9% Ro, Upper Cretaceous Austin Chalk, La Salle County, Texas. (C) Cleavage-sheet intraP pores within a clay particle; 11,209 ft. (3417 m), 1.3% Ro, Upper Jurassic Haynesville Formation, Harrison County, Texas. (D) Dissolution pores after dolomite crystals and fossils. InterP pores also exist; 15,934 ft. (4857 m), approximately 1.8% Ro, Lower Cretaceous (Pine Island Shale Member) Pearsall Formation, Bee County, Texas. (E) IntraP pores within phosphate grain; 8845 ft. (2696 m), vitrinite reflectance (Ro) = approximately 1.2%, Lower Cretaceous (lower Bexar Shale Member) Pearsall Formation, Maverick County, Texas. (F) Large OM particle with OM pores. 7625 ft. (2324 m), Ro approximately 1.6%, Mississippian Barnett Shale, Wise County, Texas (Loucks et al., 2012). Figs. (A)–(F) are from Loucks et al. (2012), which can be reprinted by permission of the AAPG whose permission is required for further use. (G) InterP pore space was filled with solid bitumen or pyrobitumen. Lower Bexar Member of Lower Cretaceous Pearsall Formation, LaSalle County, Texas, Tidewater Oil Mable Wilson No. 2, 11,755 ft., calculated Ro = ~1.2%. (H) multichambered globigerinid with early calcite crystals followed by organic-matter emplacement. The middle chamber exhibits well-developed calcite crystals that precipitated before the formation of solid bitumen. Solid bitumen contains devolatilization cracks. Upper Cretaceous Eagle Ford Formation, Zavala County, Texas, Gose Hassett No. 3, 6214 ft., calculated Ro = ~0.8%. (I) Close-up of the bottom left chamber shown in H. Solid bitumen in intraP pores after calcite crystals grow, and solid bitumen exhibits a spongy OM pore texture network. (J) Solid bitumen or pyrobitumen displays a spongy OM-pore network. Mississippian Barnett Shale, Wise County, Texas, Texas United Blakely No.1, 7111 ft., calculated Ro = ~1.9–2.2% (Loucks and Reed, 2014). Figs. (G)–(J) are from Loucks and Reed (2014), which can be republished by permission of the Gulf Coast Association of Geological Societies, whose permission is required for further publication use.

3

International Journal of Coal Geology 218 (2020) 103368

Y. Wu, et al.

cementation takes place. This step is displayed in Fig. 2B. The cements can form the pores in interP and intraP, but might have no contribution for those places where kerogen has occluded interP pores (Loucks and Reed, 2014). Besides, the generated crystals also appear in intraP pores of globigerinid chambers. After these steps, the shale materials are more deeply buried, corresponding to a Ro of ~0.75% (Fig. 2C). At this time, the shales are in the oil window, which results in converting the kerogen to petroleum products such as bitumen and oil. These petroleum products are expelled from the kerogen and flow into the interP and intraP pores (Lewan, 1991; Loucks and Reed, 2014; Milliken et al., 2014; Pommer and Milliken, 2015). It should be noted that many geologists believe that OM can remain in place and migrate. Original depositional OM is composed of kerogen, which can evolve with the increase of thermal maturation from depositional kerogen into bitumen, solid bitumen, and finally into pyrobitumen (Lewan, 1991; Loucks and Reed, 2014; Milliken et al., 2014; Pommer and Milliken, 2015). Also, the migrated petroleum products contain OM pores, demonstrated in Fig. 2C. In addition, the organic acids are released during this period, which results in the dissolution of chemically unstable minerals (Loucks et al., 2012). When the temperature increases with the burial depth (Ro: about 1.3%, Fig. 2D), the gas is generated from the depositional kerogen and will create more OM pores (Ko et al., 2018, 2017, 2016). The oil and bitumen in the interP and intraP pores will transform into new petroleum products, for instance, solid bitumen and pyrobitumen. During the period, more OM pores will be generated in the new petroleum products (Ko et al., 2018, 2016; Loucks and Reed, 2014). Note that OM pores can develop in depositional kerogen and in the generated petroleum products. More information about the evolution of components and pore systems can be found in Loucks and other co-authors (Ko et al., 2018, 2017, 2016; Loucks and Reed, 2014).

Although great quantities of shale digital images/models have been constructed, there are still some drawbacks in the previous models. One of the obvious disadvantages is that, firstly, most of the previous models do not account for the development of components and pore systems in shales with the variation of thermal and pressure maturations. Secondly, many of the previous studies have not carefully divided the iOM into specific minerals. In other words, the iOM is not fully and realistically considered. Due to this drawback, the intraP pores have usually been ignored. Furthermore, the current methods are not able to fully reproduce the processes that are taking place in shale formations, but only the final models are generated. Producing a sequence of models, starting from the initial depositional time to when the shale formations are constructed, can be very useful. In this study, the abovementioned issues are addressed by developing a more realistic and process-based modeling technique. In this method, various components and pores can be considered and reproduced in the models. Besides, the evolutions of components and pore systems with the change of thermal and pressures maturations are also considered. Aside from all these new features, the effects of cementation degree, dissolution degree, and the percentage of petroleum products and OM pores are also taken into consideration. The rest of the paper is organized as follows. Section 2 gives a detailed description of pore-system evolution with varied thermal maturation, which is the basic theory of our dynamic and process-based modeling. Next, our dynamic modeling algorithm and workflow based on quartet structure generation set (QSGS) and image processing operation are introduced. Section 3 presents the generated shale models. In this section, the generated modes will undergo an extensive sensitivity analysis by varying pore distributions, the fractions of the utilized elements, and gas flow in these models are simulated. Moreover, the degree of cementation, the fractions of interP, intraP, and OM pores are changed in the algorithm to generate new dynamic models. Finally, Section 4 summaries the advantages and limitations of the proposed algorithm.

2.2. Dynamic pore-scale modeling of shales Before describing the dynamic modeling algorithm, two algorithms, quartet structure generation set (QSGS) and dilation algorithms, are first introduced. The QSGS will be used to generate the pyrite, crystals, kerogen, petroleum products, OM pores, intrap pores, and chemically stable and unstable minerals. The dilation operation is performed to produce the cements in shales.

2. Methodology 2.1. Pore-system evolution Before describing the processes involved in pore-system evaluation in shale systems, some geological terms used in this paper are briefly introduced (Loucks and Reed, 2014). OM represents any liquid or solid composed of organic material, including kerogen and petroleum products. Kerogen is an insoluble organic material that can be converted to petroleum products through thermal maturation. Petroleum products, themselves, contain various hydrocarbons such as bitumen, solid bitumen, pyrobitumen, oil, and gas. With the aid of the experimental techniques, it has been shown that the existing components and pores in shale systems are greatly affected by the thermal maturation (i.e., temperature and pressure). Loucks and Reed (2014) summarized an idealized history of pore-system evolution with increasing thermal maturity, shown in Fig. 2. Fig. 2A displays the original and nonlithified sediment including various minerals and kerogen. When the pressure of the formation increases, the sediment will be compacted accordingly; Fig. 2(B–D). During this process, interP pore space greatly decreases and the kerogen is squeezed. In this study, we mainly study the pore-system development when Ro changes from 0.4% to 1.3%. It should be noted that Ro represents vitrinite reflectance. Vitrinite reflectance is a measure of the percentage of incident light reflected from the surface of vitrinite particles in sedimentary rocks. This measure is often given as a mean value to represent all vitrinite particles measured in an individual sample (Dembicki, 2009). During the shallow burial period (Ro: about 0.4%, see Fig. 2A), there are no OM pores due to too low thermal maturity. A large number of interP and a few of intraP pores (in globigerinid chambers) exist in shale systems. When the temperature and pressure increase, the

2.2.1. Quartet structure generation set The methods of modeling stochastic porous media have been widely used to characterize the porous materials and predict the their properties in hydrology, petroleum engineering, soil science, biology, fuel cells, and construction engineering (Blunt et al., 2013; Chengyan et al., 2018; Karsanina et al., 2018; Tahmasebi, 2018c; Tahmasebi and Sahimi, 2015, 2013, 2012; Wu et al., 2018a, 2018b; Zhu et al., 2019). The detailed literature review of modeling stochastic porous media can be found in Li et al. (2018), Wu et al. (2018a), and Tahmasebi and Sahimi (2012) (Chengyan et al., 2018; Tahmasebi and Sahimi, 2012; Wu et al., 2018a). In this study, we employed the QSGS method to generate various minerals and pores in shale samples. The QSGS algorithm was developed based on the grow-with-time model (Coveney et al., 1998) and cluster growth theories (Wang et al., 2007). Modeling a stochastic shale sample using QSGS requires four steps that are described as follows, which are conducted on a simulation grid G: (1) Classifying: The minerals and pores should be classified into several groups or phases according to their similarity of characteristics. For example, the calcite, dolomite, and feldspar belong to a phase since they are chemically unstable and easily dissolved by acid fluid (Ma et al., 2018; Yuan et al., 2019). For this modeling, the phases contain some distinct elements such as pyrite, crystals, petroleum products, kerogen, intrap pores, and chemically stable and unstable minerals. (2) Seeding: We randomly select several cells in G as the seeds of the 4

International Journal of Coal Geology 218 (2020) 103368

Y. Wu, et al.

Fig. 2. Idealized history of pore-system development with increasing thermal maturation (Loucks and Reed, 2014). In this study, we mainly consider the pore-system evolution of shales in the last four stages (A, B, C, and D).

5

International Journal of Coal Geology 218 (2020) 103368

Y. Wu, et al.

C, and D in Fig. 4 are based on stages A, B, C, and D in Fig. 2, respectively.

(3) Growth: The particles of every phase will grow from the seeds to adjacent cells along with the eight directions (Fig. 3). The growth probability Gp of each route can govern the growth directions and the shape of the particles. Thus, the produced particles can be circles, ovals, squares, diamonds, and rectangles, which all are typical in shale formations. (4) Repeating: Repeat step (3) until fi (i.e. the volume fraction of phase i) reaches the predefined value and continue with steps (2) and (3) for all remaining phases.

(1) Model A: This model corresponds to stage A in Fig. 2. At this stage, Ro is about 0.4% and the shale samples are in the pre-oil window. First, the simulation grid G is filled with interP pores. Then, the globigerinids are created using several adjacent rings. The clay minerals are similar to intersected curves. Next, the pyrite, chemically stable and unstable minerals, and depositional kerogen are produced using the QSGS algorithm. Then, the intraP pores in pyrite are generated. The fractions of these components are taken into consideration while using QSGS to generate them. After performing the above-described steps, the nonlithified sediment is generated. Subsequently, the seam carving algorithm is executed to remove a large amount of interP pore space from the nonlithified sediment. The seam carving algorithm can change the size of an image and keep the high-energy objects unchanged through removing low-energy pixels (Avidan and Shamir, 2007). InterP pores in nonlithified muddy sediment are low-energy objects. Thus, the compaction and removal of interP pores can be simulated using seam carving algorithm. This algorithm is used to mimic the gravity force and compaction in sedimentary environments. As another alternative, one can use the discrete mechanical modeling, such as a discrete element method, to perform the compaction (Tahmasebi and Kamrava, 2019; Zhang and Tahmasebi, 2019, 2018). Finally, the morphological alteration of kerogen (e.g. squeeze and expansion) can be simulated using the dilation operation. (2) Model B: From stage A to stage B, the main diagenesis is cementation. The chemically stable and unstable minerals and globigerinids overgrow along their boundaries. The overgrowth is considered using dilation operation. It should be noted that the cements cannot be formed in the locations where kerogen has occluded interP pores. Besides, several crystals are also generated in the globigerinid chambers (Fig. 4). (3) Model C: For stage C in Fig. 2, the samples are in the oil window. During this stage, the petroleum products are produced and migrate into the interP and intraP pores. OM pores exist in both deposition kerogen and migrated petroleum products (Fig. 4). Meanwhile, the intraP pores are produced into the chemically unstable minerals such as calcite and feldspar. These processes can be accomplished with the aid of the QSGS algorithm. (4) Model D: The temperature and pressure become higher for this model. The Ro value reaches 1.3% or larger. The changes lead to generating more OM pores in kerogen and petroleum products. Also, more OM pores can be produced using the QSGS method.

2.2.2. Dilation operation In mathematical morphology, the dilation is one of the commonly used operations which can probe and expand the included shapes in an image. The dilation fills the gulfs at the boundaries of the object using a structuring element, which enlarges the size of the object. One object A in an image is dilated by a structuring element B (Gonzalez et al., 2004), which is defined by:

In our study, the order of placing the components and pores are critically taken into account based on the demonstrated sequence in Fig. 2. As above mentioned, pore particle radii can be considered using the QSGS algorithm. Herein, some pre-defined size distributions for each component and pore are used. For example, OM pores are smaller than InterP pores in this algorithm. A set of dynamic shale models can be generated by following the above-mentioned steps.

Fig. 3. Eight growth directions for every cell. The red cell represents the seed of a particle. The seed can grow along with the green arrows. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

particles, which occurs when the particles in the sample are generated. The seed number can be controlled by a seed distribution probability (Sd). For phase i, Sdi can not only guarantee the statistically uniform distributions of the seeds in G (Chen et al., 2015), but also affect the average particle radius ri as follows:

Sdi =

fi Kri2

,

(1)

where fi is the volume fraction of phase i, and K is the coefficient related to the shape of the particle. K is π when the shape of the particles tends to be circles, and K is 1 when the particles are squares. Therefore, the particle and pore sizes are also taken into consideration in this modeling algorithm.

A⨁B = {z ∈ E | (B s )z ⋂A ≠ ∅}

(2) 2.3. Calculation of pore size

where B denotes the symmetric B, B = {x ∈ E| −x ∈ B}. The dilation operation has been used to generate the cementation process in rocks (Tahmasebi and Kamrava, 2018). Similarly, we will employ it to simulate the cementation in shale samples. s

s

As mentioned earlier, the pores are distributed in both OM and iOM. In the previous studies, the pores are usually considered as circles through which the equivalent pore radius r is calculated:

2.2.3. Dynamic modeling In this section, the pore-system development from stage A to stage D in Fig. 2 is studied. Based on the analysis of pore-system evolution with the variation of thermal maturation in Section 2.1, the dynamic shale models are constructed in four stages. The detailed flowchart of dynamic pore-scale modeling is demonstrated in Fig. 4. The models A, B,

re =

A , π

(3)

where A is the area of a pore. In reality, the shape of many pores in 2D images is not circular. Thus, a new technique based on the pore shape factor is presented to compute the equivalent pore size in 2D images 6

International Journal of Coal Geology 218 (2020) 103368

Y. Wu, et al.

Fig. 4. Main flowchart of dynamic pore-scale modeling of shales.

where μ is the gas viscosity, M is the molar mass, ϕ is the porosity of porous medium, R is the gas constant, T is the temperature, ρav is the average gas density, τ is the tortuosity of pore space, Df is the fractal dimension of the pore boundary, δ is the ratio of normalized molecular size to the local average pore diameter, and kD is the Darcy permeability (Akkutlu et al., 2015; Singh and Javadpour, 2014). Besides, b and Dk can be calculated by:

(Wu et al., 2019c). This technique can generate a more reasonable pore size distribution for 2D images. The detailed steps of this technique are as follows: (i) Separate the pore space into individual pores and obtain the area (A) and perimeter (P) of each pore. (ii) Calculate the shape factor Sf using (Mason and Morrow, 1991):

Sf =

A . P2

(4)

Note that the shape factor is a widely used parameter for describing the pore morphology. The larger the shape factor of a pore, the smoother the pore is and the closer the pore is to circle. As to the triangle, square, and circle, their shape factors are 0.0481, 0.071, and 0.0796, respectively. (iii) To accurately get the equivalent radius of a pore, the pore is equated with one circular, square, or triangular pore according to Sf. Then, the corresponding radius is equal to the radius of the inscribed circle of the equivalent geometry. If Sf < 0.0481, the pore is considered as a triangular object, and the radius of the pore can be obtained by re = 2A/P. If Sf > 0.071, the pore will be taken as a circle, and its radius is calculated by re = A/ π . Apart from these two cases, the radius of the pore is set as one square and its equivalent radius can be given by re = A .

8πRT 0.5 μ 2 ⎛ − 1⎞ ⎞ b=⎛ ⎝ M ⎠ Rav ⎝ α ⎠

(6)

2Rav 8RT 0.5 ⎛ ⎞ , 3 ⎝ πM ⎠

(7)

Dk =

where Rav is the weighted arithmetic average pore radius. Ghanbarian and Javadpour (2017) proposed a more effective method to evaluate Rav via the effective approximation (Ghanbarian and Javadpour, 2017). The tangential momentum accommodation coefficient α can be obtained by (Agrawal and Prabhu, 2008):

α = 1 − log(1 + Kn0.7).

It should be noted that α is negative when Knudsen number Kn is larger than 23.1(Ghanbarian, 2018). Kn can be calculated by:

2.4. Simulation of gas flow

Kn = The simulation of the gas flow in shales is a complex process since the gas flow in nanometer pores does not abide by the conventional Darcy's law. Furthermore, gas flow mechanisms in OM and iOM pores are not the same (Javadpour, 2009; Song et al., 2018b). In this study, different control equations are used to simulate gas flow in iOM and OM pore spaces. The apparent permeability function (Naraghi and Javadpour, 2015; Singh and Javadpour, 2016) and Langmuir sorption permeability (Darabi et al., 2012; Naraghi et al., 2018) are employed to calculate the permeability of gas flow in iOM and OM pores, respectively. The apparent permeability accounting for slip flow and Knudsen diffusion was originally defined by Javadpour (2009) and then modified as (Darabi et al., 2012):

K app =

μM ϕ b (δ ) Df − 2Dk + kD ⎜⎛1 + ⎟⎞, RTρav τ p⎠ ⎝

(8)

λ KB T = , Rav Rav 2 πσ 2P

(9)

where KB is the Boltzmann constant in J/K, σ is the collision diameter of molecules, and P is the pressure of the system (Naraghi and Javadpour, 2015). Thus, one can use Eq. (5) to get the apparent permeability. Besides, the average permeability of the shale models can be obtained using the geometric average method (Naraghi and Javadpour, 2015), which can be expressed by:

KLB =

KUB =

K app =

(5) 7

Nc Nr

Nr

∑

1 N

1 ki, j

(10)

Nc 1 Nr ∑Nc N 1 j=1 r

∑i = 1 ki, j

(11)

KUB KLB ,

(12)

i=1

∑ j =c 1

International Journal of Coal Geology 218 (2020) 103368

Y. Wu, et al.

Fig. 5. The generated dynamic shale models using the dynamic modeling algorithm. These models are generated based on the analysis of the processes in Fig. 2. The sizes of the first image and other images are 168 and 144 μm2, respectively.

8

International Journal of Coal Geology 218 (2020) 103368

Y. Wu, et al.

Fig. 6. 16 generated models of four cases using the dynamic modeling algorithm. These models correspond to the schematic diagrams in Fig. 1. The sizes of the models for cases 1 and 2 are 144 μm2, and the sizes of the models of cases 3 and 4 are 196 μm2.

where Ds is the surface diffusion coefficient, Γmax is the maximum Langmuir adsorption, ρgrain is the grain density of the sample, Bg is the gas formation volume factor, fTOC is the fraction of the OM, and β is the Langmuir constant (Naraghi et al., 2018). If the pores are much larger than 2 nm, the effect of surface diffusion is trivial (Wu et al., 2016). More details about calculating the apparent permeability based on 2D shale models can be found in Naraghi and Javadpour (2015) and Naraghi et al. (2018) (Naraghi and Javadpour, 2015; Naraghi et al., 2018).

where KLB and KLB are the permeability of lower and upper bounds, respectively, Nr and Nc are the number of grids in the vertical and horizontal directions, respectively, and ki, j is the permeability of the grids in row i and column j (Naraghi and Javadpour, 2015). On the other hand, the surface diffusion should be accounted for when gas flows in OM pores. The contribution of surface diffusion to the permeability can be expressed by (Naraghi et al., 2018):

Ks = μDs

Γmax ρgrain Bg (1 − ϕ) ϕ fTOC

1 τ β P+

(

1 2 β

)

, (13) 9

International Journal of Coal Geology 218 (2020) 103368

Y. Wu, et al.

Table 2 Components percentage and porosity in dynamic shale models. Case 1 Percentage (%) Pyrite Clay minerals Globigerinids Kerogen Chemically-unstable minerals Chemically-stable minerals InterP pores IntraP pores in pyrite IntraP pores in globigerinids Crystals in globigerinids Petroleum products generated in interP pores Petroleum products generated in globigerinids OM pores in kerogen IntraP pores in chemically unstable minerals OM pores in petroleum products

1A 1.65 2.44 1.38 16.46 28.53 28.59 17.36 1.42 2.17

Case 2 1B 1.65 2.44 1.54 16.46 29.42 31.25 13.65 1.42 1.92 0.25

1C 1.65 2.44 1.54 14.54 27.03 31.25 2.48 1.42

1D 1.65 2.44 1.54 13.20 27.03 31.25 2.48 1.42

0.25 9.29 1.35 1.92 2.39 2.45

0.25 8.85 1.35 3.26 2.39 2.89

2A 2.19 1.73 1.89 28.30 23.23 24.99 12.83 2.04 2.80

Case 3 2B 2.19 1.73 2.17 28.30 24.68 25.77 10.33 2.04 2.32 0.49

2C 2.19 1.73 2.17 25.09 22.85 25.77 1.96 2.04

2D 2.19 1.73 2.17 23.00 22.85 25.77 1.96 2.04

0.49 7.03 1.69 3.21 1.83 1.96

0.49 6.56 1.69 5.3 1.83 2.43

3A 0.95 2.31 1.11 32.92 22.14 27.4 10.98 0.88 1.30

Case 4 3B 0.95 2.31 1.22 32.92 23.37 28.23 8.80 0.88 1.24 0.06

3C 0.95 2.31 1.22 28.71 21.62 28.23 1.80 0.88

3D 0.95 2.31 1.22 25.95 21.62 28.23 1.80 0.88

0.06 5.94 0.86 4.21 1.75 1.44

0.06 5.69 0.86 6.97 1.75 1.69

4A 0.76 2.25 2.77 29.73 25.70 20.54 12.14 0.64 5.48

4B 0.76 2.25 3.14 29.73 26.62 22.03 9.35 0.64 4.51 0.97

4C 0.76 2.25 3.14 26.27 24.63 22.03 1.45 0.64

4D 0.76 2.25 3.14 24.05 24.63 22.03 1.45 0.64

0.97 6.80 3.28 3.46 2.00 2.33

0.97 5.30 3.28 5.68 2.00 2.83

Fig. 7. Size distributions of all pores (including iOM and OM pores) of 16 models from four cases.

3. Results and discussion

generated accordingly.

Following the flowchart of dynamic modeling of shales, 21 models are generated to demonstrate the performance of the proposed algorithm, shown in Figs. 5 and 6. In addition, the geometrical properties and two-point correlation functions of pore systems in 16 models from the designated four cases are characterized. Moreover, the gas flow in these models is simulated. Finally, we considered the percentage variations of components and pores in shales and various models are

3.1. Dynamic models To represent the taken steps for building a stochastic shale sample and also the correlation between our proposed method and the previously described geological processes, five dynamic shale models are displayed in Fig. 5. In other words, these five models are exactly the five stages in Fig. 2. It can be seen from Fig. 5 that various components 10

International Journal of Coal Geology 218 (2020) 103368

Y. Wu, et al.

Fig. 8. Size distributions of OM pores of models C and D from four cases.

pores when the thermal maturation increases. In the models C and D, the OM pores occupy a larger fraction than those of interP and intraP pores. The variations of components and pores also indicate that various shale models can be generated using the proposed algorithm.

including clay minerals, pyrite, globigerinids, chemically unstable and stable minerals, and kerogen are available in the nonlithified sediment. Their morphological features are extremely similar to the realistic characteristics. Model A can be obtained after the interP pores are compacted and the kerogen is squeezed, which represents the compaction stage. From model A to model B, the chemically unstable and stable minerals, and globigerinids overgrow. The differences in minerals before and after cementation are shown in the white circles in Fig. 5. The cementation causes a reduction in the porosity of interP pores. Besides, the crystals appear in the globigerinid chambers in model B, which corresponds to stage B in Fig. 2. Next, the petroleum products (light red) are generated from the kerogen and migrated to interP and intraP pores in model C, which results in a significant reduction of interP pores. Meanwhile, numerous OM and intraP pores are also produced in OM and chemically unstable minerals, respectively. In terms of model D in Fig. 5, more OM pores are generated in OM compared to model C. The other four different cases manifested in Fig. 6 are generated to exhibit the performance of our dynamic modeling algorithm. Herein, we only exhibit the shale models from stage A to stage D since we mainly investigate the variations of pore-system properties and permeability. The size of the models for cases 1 and 2 is 144 μm2, and the size of the models of cases 3 and 4 is 196 μm2. The fractions of the components and porosity of all produced models are computed and listed in Table 2. Four cases with the fractions of components represent different shales in the realistic formations. For instance, there are more OM in the models from cases 3 and 4, compared with the models from cases 1 and 2. Generally speaking, more OM may generate more OM

3.2. Characterization of pore systems The size distributions of all pores, including iOM and OM pores, in the produced models are calculated through which one can achieve a quantitative way to study the differences between the models. The size distributions of all pores in 16 models from four cases are displayed in Fig. 7. As can be observed from Fig. 7, the pore size distributions have multiple peaks, which is close to the pore sizes in the realistic shales (Tahmasebi et al., 2015b). In other words, two dominate distributions are recognized for the shale samples which represent the pore size distributions for OM and iOM. The pore size distribution in model A has a wider distribution, compared with the curves of other models of the same case. The pore sizes of models C and D become very small due to the fact that the generated petroleum products occupy most of interP pore spaces in shales. The pore size distributions of models C and D from the same case are similar, which indicates that only more OM pores are generated in model D. To make the comparisons clearer, the size distributions of OM pores of models C and D from four cases are also obtained and the results are demonstrated in Fig. 8. It is clear that the pore size distribution of model C is similar to that of model D for the same case. Compared the size distributions of all pores and OM pores in models C and D, it can be observed that the curves of OM pores are narrower than those of all pores, which reveals that OM pores are 11

International Journal of Coal Geology 218 (2020) 103368

Y. Wu, et al.

Fig. 9. Two-point correlation functions (TPCFs) of cases 1, 2, 3 and 4.

P (i, i + h ) = E {I (i) I (i + h )} = Prob {I (i) = 1, I (i + h ) = 1},

(14)

where h is the distance between two pixels in any direction. I(i) is an indicator function. For example, I(i) = 1 represents that i fall into the void space. P(i, i + h) can be interpreted as the probability of two pixels to be simultaneously in pore space (Krishnan and Journel, 2003; Tahmasebi and Sahimi, 2012). The TPCFs for 16 models from four cases are computed and the results are shown in Fig. 9. As can be seen, the larger the porosity, the larger the TPCF. It can be found that the TPCF of model A is larger than the other models, which is due to containing a large number of pore spaces. 3.3. Simulation of gas flow Since the pores in iOM and OM are inserted separately, thus, one can solve the gas flow in different pores. The pore size, porosity, fractal dimension, tortuosity and the fraction of OM can be obtained from the shale models. Nc and Nr are set as 500, temperature and pressure of simulation are set as 300 K and 2000 psi, respectively. The Maximum Langmuir adsorption capacity (Γmax) is set as 110 scf/ton, and the Langmuir adsorption coefficient (β) is set as 0.165 MPa−1. The permeability values of 16 shale models of four cases are calculated and presented in Fig. 10. As can be observed, the permeability from model A to model C decreases rapidly. This can be explained by the pore-system development from model A to model C. In fact, the gradual decrease in the porosity and pore size from model A to model C results in the drop of the permeability. To verify this, more OM pores are inserted into model D compared to model C, which results in a larger permeability in model D. This brings us to this conclusion that the permeability of the shale models depends on the pore size and porosity of the models. In

Fig. 10. Permeability variations of the sixteen shale models of four cases.

smaller than iOM pores in our shale models. It should be mentioned that the apertures of OM and iOM pores in our shale models can vary by using different parameters in the dynamic modeling algorithm. The statistical correlation function is also a useful parameter for characterizing the pore structures (Ding et al., 2018; Krishnan and Journel, 2003; Tahmasebi, 2018d; Tahmasebi and Sahimi, 2012; Wu et al., 2018a, 2019b, 2019a). For example, the two-point correlation function (TPCF) describes the probability P(i, i + h) of two pixels in the image with distance h, where both laid in a certain phase (pore or solid space). TPCF can be calculated by: 12

International Journal of Coal Geology 218 (2020) 103368

Y. Wu, et al.

Fig. 11. The models before and after enhancing cementation for case 1.

Fig. 12. The models before and after increasing interP pores for case 1.

Fig. 13. The models before and after increasing intraP pores for case 1.

3.4. Dynamic variation

addition, the permeablity varies with the size when the sample's dimension is smaller than the REV (Representative Elementary Volume). If the model is larger than the REV, the simulation of the permeability does not anymore rely on its size. In this study, the permeability of the models tends to be stable when the they are larger than 121 μm2.

Apart from the dynamic models shown in Fig. 6, more variations in the generated models in any stage can also be simulated using the proposed method. The dynamic variations represent that all the components (e.g., various minerals and OM) and pores (e.g., interP and intraP pores) in the models can be changed. As such, the cementation degree, interP pores, and intraP pores in shale models for four cases 13

International Journal of Coal Geology 218 (2020) 103368

Y. Wu, et al.

Fig. 14. The models before and after increasing the cementation, interP pores, and intraP pores for cases 2, 3, and 4. One can compare the areas in the white circles in the same locations of reference and new models to observe the differences.

(Figs. 11–14) are varied to demonstrate the performance of the presented method.

there are more interP pores in the second model than those in the first model.

3.4.1. Variation of cementation As discussed earlier, cementation can be simulated using the dilation morphologic operation. In the dilation operation, the size of the structuring element can affect the degree of cementation. Therefore, one can increase the size of this structuring element to enhance the cementation in shales. As can be observed from Fig. 11, the larger structuring element was used to enhance the cementation, which results in that the stronger cementation occurred in shales. The models before and after enhancing cementation are manifested in Fig. 11. One can observe the differences before and after enhancing cementation in the white circles. The cementation is stronger in the second model then the first model in Fig. 11.

3.4.3. Variation of intraP pores The dissolution is also one of the common diageneses in the sedimentary rocks. For different shales, the degree of the dissolution maybe not identical, which leads to variation of intraP pores. Herein, the variations of dissolution and intraP pores are also taken into account. Two models from case 1 with different interP pores are generated, exhibited in Fig. 13. As can be seen, there are more intraP pores in the second model than the first model. 3.4.4. Variation of OM pores OM pores can make a great contribution to the permeability of shales. The fractions of OM pores may greatly vary in different shale formations. As such, the OM pores should be easily changed according to different scenarios. Using the above-proposed algorithm, one can produce models with different percentages of OM pores. The variations can be seen from models C and D for four cases in Fig. 6. Moreover, the dynamic models with different cements, interP pores, and intraP pores of cases 2, 3, and 4 are generated and exhibited in Fig. 14. In addition, the proposed algorithm can generate an extremely huge shale model without any limitations in terms of computational time and memory. The dynamic variations of these components allow us to generate any shale models with various combinations and sizes. Furthermore, these dynamic models with different components can be used to investigate the effects of every component (e.g., OM and intraP

3.4.2. Variation of interP pores Under some specific circumstances of temperature and pressure, petroleum products with fewer fractions may be generated in shales during the oil window. In other words, the presence of fewer petroleum products indicates that there are more reserved interP pores. To simulate this, one can adjust the percentage of petroleum products or interP pores in shale models. In our algorithm, the fractions of petroleum products and interP pores can be adjusted to achieve this goal. Let us assume that fewer petroleum products and more interP pore must be generated. The variations of the models are simulated using the proposed algorithm and the results are displayed in Fig. 12. It is clear that 14

International Journal of Coal Geology 218 (2020) 103368

Y. Wu, et al.

pores) on shale properties.

Avidan, S., Shamir, A., 2007. Seam carving for content-aware image resizing. In: ACM SIGGRAPH 2007 Pap.- SIGGRAPH ‘07. 26. pp. 10. https://doi.org/10.1145/ 1275808.1276390. Blunt, M.J., Bijeljic, B., Dong, H., Gharbi, O., Iglauer, S., Mostaghimi, P., Paluszny, A., Pentland, C., 2013. Pore-scale imaging and modelling. Adv. Water Resour. 51, 197–216. https://doi.org/10.1016/j.advwatres.2012.03.003. Cao, G., Lin, M., Jiang, W., Zhao, W., Ji, L., Li, C., Lei, D., 2018. A statistical-coupled model for organic-rich shale gas transport. J. Pet. Sci. Eng. 169, 167–183. https://doi. org/10.1016/J.PETROL.2018.05.033. Chen, L., Kang, Q., Dai, Z., Viswanathan, H.S., Tao, W., 2015. Permeability prediction of shale matrix reconstructed using the elementary building block model. Fuel 160, 346–356. https://doi.org/10.1016/J.FUEL.2015.07.070. Chengyan, L., Yuqi, W., Lihua, R., Yang, W., Weichao, Y., Xiaolong, S., Xianguo, Z., Yimin, Z., 2018. Review of digital core modeling methods. Prog. Geophys. 33, 0679–0689. Coveney, P.V., Maillet, J.-B., Wilson, J.L., Fowler, P.W., Al-Mushadani, O., Boghosian, B.M., 1998. Lattice-gas simulations of ternary amphiphilic fluid flow in porous media. Int. J. Mod. Phys. C 9, 1479–1490. https://doi.org/10.1142/s0129183198001345. Cui, X., Bustin, A.M.M., Bustin, R.M., 2009. Measurements of gas permeability and diffusivity of tight reservoir rocks: different approaches and their applications. Geofluids 9, 208–223. https://doi.org/10.1111/j.1468-8123.2009.00244.x. Curtis, M.E., Cardott, B.J., Sondergeld, C.H., Rai, C.S., 2012. Development of organic porosity in the Woodford Shale with increasing thermal maturity. Int. J. Coal Geol. 103, 26–31. https://doi.org/10.1016/J.COAL.2012.08.004. Daigle, H., Johnson, A., Thomas, B., 2014. Determining fractal dimension from nuclear magnetic resonance data in rocks with internal magnetic field gradients. Geophysics. https://doi.org/10.1190/Geo2014-0325.1. Daigle, H., Hayman, N.W., Jiang, H., Tian, X., Jiang, C., 2017. Multiscale pore networks and their effect on deformation and transport property alteration associated with hydraulic fracturing. Energy Procedia 125, 71–79. https://doi.org/10.1016/J. EGYPRO.2017.08.067. Darabi, H., Ettehad, A., Javadpour, F., Sepehrnoori, K., 2012. Gas flow in ultra-tight shale strata. J. Fluid Mech. 710, 641–658. https://doi.org/10.1017/jfm.2012.424. Dembicki, H., 2009. Three common source rock evaluation errors made by geologists during prospect or play appraisals. Am. Assoc. Pet. Geol. Bull. 93, 341–356. https:// doi.org/10.1306/10230808076. Ding, K., Teng, Q., Wang, Z., He, X., Feng, J., 2018. Improved multipoint statistics method for reconstructing three-dimensional porous media from a two-dimensional image via porosity matching. Phys. Rev. E 97, 1–10. https://doi.org/10.1103/PhysRevE.97. 063304. Gerke, K.M., Karsanina, M.V., Mallants, D., 2015. Universal stochastic multiscale image fusion: An example application for shale rock. Sci. Rep. 5, 1–12. https://doi.org/10. 1038/srep15880. Ghanbarian, B., 2018. Estimating gas relative permeability of shales from pore size distribution. In: SPE Annu. Tech. Conf. Exhib. September, SPE-191878-MS, https://doi. org/10.2118/191878-ms. Ghanbarian, B., Javadpour, F., 2017. Upscaling pore pressure-dependent gas permeability in shales. J. Geophys. Res. Solid Earth 122, 2541–2552. https://doi.org/10.1002/ 2016JB013846. Ghanbarian, B., Torres-Verdín, C., Lake, L.W., Marder, M., 2019. Gas permeability in unconventional tight sandstones: Scaling up from pore to core. J. Pet. Sci. Eng. 173, 1163–1172. https://doi.org/10.1016/J.PETROL.2018.10.057. Gonzalez, R.C., Woods, R.E., Eddins, S.L., 2004. Digital Image Processing Using Matlab. Pearson Prentice Hallhttps://doi.org/10.1117/1.3115362. Heller, R., Vermylen, J., Zoback, M., 2014. Experimental investigation of matrix permeability of gas shales. Am. Assoc. Pet. Geol. Bull. 98, 975–995. https://doi.org/10. 1306/09231313023. Houben, M.E., Desbois, G., Urai, J.L., 2014. A comparative study of representative 2D microstructures in Shaly and Sandy facies of Opalinus Clay (Mont Terri, Switzerland) inferred form BIB-SEM and MIP methods. Mar. Pet. Geol. 49, 143–161. https://doi. org/10.1016/J.MARPETGEO.2013.10.009. Hu, Q., Zhang, Y., Meng, X., Li, Z., Xie, Z., Li, M., 2017. Characterization of micro-nano pore networks in shale oil reservoirs of Paleogene Shahejie Formation in Dongying Sag of Bohai Bay Basin, East China. Pet. Explor. Dev. 44, 720–730. https://doi.org/ 10.1016/S1876-3804(17)30083-6. Javadpour, F., 2009. Nanopores and apparent permeability of gas flow in mudrocks (shales and siltstone). J. Can. Pet. Technol. 48, 16–21. https://doi.org/10.2118/0908-16-DA. Ji, L., Lin, M., Cao, G., Jiang, W., 2019. A core-scale reconstructing method for shale. Sci. Rep. 9, 4364. https://doi.org/10.1038/s41598-019-39442-5. Kamrava, S., Tahmasebi, P., Sahimi, M., 2019a. Enhancing images of shale formations by a hybrid stochastic and deep learning algorithm. Neural Netw. 118, 310–320. https://doi.org/10.1016/J.NEUNET.2019.07.009. Kamrava, S., Tahmasebi, P., Sahimi, M., 2019b. Linking morphology of porous media to their macroscopic permeability by deep learning. Transp. Porous Media. https://doi. org/10.1007/s11242-019-01352-5. Karacan, C.Ö., 2003. An effective method for resolving spatial distribution of adsorption kinetics in heterogeneous porous media: application for carbon dioxide sequestration in coal. Chem. Eng. Sci. 58, 4681–4693. https://doi.org/10.1016/J.CES.2003.05. 002. Karacan, C.Ö., Okandan, E., 2000. Assessment of energetic heterogeneity of coals for gas adsorption and its effect on mixture predictions for coalbed methane studies. Fuel 79, 1963–1974. https://doi.org/10.1016/S0016-2361(00)00055-7. Karsanina, M.V., Gerke, K.M., Skvortsova, E.B., Ivanov, A.L., Mallants, D., 2018. Enhancing image resolution of soils by stochastic multiscale image fusion. Geoderma 314, 138–145. https://doi.org/10.1016/j.geoderma.2017.10.055. Kelly, S., El-Sobky, H., Torres-Verdín, C., Balhoff, M.T., 2016. Assessing the utility of FIB-

4. Conclusion Constructing stochastic shale models can be used to not only predict the intrinsic properties (e.g., permeability and formation factor) of shales, but also study the effects of various components and pores on rock properties. Building such models is also necessary for the situations where the comprehensive data are not available as it poses considerable uncertainty. In this study, a dynamic modeling algorithm of shales considering pore-system evolution based on quartet structure generation set and morphological operation was developed. We firstly analyzed the patterns of the pore-system evolution with the change of the thermal and pressure. Based on these patterns, dynamic shale models of multiple cases with the pore-system evolution were constructed. In these dynamic models, various minerals, OM, cements, globigerinids, and three types of pores (i.e., interP, intraP, and OM pores) are all taken into account. Besides, the fractions of every component, the average sizes of all and OM pores, the size distributions of all and OM pores, and two-point correlation functions of these models were also characterized. The simulation of the gas flow in OM and iOM pores were also conducted in order to reveal the permeability variations of these models. These characterizations indicate that various shale models can be generated using the proposed algorithm. Furthermore, dynamic variations of cementation, interP, intraP, and OM pores in shale models indicated to exhibit the performance of the presented algorithm. In fact, any components and pores in shale models can be changed using the proposed dynamic modeling algorithm. As a result, it was shown extensively that the proposed modeling technique in this study embraces an excellent performance to construct any dynamic shale models. As one of the drawbacks of this study is that not all the components in the realistic shales are considered into our models, but the main representative minerals have been taken into account. In terms of future research, this method will be extended to 3D shale models. Further, the effects of every component and pore on intrinsic properties will be also investigated based on these dynamic shale models. Declaration of Competing Interests The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Author Contributions Yuqi Wu: Methodology, Data curation, Writing- Original draft preparation, Visualization, Investigation, Validation; Pejman Tahmasebi: Conceptualization, Methodology, Validation, Writing Review & Editing, Supervision; Chengyan Lin: Validation, Writing Review & Editing, Supervision, Funding acquisition; Chunmei Dong: Investigation. Acknowledgments This work was funded by the Fundamental Research Funds for the Central Universities (18CX06024A) and Technology Major Project, P.R. China (2016ZX05054012, 2017ZX05009001). References Agrawal, A., Prabhu, S.V., 2008. Survey on measurement of tangential momentum accommodation coefficient. J. Vac. Sci. Technol. A Vacuum Surf. Film. 26, 634–645. https://doi.org/10.1116/1.2943641. Akkutlu, I.Y., Efendiev, Y., Savatorova, V., 2015. Multi-scale asymptotic analysis of gas transport in shale matrix. Transp. Porous Media 107, 235–260. https://doi.org/10. 1007/s11242-014-0435-z.

15

International Journal of Coal Geology 218 (2020) 103368

Y. Wu, et al.

Contribution to understanding gas storage and migration pathways in fine-grained rocks. Am. Assoc. Pet. Geol. Bull. 95, 2017–2030. https://doi.org/10.1306/ 03301110145. Song, W., Yao, J., Ma, J., Couples, G.D., Li, Y., Sun, H., 2018a. Pore-scale numerical investigation into the impacts of the spatial and pore-size distributions of organic matter on shale gas flow and their implications on multiscale characterisation. Fuel 216, 707–721. https://doi.org/10.1016/j.fuel.2017.11.114. Song, W., Yao, J., Ma, J., Sun, H., Li, Y., Yang, Y., Zhang, L., 2018b. Numerical simulation of multiphase flow in nanoporous organic matter with application to coal and gas shale systems. Water Resour. Res. 54, 1077–1092. https://doi.org/10.1002/ 2017WR021500. Sun, M., Yu, B., Hu, Q., Yang, R., Zhang, Y., Li, B., Melnichenko, Y.B., Cheng, G., 2018. Pore structure characterization of organic-rich Niutitang shale from China: small angle neutron scattering (SANS) study. Int. J. Coal Geol. 186, 115–125. https://doi. org/10.1016/j.coal.2017.12.006. Tahmasebi, P., 2018a. Nanoscale and multiresolution models for shale samples. Fuel 217, 218–225. https://doi.org/10.1016/j.fuel.2017.12.107. Tahmasebi, P., 2018b. Multiple point statistics: A review. In: Handbook of Mathematical Geosciences. Springer International Publishing, Cham, pp. 613–643. https://doi.org/ 10.1007/978-3-319-78999-6_30. Tahmasebi, P., 2018c. Packing of discrete and irregular particles. Comput. Geotech. 100, 52–61. https://doi.org/10.1016/J.COMPGEO.2018.03.011. Tahmasebi, P., 2018d. Accurate modeling and evaluation of microstructures in complex materials. Phys. Rev. E 97, 023307. https://doi.org/10.1103/PhysRevE.97.023307. Tahmasebi, P., Kamrava, S., 2018. Rapid multiscale modeling of flow in porous media. Phys. Rev. E 98, 052901. https://doi.org/10.1103/PhysRevE.98.052901. Tahmasebi, P., Kamrava, S., 2019. A pore-scale mathematical modeling of fluid-particle interactions: Thermo-hydro-mechanical coupling. Int. J. Greenh. Gas Control. 83, 245–255. https://doi.org/10.1016/j.ijggc.2018.12.014. Tahmasebi, P., Sahimi, M., 2012. Reconstruction of three-dimensional porous media using a single thin section. Phys. Rev. E - Stat. Nonlinear, Soft Matter Phys. 85, 1–13. https://doi.org/10.1103/PhysRevE.85.066709. Tahmasebi, P., Sahimi, M., 2013. Cross-correlation function for accurate reconstruction of heterogeneous media. Phys. Rev. Lett. 110, 078002. https://doi.org/10.1103/ PhysRevLett.110.078002. Tahmasebi, P., Sahimi, M., 2015. Reconstruction of nonstationary disordered materials and media: Watershed transform and cross-correlation function. 91. pp. 032401. https://doi.org/10.1103/PhysRevE.91.032401. Tahmasebi, P., Javadpour, F., Sahimi, M., 2015a. Multiscale and multiresolution modeling of shales and their flow and morphological properties. Sci. Rep. 5, 16373. https://doi.org/10.1038/srep16373. Tahmasebi, P., Javadpour, F., Sahimi, M., 2015b. Three-Dimensional Stochastic Characterization of Shale SEM Images. Transp. Porous Media 110, 521–531. https:// doi.org/10.1007/s11242-015-0570-1. Tahmasebi, P., Javadpour, F., Sahimi, M., 2016. Stochastic shale permeability matching: Three-dimensional characterization and modeling. Int. J. Coal Geol. 165, 231–242. https://doi.org/10.1016/j.coal.2016.08.024. Wang, M., Wang, J., Pan, N., Chen, S., 2007. Mesoscopic predictions of the effective thermal conductivity for microscale random porous media. Phys. Rev. E - Stat. Nonlinear, Soft Matter Phys. 75, 1–10. https://doi.org/10.1103/PhysRevE.75. 036702. Wu, K., Chen, Z., Li, X., Guo, C., Wei, M., 2016. A model for multiple transport mechanisms through nanopores of shale gas reservoirs with real gas effect–adsorptionmechanic coupling. Int. J. Heat Mass Transf. 93, 408–426. https://doi.org/10.1016/ J.IJHEATMASSTRANSFER.2015.10.003. Wu, Y., Lin, C., Ren, L., Yan, W., An, S., Chen, B., Wang, Y., Zhang, X., You, C., Zhang, Y., 2018a. Reconstruction of 3D porous media using multiple-point statistics based on a 3D training image. J. Nat. Gas Sci. Eng. 51, 129–140. https://doi.org/10.1016/j. jngse.2017.12.032. Wu, Y., Lin, C., Ren, L., Yan, W., Wang, Y., Chen, S., You, C., Zhang, L., 2018b. Digital core modeling based on multiple-point statistics. J. China Univ. Pet. Nat. Sci. 42, 12–21. Wu, Y., Tahmasebi, P., Lin, C., Munawar, M.J., Cnudde, V., 2019a. Effects of micropores on geometric, topological and transport properties of pore systems for low-permeability porous media. J. Hydrol. 575, 327–342. https://doi.org/10.1016/J.JHYDROL. 2019.05.014. Wu, Y., Tahmasebi, P., Yu, H., Lin, C., Wu, H., Dong, C., 2019b. Pore-scale 3D Dynamic Modeling and Characterization of Shale Samples: Considering the Effects of Thermal Maturation. Journal of Geophysical Research : Solid Earth. https://doi.org/10.1029/ 2019JB018309. Wu, Y., Tahmasebi, P., Lin, C., Zahid, M.A., Dong, C., Golab, A.N., Ren, L., 2019c. A comprehensive study on geometric, topological and fractal characterizations of pore systems in low-permeability reservoirs based on SEM, MICP, NMR, and X-ray CT experiments. Mar. Pet. Geol. 103, 12–28. Yang, Y., Yao, J., Wang, C., Gao, Y., Zhang, Q., An, S., Song, W., 2015. New pore space characterization method of shale matrix formation by considering organic and inorganic pores. J. Nat. Gas Sci. Eng. 27, 496–503. https://doi.org/10.1016/J.JNGSE. 2015.08.017. Yang, R., He, S., Hu, Q., Hu, D., Zhang, S., Yi, J., 2016. Pore characterization and methane sorption capacity of over-mature organic-rich Wufeng and Longmaxi shales in the Southeast Sichuan Basin, China. Mar. Pet. Geol. 77, 247–261. https://doi.org/10. 1016/j.marpetgeo.2016.06.001. Yao, S., Wang, X., Yuan, Q., Zeng, F., 2018. Estimation of shale intrinsic permeability with process-based pore network modeling approach. Transp. Porous Media 125, 127–148. https://doi.org/10.1007/s11242-018-1091-5. Yu, H., Chen, J., Zhu, Y., Wang, F., Wu, H., 2017. Multiscale transport mechanism of shale

SEM images for shale digital rock physics. Adv. Water Resour. 95, 302–316. Ko, L.T., Loucks, R.G., Zhang, T., Ruppel, S.C., Shao, D., 2016. Pore and pore network evolution of Upper Cretaceous Boquillas (Eagle Ford–equivalent) mudrocks: Results from gold tube pyrolysis experiments. Am. Assoc. Pet. Geol. Bull. 100, 1693–1722. https://doi.org/10.1306/04151615092. Ko, L.T., Loucks, R.G., Ruppel, S.C., Zhang, T., Peng, S., 2017. Origin and characterization of Eagle Ford pore networks in the south Texas Upper Cretaceous shelf. Am. Assoc. Pet. Geol. Bull. 101, 387–418. https://doi.org/10.1306/08051616035. Ko, L.T., Ruppel, S.C., Loucks, R.G., Hackley, P.C., Zhang, T., Shao, D., 2018. Pore-types and pore-network evolution in Upper Devonian-Lower Mississippian Woodford and Mississippian Barnett mudstones: Insights from laboratory thermal maturation and organic petrology. Int. J. Coal Geol. 190, 3–28. https://doi.org/10.1016/j.coal.2017. 10.001. Krishnan, S., Journel, A.G., 2003. Spatial connectivity: from variograms to multiple-point measures. Math. Geol. 35, 915–925. https://doi.org/10.1023/B:MATG.0000011585. 73414.35. Lewan, M.D., 1991. Primary Oil Migration and Expulsion as Determined by Hydrous Pyrolysis. In: 13th World Pet. Congr. 20-25 October, Buenos Aires, Argentina, pp. 215–223. Li, J.J., Yin, J., Zhang, Y., Lu, S., Wang, W., Li, J.J., Chen, F., Meng, Y., 2015. A comparison of experimental methods for describing shale pore features - A case study in the Bohai Bay Basin of eastern China. Int. J. Coal Geol. 152, 39–49. https://doi.org/ 10.1016/j.coal.2015.10.009. Li, C., Lin, M., Ji, L., Jiang, W., Cao, G., 2018. Rapid evaluation of the permeability of organic-rich shale using the 3D intermingled-fractal model. SPE J. 23, 2175–2187. https://doi.org/10.2118/191358-pa. Liu, B., Schieber, J., Mastalerz, M., 2017. Combined SEM and reflected light petrography of organic matter in the New Albany Shale (Devonian-Mississippian) in the Illinois Basin: A perspective on organic pore development with thermal maturation. Int. J. Coal Geol. 184, 57–72. https://doi.org/10.1016/J.COAL.2017.11.002. Loucks, R.G., Reed, R.M., 2014. Scanning-electron-microscope petrographic evidence for distinguishing organic-matter pores associated with depositional organic matter versus migrated organic matter in mudrock. Gulf Coast Assoc. Geol. Soc. 3, 51–60. Loucks, R.G., Reed, R.M., Ruppel, S.C., Jarvie, D.M., 2009. Morphology, genesis, and distribution of nanometer-scale pores in siliceous mudstones of the Mississippian Barnett Shale. J. Sediment. Res. 79, 848–861. https://doi.org/10.2110/jsr.2009.092. Loucks, R.G., Reed, R.M., Ruppel, S.C., Hammes, U., 2012. Spectrum of pore types and networks in mudrocks and a descriptive classification for matrix-related mudrock pores. Am. Assoc. Pet. Geol. Bull. 96, 1071–1098. https://doi.org/10.1306/ 08171111061. Ma, J., Sanchez, J.P., Wu, K., Couples, G.D., Jiang, Z., 2014. A pore network model for simulating non-ideal gas flow in micro- and nano-porous materials. Fuel 116, 498–508. https://doi.org/10.1016/j.fuel.2013.08.041. Ma, P., Lin, C., Zhang, S., Dong, C., Zhao, Y., Dong, D., Shehzad, K., Awais, M., Guo, D., Mu, X., 2018. Diagenetic history and reservoir quality of tight sandstones: A case study from Shiqianfeng sandstones in upper Permian of Dongpu Depression, Bohai Bay Basin, eastern China. Mar. Pet. Geol. 89, 280–299. https://doi.org/10.1016/J. MARPETGEO.2017.09.029. Mason, G., Morrow, N.R., 1991. Capillary behavior of a perfectly wetting liquid in irregular triangular tubes. J. Colloid Interface Sci. 141, 262–274. https://doi.org/10. 1016/0021-9797(91)90321-X. Mehmani, A., Prodanović, M., 2014. The application of sorption hysteresis in nano-petrophysics using multiscale multiphysics network models. Int. J. Coal Geol. 128–129, 96–108. https://doi.org/10.1016/j.coal.2014.03.008. Milliken, K.L., Rudnicki, M., Awwiller, D.N., Zhang, T., 2013. Organic matter-hosted pore system, Marcellus Formation (Devonian), Pennsylvania. Am. Assoc. Pet. Geol. Bull. 97, 177–200. https://doi.org/10.1306/07231212048. Milliken, K.L., Ko, L.T., Pommer, M., Marsaglia, K.M., 2014. Sem petrography of eastern Mediterranean sapropels: Analogue data for assessing organic matter in oil and gas shales. J. Sediment. Res. 84, 961–974. https://doi.org/10.2110/jsr.2014.75. Munawar, M.J., Lin, C., Cnudde, V., Bultreys, T., Dong, C., Zhang, X., De Boever, W., Zahid, M.A., Wu, Y., 2018. Petrographic characterization to build an accurate rock model using micro-CT: Case study on low-permeable to tight turbidite sandstone from Eocene Shahejie Formation. Micron 109, 22–33. https://doi.org/10.1016/j.micron. 2018.02.010. Naraghi, M., Javadpour, F., 2015. A stochastic permeability model for the shale-gas systems. Int. J. Coal Geol. 140, 111–124. https://doi.org/10.1016/j.coal.2015.02. 004. Naraghi, M.E., Javadpour, F., Ko, L.T., 2018. An object-based shale permeability model: Non-darcy gas flow, sorption, and surface diffusion effects. Transp. Porous Media 125, 23–39. https://doi.org/10.1007/s11242-017-0992-z. Okolo, G.N., Everson, R.C., Neomagus, H.W.J.P., Roberts, M.J., Sakurovs, R., 2015. Comparing the porosity and surface areas of coal as measured by gas adsorption, mercury intrusion and SAXS techniques. Fuel 141, 293–304. https://doi.org/10. 1016/j.fuel.2014.10.046. Øren, P.-E., Bakke, S., 2002. Process based reconstruction of sandstones and prediction of transport properties. Transp. Porous Media 46, 311–343. https://doi.org/10.1023/ A:1015031122338. Pommer, M., Milliken, K., 2015. Pore types and pore-size distributions across thermal maturity, Eagle Ford Formation, southern Texas. Am. Assoc. Pet. Geol. Bull. 99, 1713–1744. https://doi.org/10.1306/03051514151. Singh, H., Javadpour, F., 2014. Nonempirical apparent permeability of shale. SPE Reserv. Eval. Eng 17https://doi.org/10.1190/urtec2013-127. SPE-170243-PA. Singh, H., Javadpour, F., 2016. Langmuir slip-Langmuir sorption permeability model of shale. Fuel 164, 28–37. https://doi.org/10.1016/J.FUEL.2015.09.073. Slatt, R.M., O’Brien, N.R., 2011. Pore types in the Barnett and Woodford gas shales:

16

International Journal of Coal Geology 218 (2020) 103368

Y. Wu, et al.

2018.07.018. Zhang, X., Tahmasebi, P., 2019. Effects of grain size on deformation in porous media. Transp. Porous Media 1–21. https://doi.org/10.1007/s11242-019-01291-1. Zhao, T., Zhao, H., Ning, Z., Li, X., Wang, Q., 2018. Permeability prediction of numerical reconstructed multiscale tight porous media using the representative elementary volume scale lattice Boltzmann method. Int. J. Heat Mass Transf. 118, 368–377. https://doi.org/10.1016/J.IJHEATMASSTRANSFER.2017.11.004. Zhu, L., Zhang, C., Zhang, C., Zhou, X., Zhang, Z., Nie, X., Liu, W., Zhu, B., 2019. Review article challenges and prospects of digital core-reconstruction research. Geofluids 2019, 1–29. https://doi.org/10.1155/2019/7814180.

gas in micro/nano-pores. Int. J. Heat Mass Transf. 111, 1172–1180. https://doi.org/ 10.1016/J.IJHEATMASSTRANSFER.2017.04.050. Yu, H., Zhu, Y., Jin, X., Liu, H., Wu, H., 2019. Multiscale simulations of shale gas transport in micro/nano-porous shale matrix considering pore structure influence. J. Nat. Gas Sci. Eng. 64, 28–40. https://doi.org/10.1016/j.jngse.2019.01.016. Yuan, G., Cao, Y., Schulz, H.-M., Hao, F., Gluyas, J., Liu, K., Yang, T., Wang, Y., Xi, K., Li, F., 2019. A review of feldspar alteration and its geological significance in sedimentary basins: From shallow aquifers to deep hydrocarbon reservoirs. Earth-Sci. Rev. 191, 114–140. https://doi.org/10.1016/J.EARSCIREV.2019.02.004. Zhang, X., Tahmasebi, P., 2018. Micromechanical evaluation of rock and fluid interactions. Int. J. Greenh. Gas Control 76, 266–277. https://doi.org/10.1016/J.IJGGC.

17