A novel hybrid thermodynamic model for pore size distribution characterisation for shale

Ze-Zhang Song ,Abie Abula ,Jun-Yi Zhao ,Guang-Di Liu *,Ming-Rui Li ,Dai-Lin Yang ,Yun-Long Wang

a State Key Laboratory of Petroleum Resources and Prospecting,China University of Petroleum,Beijing,102249,China

b College of Geosciences,China University of Petroleum,Beijing,102249,China

c The Exploration Department of the Changqing Oilfield Company,PetroChina,Xi'an,Shaanxi,710021,China

d Research Institute of Petroleum Exploration and Development,Southwest Oil and Gas Field Company,PetroChina,Chengdu,Sichuan 610041,China

Keywords:

ABSTRACT

1.Introduction

Despite the steady and rapid growth of non-fossil fuels in the last decade,oil and gas remained the world's dominant energy,accounting for 33.1%and 24.2%of primary energy,respectively(BP statistical review(Ersoy et al.,2019),).China was one of the biggest individual drivers of primary energy consumption.Given the increasing demand for oil and gas,unconventional resources,especially oil and gas from fractured shale formation,attracted much attention(Chen et al.,2017,Hu et al.,2021a,Hu et al.2021b,Hu et al.2018,Li et al.,2018,Wang S.et al.,2021,Wang et al.,2021).In the present research,we chose one of the typical shale formations-the Chang7 shale formation of the Upper Triassic in the Ordos Basin as our research target.

As an ultra-low permeable porous media(whose permeability usually＜0.1 mD),the shale is characterised by complex and heterogeneous pore structure,which is the primary focus of shale reservoir research.In recent years,numerous investigations were carried out to characterise the pore structure quantitatively.Specific surface area(SSA),specific pore volume(PV)and pore size distribution(PSD)are the utmost frequently adopted parameters.However,compared to SSA and PV,the studies on PSD are still not comprehensive,in-depth.Most scholars concentrated on the PSD pattern and prominent PSD peaks(Chen et al.,2016;Labani et al.,2013).

The importance and the urgency of in-depth research of PSD in shale reservoir research reside in the following aspects.(1)PSD determines the ratio of gases of various occurrences.For shale within oil window,shale gas could exist in pore networks in three phases:the free gas in micro-fractures and intergranular pores,the adsorbed gas mainly on the surface of the kerogen and clay minerals and the dissolved gas in the kerogen and bitumens(Curtis,2002).Numerous experiments have claimed that adsorbed gas in organic-rich shales is mainly associated with micropores(Hao et al.,2013),while free gas is more dependent on mesopores and macropores.PSD determines the ratio of micropores,mesopores and macropores;thus,it further decides the ratio of gases of various occurrences.(2)PSD dominates the adsorption mechanism of adsorbed gas.Adsorption of shale gas in the micropore range is complemented under the“pore-filling”mechanism,while that in mesopore range is dominated by the combination of multilayer adsorption on mesopore walls with physical condensation of adsorbate in pores.(3)PSD determines the adsorption capacity of the shale reservoir.It is well documented that micropore-ranged pore networks have a prominent contribution to adsorption capacity.So,a micropore-biased PSD is more beneficial for a greater adsorption capacity.(4)PSD is crucial for the fluid flow in nanoscaled pore networks.Heller et al.(2014)pointed out that pores ranging from tens of nanometers to 100-200 nm are more effective fluid flow paths.Similarly,Javadpour(2009)reported that when the pore size of shale is reduced to less than 100 nm,the apparent permeability significantly deviates from the Darcy permeability as a result of flowing model transformation from conventional Darcy flow to diffusive transport.From this perspective,PSD is an important parameter reflecting the seepage capacity of reservoir fluid in a nano-scaled pore network.Therefore,the research on PSD is also critical for the gas producibility and economic feasibility assessment(Xiong et al.,2015).(5)PSD even affects reservoir properties,such as elasticity and mechanical(Kuila and Prasad,2013).

Given the significance of PSD for the shale reservoir characterisation,multi-methods have been developed for its qualitative and quantitative evaluations,which could be categorised into two groups-direct and indirect methods.Direct methods,including Field Emission Scanning Electron Microscopy(FE-SEM),Focused Ion Beam-Scanning Electron Microscopy(FIB-SEM),CT-scanning,are capable of giving a direct measure of PSD down to nanometer scale.However,direct methods are usually limited by equipment resolution(Yang et al.,2016).Additionally,the accuracy of PSD is greatly affected by pre-sampling processes,which may significantly alter the pore network.Indirect methods consist of High-pressure Mercury Intrusion Porosimetry(MIP),Constant-rate Mercury Intrusion Porosimetry,Nuclear Magnetic Resonance(NMR),Raybased techniques(Small-angle X-ray scattering and Ultra-smallangle neutron scattering),Low-pressure Gas Adsorption(nitrogen,carbon dioxide),et al.Based on specific physical models,data obtained by these indirect methods could be interpreted to PSD.However,the deficiencies of these indirect methods are also noticeable.Under high pressure,MIP may alter the original pore structure.During the pre-sampling process of NMR,saturating the sample with liquid may lead to clay hydration and damage the shale's weak pore system.Ray-based techniques are not yet widely used in China(Lu et al.,2018).Besides,most of these indirect methods can only analyse the interconnected pore networks.

Although various methods have their advantages and shortcomings,SEM and gas adsorption techniques have shown their operating convenience,wide measuring range of pore size and are usually the combination of the first choice to investigate the PSD of shale.IUPAC categorised three types of nano-scaled pores(Sing 1985;Thommes et al.,2015)-micropores(＜2 nm),mesopores(2-50 nm)and macropores(>50 nm),based on pore size obtained by gas adsorption techniques,and this classification is widely adopted in the research of shale reservoir characterisation.Through FE-SEM studies on North American shale gas samples,Loucks et al.(2012)proposed a pore type classification scheme for shale--intergranular pores,intragranular pores and organic pores.Chen et al.(2016)clarified the dominance of pore types on PSD of marine shales in South China--intraparticle pores are larger while OM pores and interparticle pores associated with OM and clay are smaller.

Though scholars extensively utilised the gas adsorption technique and FE-SEM for shale pore network characterisation,some works could still be further refined.Scholars often see the gas adsorption technique as a straight-to-interpret technique.Hence,parameters given by the gas adsorption technique are usually directly adopted to explain pore-structure-related issues.In this process,pore networks are generally simplified as monogeometric,most frequently cylinder-shaped or slit-shaped.However,the same experimental data interpreted by different geometric models may give very different results regarding the PSD characterisation.Apparently,a monogeometric model does not provide a realistic explanation of PSD.Similarly,scholars often use FE-SEM to describe the location,shape,types of pores,and even quantitative characterise PSD.However,image-processing techniques usually measure PSD using equivalent pore width,and the pores'geometry is not seriously taken into account.

Different PSD models are developed based on different pore geometry models.For the PSD calculation of micropores,the Horvath-Kawazoe(HK)model and the Saito-Foley(SF)model are preferred.However,the former model is based on slit-shaped pore geometry,while the latter one,on cylinder-shaped.For PSD calculation over mesopore and part of macropore range,the classical Barett-Joyner-Halend(BJH)and the Dollimore-Heal(DH)models are the most commonly used models in the literature.Both BJH and DH models are based on the Kelvin equation and modified for multilayer adsorption,assuming that pores are cylinder-shaped with open ends.The PSD obtained from these two models are very close to each other.The accuracy of the macroscopic thermodynamic models(BJH,DH,HK,SF)is limited because they assume that the pores are monogeometric and the fluid in the pore is free,with similar thermophysical properties.However,the thermodynamic properties of confined fluids in complex pore networks are quite different from those of free fluids.

DFT(Density Functional Theory)model provides a microscopic interpretation of gas adsorption in micro and mesopores on a molecular level.Complex mathematical modelling of gas-solid and gas-gas(gas-liquid)interactions plus geometrical considerations(pore geometry)leads to the capability to reflect the confined fluid's thermodynamic properties in the pore more realistically.A new DFT model for a specific porous material is established based on a specific kernel-a set of theoretical adsorption isotherms in a wide range of pore sizes,considering pore geometry.Earlier,there were only monogeometric kernels.With the development of technology,the Hybrid Kernel,for example,the slit-cylindrical adsorption kernel,has emerged(Gor et al.,2012;Neimark et al.,2009).From the perspective of geometric heterogeneity,the DFT model has evolved into a model that can simultaneously consider the influence of multiple pore geometry on gas adsorption.

Though the DFT model has various advantages,it is ordinarily inaccessible for individual researchers.All DFT models are based on Kernels or core programs.The Kernels'establishment demands a rich supply of theoretical isotherms obtained from extensive experiments on uniform and regular porous materials.Gasadsorption-equipment manufacturers have commercialised the DFT model to interpret experimental data and calculate PSD from adsorption isotherms.Given this,individual researchers need a method,which could be utilised to gain a more realistic interpretation of nitrogen adsorption using macroscopic thermodynamic models urgently.

Given these challenges,this work aims to move a step further to establish a new model to obtain a more realistic PSD for the shale pore network.We simplified pore space into two geometric types-cylindrical-and slit-shaped.Using Watershed Segmentation by flooding,we separated the two types of pores and gave a quantitative assessment of their contribution to PSD.On this basis,we integrated the typical monogeometric(cylinder-shaped and slitshaped)thermodynamic models to form a novel hybrid model for PSD interpretation using nitrogen adsorption data.By integrating geochemical analysis,mineral composition analysis,SEM observation,quantitative comparison with the results obtained by the DFT model,and fractal analysis,we clarified the validity of the hybrid model.The hybrid model proposed in this paper can better reflect the contribution of pores with different geometry to PSD.The PSD obtained by the hybrid model is closer to the actual situation of pore space.This model has essential application prospects for a better interpretation of shale pore space.

2.Geologic background

The Ordos Basin is located in the western part of the North China Block(Fig.1(a))and has six fundamental tectonic units(Fig.1(b)).Our research area-Longdong District,located in the southwest of the Shanbei Slope,is a frontier area for Changqing Oil Company in continental shale oil research.The Upper Triassic Yanchang Series Chang7 formation(Fig.2)is characterised by well-developed dark/black shale with a great thickness(50-120 m),stable distribution,high content of organic matter(TOC could exceed 10%),medium maturity(Ro-0.82%-1.2%,early mature to mature stage).To date,Chang7 formation is still in the oil window;a large amount of generated hydrocarbon is retained in the Chang7 formation and leads to excellent hydrocarbon potential.Thus,it became a frontier object for shale oil research.The sedimentary micro-facies of Chang7 are mainly semi-deep or deep-lacustrine.The lithological column and our coring photo suggest that the black shale and carbonaceous shale dominate.The logging responses of Chang7 are typical-high resistivity(RT),high acoustic slowness,high CNL porosity,high gamma-ray,low density(see Fig.2).

3.Experiments and methodology

3.1.Sampling,pyrolysis and XRD analysis

Five core samples originating from the Middle and the Lower Chang7 Formation of the Upper Triassic were obtained from five exploration wells in the research area.To ensure sufficient differences in pore space of different samples for comparison,the five samples collected in this work are from wells with sufficient plane distance(one sample is from the northernmost,one is from the southernmost,and the other three samples are from the center of the study area).Sufficient plane distance ensures noticeable differences in sedimentary microfacies,mineral compositions and the abundance of organic matter.Moreover,the depth of different samples varies to ensure that the maturity of the samples is different so that the development degree of organic-matter-hosted pores(cylinders)is different.We crushed the studied core samples to yield particle sizes between 60 and 80 mesh sizes using the Spherical Grinder in the State Key Laboratory of Petroleum Resources and Prospecting.Then,we performed pyrolysis using a CS-230HC carbon/sulfur determinator.After that,we carried out the Soxhlet Extraction on these powder samples.In this process,taking safety into consideration,dichloromethane,which is less toxic than other kinds of chloro-hydrocarbons,was selected as the organic solvent.Each sample had been undergone no less than 4-days'(96-h')extraction to extract the residual hydrocarbon from the samples as entirely as possible.After that,part of the sample was used for the X-ray diffraction analysis(XRD)experiment,which was carried out in Beijing Research Institute of Uranium Geology following the Chinese Oil and Gas Industry Standard(SY/T)5163-2010(Analysis method for clay minerals and common non-clay minerals in sedimentary rocks by the X-ray diffraction).Mineral compositions,such as quartz,clay mineral,feldspar,dolomite,orthoclase,pyrite,are estimated with Panalytical X'Pert PRO.Also,part of the powder samples was oven-dried in Beijing Centre for Physical ＆ Chemical Analysis at 65?C(149?F)for approximately 72 h until a constant weight was achieved.The weight changes were monitored/recorded every 6 h.This portion of the sample was prepared for the LTNA experiment.

Fig.1.Geological background of the present work,showing(a)the location of the Ordos Basin in China;(b)the location of the research area-the Longdong District in the Ordos Basin(modified from Song et al.(2018))and the distribution of wells,from which shale samples were obtained.

Fig.2.The typical logging responses and lithology of Chang7 formation of Upper Triassic in Longdong District(Taking the well Z22 as an example).

3.2.Low-temperature nitrogen adsorption experiment(LTNA)

In this research,the physisorption of N2on shale samples was carried out with the help of Quantachrome Autosorb-iQ MP-a highly sophisticated automatic gas adsorption analyser,in the Beijing Centre for Physical ＆ Chemical Analysis.This set of adsorption analysers has been designed to be the most flexible,versatile,and modular surface area and pore size distribution analyser on the market.Autosorb-iQ is capable of characterising the pore structure of shale samples precisely to nano-scale,providing specific surface area analysis from 0.0005 m2/g to known upper limit,pore size distribution analysis from 0.35 nm to 500 nm(using nitrogen as adsorbate)and pore volume analysis as precise as 0.0001 cubic centimetres.Prior to the analysis,shale samples were outgassed under a turbomolecular pump vacuum for 8 h under 110?C.The outgassing process was completed when the pressure and the weight of the sample in the sample holder were stabilised.Then nitrogen with a purity higher than 99.999% was used as adsorbate,and the physisorption was carried out at liquid nitrogen temperature(-195.85?C/77.3 K).The equilibration time was set to 90 s,and the sample weight ranged from 2 to 5 g.The Autosorb-iQ automatically recorded adsorbed gas volume with increasing/decreasing relative pressure by applying the gravimetric method.

3.3.A novel hybrid model for shale PSD description

3.3.1.Step 1:analysis using monogeometric thermodynamic models

In this study,we simplified the pore network into a dualgeometric model-a combination of cylinder-shaped and slitshaped pores.Furthermore,we assume that all the cylindershaped pores are open on two ends,and the slit-shaped pores are open on all sides.

Like the BJH model,PSD calculation from nitrogen adsorption/desorption data utilising classical thermodynamic theory is implemented based on an assumption:two mechanisms-physical adsorption on the pore walls and capillary condensation in the inner capillary volume,simultaneously determine the equilibrium during adsorption/desorption(Barrett et al.,1951).Thus,the pore radius calculated from the monogeometric thermodynamic model is the sum of the physically adsorbed layer's thickness and the capillary radius calculated by the classical Kelvin equation.

For cylinder-shaped pores,in the area where capillary condensation is in presence,pore radius(r)is the sum of the thickness of adsorption layer(t)at the arbitrary pressure and Kelvin radius(rk)of the meniscus(Fig.3).The thickness of the adsorption layer(t)can be calculated by Equation(2),while the Kelvin radius can be calculated by Equation(3).

Fig.3.Geometric parameters for PSD calculation of cylinder-shaped and slit-shaped pores.

where x is the relative pressure,the boiling temperature of nitrogen T=77.3 K,the molar volume of liquid adsorptive vm=34.65 mL/mol,the surface tension of nitrogen at its boiling pointγ=8.85×10-3J/m2,φ=0,the universal gas constant R=8.314×10-23J/mol/K.

The adsorption and desorption processes in slit-shaped pores are different from those in cylinder-shaped pores.Firstly,capillary condensation will not occur in the adsorption process until relative pressure reaches 1.Secondly,the meniscus generated in slit-shaped mesopores is hemicylindrical,while in cylinder-shaped pores is hemisphere surface.Finally,in the desorption process,when capillary condensation is in presence,pore width(r)is the sum of the thickness of the adsorption layer(t)at the arbitrary pressure and core radius(rk)of the meniscus(Fig.3).

When the gas-liquid equilibrium occurs in parallel slit-shaped pores with a diameter of d(d=2×(r+tk)),the Kelvin radius could be obtained by Equation(4)while t and r can still be calculated using Equations(2)and(1):

Hence,under the same relative pressure,the Kelvin radius of the cylinder-shaped pore is twice as large as that of the parallel slitshaped pore.That is to say,the capillary effect of the cylindershaped pore is more robust than that of the parallel slit-shaped pore.

Also,the pore geometry affects the calculation of pore structure parameters,like SSA and PV.

As in the parallel slit-shaped pore,capillary condensation does not occur during the adsorption process;we select the desorption branch for the PSD calculation.Let V(r)be the distribution function of pore volume,calculating PSD is to find the pore volume in every step of pressure decrease during desorption,in the range of(ri,ri-1)accordingly,which in integral forms gives:

Let v(r)be the total desorption amount in a specific pressure range composed of two parts:de-condensation and adsorption layers'thinning(the desorbed amount).Then,the desorption amount of gasΔviobtained from the experiment could be expressed like:

So,when relative pressure decreases from 1 to xi,the evaporated amount of condensed liquid equals to:

As the evaporated amount of condensed liquid equals the released pore volume,Equation(7)equals Equation(8).

Thus,if we consider two steps of pressure-decrease,Equation(9)could be obtained.

By solving the above equations,the following results are obtained:

Let ri=(ri-1+ri)andti=(ti-1+ti),then,apply the integral mean value theorem to the first part on the right side of the above formula:

LetΔti=(ti-1-ti),and apply the integral mean value theorem to the second part of Equation(10):

Integrating Equation(8-11),the following equation is obtained:

Thus,the obtained recursive Eq.(13)could be used for PSD calculation for the ideal cylinder-shaped pore network.

When it comes to the ideal slit-shaped pore network,define the distance between parallel slits as d;similarly,the recursive Eq.(14)for PSD calculation could be obtained:

The detailed mathematical derivation of the PSD calculation of monogeometric cylinder-shaped and slit-shaped pore networks could refer to the work of Liu et al.(2017);Jin and Huang(2015).

3.3.2.Step 2:pore geometric segmentation using watershed by flooding

In geology,a watershed is a dividing ridge between adjacent catchment basins(drainage areas).In the study of image processing,Watershed Algorithms treats the grey-scale image as a topographic map with the grey-scale value of each pixel in the image representing its elevation,primarily for segmentation purposes.There are different kinds of Watershed Algorithms,like Watershed by flooding(Beucher and Lantuejoul,1979),Watershed by topographic distance(Meyer,1994),Watershed by the Drop of Water Principle(Cousty et al.,2009),Inter-pixel Watershed(Beucher and Meyer,1992)and Meyer's Flooding Algorithm(Barnes et al.,2014).The idea of Watershed by flooding was first introduced by(Beucher and Lantuejoul,1979),and it is the most frequently utilised Watershed Algorithm.This algorithm floods catchment basins from user-defined markers until basins attributed to different markers meet on watershed lines.In many cases,markers are chosen as the local minima of the image,from which basins are flooded.The formed watershed lines are the boundaries between catchment basins.

Watershed Algorithms have an excellent response to weak edges.Thus,this algorithm can extract pore structure features from SEM images with high precision and clarity.The closed catchment basin obtained by Watershed Algorithm makes it possible to analyse the image's regional characteristics.When it comes to pore structure analysis on SEM images,the formed catchment basins represent pores and throats.By analyzing each catchment basin's geometric characteristics,we obtain the geometric information of every pore and throat,making it possible for the geometric analysis of pore structure.

The high sensitivity of the Watershed Algorithm to edges also brings some problems,one of them-over-segmentation.Thus,before analyzing segments obtained from Watershed Algorithm,some pre-processing(the first four steps shown in Table 1)should be done for the SEM image.The entire process for pore structure analysis on SEM image using Watershed Algorithm is shown in Table 1 and Fig.2.

Table1 Pore structure analysis on SEM image using Watershed by flooding.

3.3.3.Step 3:creation of the novel hybrid model

Based on the monogeometric thermodynamic models and image processing utilizing Watershed Algorithms on typical SEM images,we combined the results of the two to form a new PSD model for shale pore structure investigation-the Hybrid Model,shown as follows:

where the?cyl?and?slitare the ratio of cylinder-shaped and slitshaped pores obtained using Watershed Algorithms on typical SEM images.

Before further interpretation,we often encounter three problems.

Which branch to utilise,the adsorption or the desorption branch?Theoretically,the Kelvin equation is based on a cylindrical gas-liquid interface;thus,the adsorption branch better reflects the real cylinder-shaped pore network.It is the same for the ink-bottle pores,whose meniscus in the thin neck determines the balance between capillary condensation and evaporation.However,the desorption branch is better for slit-shaped pores.Because only indesorption there will be a meniscus corresponding to the Kelvin equation.In reality,especially when it comes to the complex pore network of shale formation,the pore geometry consists of not only cylinder-shaped but also slit-shaped pores.Therefore,we selected the desorption branch for the PSD calculation in this paper.

Fig.4.Graphic flow for pore structure geometric analysis using Watershed Algorithm by flooding(Sample ID:Z22-1584.65 m).

In which form to better present PSD?For shale formation,large amounts of gas exist in the pore network in the adsorbed state,which is determined by the specific surface area(SSA).Hence,in the presentation of PSD,we chose the plot of differential pore volume(dV/dD)versus pore size(d)for its unique advantages in showing the contribution of pores of different scales to specific surface area(Meyer and Klobes,1999).

The range of desorption data for PSD interpretation?As many researchers recommended,the thermodynamic models should be utilized primarily for mesopores but not micropores.Furthermore,we often encounter a fake peak at 4 nm in the PSD obtained from the desorption branch(Li et al.,2015).As Groen et al.(2003)explained,the tensile strength effect arouses the fake peak.So we took 5 nm as the lower limit for investigation.Also,Heller et al.(2014),Javadpour(2009)recommended that when the pore size is greater than 100 nm,the shale pore network's gas flow starts to behave like Darcy flow.Furthermore,the research on the shale pore network primarily focuses on the non-Darcy region.Considering these,we chose the range from 5 nm to 100 nm for the PSD interpretation.

3.4.Fractal characterisation using LTNA data

Quantifying the heterogeneity of pore structure of shale using the low-pressure adsorption/desorption data is an essential and well-established process(Liang et al.,2015;LIU et al.,2018;Liu et al.,2015;Yang et al.,2014).In this work,we introduce fractal characterization to verify if the hybrid model could better understand pore space heterogeneity.

Although various models are available to calculate the fractal dimension of pore structures of shale,the FHH model(Avnir and Jaroniec,1989;Pfeifer et al.,1989)is most commonly and widely used.The FHH model could be represented using Equation(16).

V is the volume of adsorbed gas at equilibrium pressure,cc/g;P0and P represent the adsorbate's saturated vapour pressure and the equilibrium pressure,respectively,MPa.K is the power-law exponent related to heterogeneity and the adsorption mechanism.The relationship between K and the fractal dimension(D)could be described as follows(Pfeifer et al.,1989).

Fig.5.The application of the FHH model in adsorption and desorption branches(Sample ID:Z22-1584.65 m).

As shown in Fig.5,when ln(ln(P0/P))=0.91,the fractal of the adsorption section usually breaks into two parts:the corresponding relative pressure and pore size in the dividing point are about 0.05 and 1.6 nm,respectively.It shows distinct fractal features in these two different pore size ranges.D1(P0/P＜0.05,d＜1.6 nm,whose adsorption behaviour is dominated by van der Waals force)and D2(P0/P>0.05,d>1.6 nm,whose behaviour is dominated by capillary condensation)represent the surface fractal dimension and pore structure fractal dimension,respectively(Hazra et al.,2018).The fractal dimension calculated by D=3K+3 is less than 2,deviating from the normal gas adsorption fractal dimension(2＜D＜3).Therefore,D=3+K is selected to calculate the fractal dimension.Furthermore,the obtained D1of the sample is less than 2,which demonstrates the inadequacy of the FHH model for pore structure characterisation in the corresponding range.Hence,we only use the mesopores-macropores'section(D2)of fractal distribution for pore structure characterisation.When the pore size is larger than 1.6 nm,the generalised fractal dimensions from adsorption and desorption branches are very close(2.591 and 2.586,respectively,with a difference of 0.005)to each other.The goodness of the adsorption branch's linear fitting is slightly higher than that of the desorption branch.Thus,this study selected the adsorption branches(except the microporous range)for FHH fractal analysis.

4.Results

This paper selected five typical samples of the Chang7 formation from five exploration wells for shale pore structure investigation.Geochemical analysis,mineral composition analysis,image processing using Watershed by flooding for SEM images,PSD calculation using novel Hybrid model and pore structure heterogeneity analysis were carried out.

Fig.6.Typical isotherms obtained in the research area with linear(left),semi-logarithmic coordinates(right),sample A-Z22 1584.65 m.According to IUPAC classification,the shapes of isotherms obtained in this paper are all Type II or Type IV isotherms.The shape of the adsorption equilibrium isotherm is related to the pore structure of the material.

Fig.7.Five typical adsorption-desorption isotherms in the research area.

4.1.Geochemical parameters and mineral composition

As shown in Table 2,all the selected samples are typical organicrich shale samples.The TOC of the samples ranges from 2.19%to as high as 16.36%,with an average value of 6.69%.The Tmaxobtained from pyrolysis varies from 445?C to 456?C,indicating that all these samples are located in the oil window-mature stage.

Table2 Geochemical parameters for the five selected samples.

Table 3 shows that all samples are dominated by quartz and clay minerals while free of carbonate minerals.Two(sample A and C)out of five samples are rich in clay minerals,reaching 71%and 51.8%,respectively.In contrast,sample E is dominated by quartz with the lowest clay content(20.2%).Sample B and sample D are abundant in pyrite,up to 17.30% and 14.70%,respectively.

Table3 Mineral composition of the five selected samples.

Fig.8.Comparison of different models for PSD calculation(the desorption data were utilized,Sample ID:Z22-1584.65 m).

4.2.Characteristics of isotherms and hysteresis loop

The isotherm shape reveals the surface properties,the pore size distribution of adsorbent and the interactive properties between the adsorbent and adsorbate.Though there are some differences between different samples regarding the shape of isotherms in the research area,the isotherms are almost all anti-S shaped(Fig.6).In general,all these isotherms reflect the following three general characteristics.

All the isotherms resemble the Type IV(a)isotherm with hysteresis loop according to IUPAC recommendations(2015),indicating the existence of both mesopores and macropores,with a dominance of mesopores.Low-pressure zone(P/P0＜0.1).The adsorption volume increased rapidly(approximately exponential increase),whether in the linear or semi-logarithmic coordinates,as we can see from the partial enlarged image(upper-right of the isotherms).The first inflexion point occurs when P/P0ranges from 0.02 to 0.05.The steep uptake of adsorbate is due to enhanced intermolecular force from narrow pore walls.In this zone,micropore-filling or monolayer adsorption should be completed.Intermediate pressure zone(0.1＜P/P0＜0.8).The volume of adsorbed gas increases slowly and steadily,resulting in a plateau.In the relative pressure zone 0.4＜P/P0＜0.8,for some samples,the adsorption and the desorption curves do not coincide with each other,forming a hysteresis loop.Capillary condensation in mesopores is the leading cause of the hysteresis loop.High-pressure zone(P/P0>0.8).We can observe a steep uptake of adsorbates,and the hysteresis loop gradually close.The second inflexion point occurs when P/P0ranges from 0.85 to 0.95.

The hysteresis loop aroused by metastable meniscus during capillary condensation and pore blockage is related to the pore geometry,PSD,and pore network heterogeneity.We further elaborated on the similarities and differences of those five adsorptiondesorption isotherms,shown as follows.Sample A is characterised by a sharp step-down of the desorption branch,representing the Type H2 hysteresis loop.The forced closure is called the“tensile strength effect”(Tripathy et al.,2019)and is brought about by mainly ink-bottle-shaped venting pores with a smaller amount of slit-shaped pores(Liu et al.,2015).Sample B,C,D also have a significant hysteresis loop without a sharp step-down of the desorption branch,resembling the Type H3 hysteresis loop,suggesting that the ink-bottle-shaped pores along with slit-shaped pores dominate among mesopores.The hysteresis loop of sample E is relatively insignificant.The adsorption and desorption curves almost coincide with each other in most parts of the isotherm.

It is worth noting that in some samples,the adsorption and desorption branches remain unclosed(almost parallel)at the lowpressure range.This phenomenon resulting from swelling,namely low-pressure hysteresis,was well interpreted by(Gregg and Sing,1982).The closure of the adsorption and desorption branches at high P/P0in Sample-A,B,C and D isotherms indicates that the large pores in these samples are mainly non-venting pores closed at one end(Sing,1985;Thommes et al.,2015).

4.4.PSD interpretation based on the hybrid model

As shown in Fig.4,after pore geometric segmentation using Watershed by flooding,we concluded that the pore network of sample A consists of 54.75%cylinder-shaped pores and 45.25%slitshaped pores.Based on the monogeometric thermodynamic models(slit and cylindrical),we created a specific hydrid model for sample A and interpreted its PSD.

As we can see from Fig.8,in the range from 5 nm to 100 nm,the PSD obtained from classical thermodynamic models,no matter cylindrical(purple,PSDCyl.)or slit model(green,PSDSlit),deviate from the PSD from the DFT model(pink,PSDDFT).In contrast,the PSD obtained from the novel hybrid model(PSDHybrid)seems to coincide with the PSDDFTvery well.It indicates that the hybrid model is more conducive to explaining the contribution of pores with different sizes to the specific surface area.

Moreover,a closer comparison of various PSDs will notice that the PSDs obtained from the classical thermodynamic models(PSDCyl.,PSDSlit)deviate from the PSDDFTprimarily in the mesopore range(5 nm-50 nm).When pore size decreases,the difference between thermodynamic models and the DFT model increases.

Quantitatively,we compared the PSDs from the two thermodynamic models(PSDCyl.,PSDSlit)and the hybrid model(PSDHybrid)with PSDDFT.Here comes a question:how shall we compare them?The ratios of cylinder-shaped and slit-shaped pores are obtained by Watershed-image-processing,which is the ratios of cross-sectional areas of the cylinder-shaped and slit-shaped pores.If we assume the height of cylinders and slits are the same,we could convert the ratio of the cross-sectional area into the ratio of volume.Therefore,to compare the three PSDs with PSDDFT,we should compare the area enclosed by different PSD curves and the X-coordinate axis.

As we see from Fig.4,the cylinder-shaped pores are more dominant than slit-shaped pores in the pore network of the sample Z22-1584.65 m.Hence,from Fig.4 and Table 4,we can see that the PSD obtained from the cylindrical monogeometric thermodynamicmodel(PSDcyl.,deviation from PSDDFT:36.74%)is more accurate than that from the slit monogeometric model(PSDslit,deviation from PSDDFT:89.07%).On this basis,the hybrid model further reduced the deviation by 16.55%,which further enhanced PSD interpretation precision.

Table4 Quantitative comparison of PSDs obtained from different thermodynamic models with that from the DFT model(Sample A:Z22-1584.65 m).

Fig.9.Watershed Segmentation by flooding results on SEM images for another four representative samples,showing the ratio of cylinder-shaped and slit-shaped pores in each sample.?c and?s represent the ratio of cylinder-shaped pores and slit-shaped pores in pore network obtained by image processing using Watershed Algorithm by flooding.

It proves three things.First of all,it is of great significance considering pore geometry when interpreting PSD.Different thermodynamic models may give very different results.Secondary,when monogeometric models are used,the model corresponding to the pore space's dominating pore geometry will give more realistic results.Last,the hybrid model could further enhance the precision of PSD interpretation over monogeometric thermodynamic models.

5.Discussion

5.1.Watershed Segmentation help understand the pore geometry realistically and quantitatively

With Watershed Segmentation by flooding,we further analyzed the pore geometry of the other four representative samples,the results were shown in Fig.9.We learned that the cylinder-shaped and slit-shaped pores are almost equally distributed in the pore network with a slight bias to cylindrical pores in sample A.We also found an apparent dominance of cylinder-shaped pores in samples B(85.6%)and D(72.78%),while slit-shaped pores dominate samples C(73.76%)and E(85.68%).

Does the ratio obtained by Watershed algorithms realistically reflect the pore space of these samples?To verify it,we studied the pore structure of the five representative samples directly with the help of FE-SEM(Fig.10)and summarized these samples'microscopic characteristics,shown as follows.

In sample E,there are widely distributed micro-fractures on a mesoporous scale with few pores.The micro-fractures are mainly developed at the contact interface between organic matter and minerals.Micro-fractures primarily exist slit-shaped.This is why the slit-shaped dominance in sample E(85.68%).The organicmatter-hosted(OM-hosted)pores are poorly developed in samples D and E,corresponding to the low adsorption volume in the low P/P0range(Fig.7).Compared with sample E,micro-fractures are undeveloped in sample D.In contrast,the dissolution pores,which exist cylinder-shaped,are very developed.From this perspective,cylinder-shaped pores dominate in sample D due to the developed dissolution pores(72.78%).

The OM-hosted pores(mainly cylinder-shaped),intercrystalline pores of clay minerals(cylinder-shaped and slit-shaped),dissolution pores(cylinder-shaped)and micro-fractures are all developed in sample A.Furthermore,most of the pores and micro-fractures are on the scale of mesopores.Their co-existence leads to an almost equal ratio of cylinder-shaped and slit-shaped pores.

Fig.10.Pore structure characteristics under SEM for the five types of representative samples.

The OM-hosted pores(cylinder-shaped)dominate in sample B.We observed not only mesopore-scaled OM-hosted pores but also micropore-scaled,which result in steeper adsorption in the low P/P0range(Fig.7).It is worth noting that we also observed a small number of micro-fractures in organic matter in sample B.In contrast,the clay-hosted pores are not so developed in sample B.The developed OM-hosted pores result in a dominance of cylindershaped pores(85.60%).Unlike the situation in sample B,clayhosted pores provide the most pore space for sample C.The mesoporous and macroporous clay-hosted pores exist mainly in slitshaped(Yang et al.,2017)and mesopore scale,which leads to a slit-shaped dominance(73.63%)of the pore network.

Through microscopic observation of pore space,we found the relationship between pore types and pore geometry.OM-hosted pores and dissolution pores are mainly cylinder-shaped,while clay-hosted pores are primarily slit-shaped.Micro-fractures,whether developed in OM or minerals,occur in slit-shaped.The ratio obtained by Watershed Algorithms is consistent with the qualitative observation of pore geometry under FE-SEM.

Though microscopic observation helps understand why a specific pore shape dominates a specific sample's pore space,it has a limitation-it only reflects pore geometric characteristics in specific fields of view,but not that of the whole rock.Therefore,it is necessary to verify the findings above with more evidence further.Thus,we further analyzed the results of XRD and pyrolysis.

To some extent,maturity determines the degree of the development of OM-hosted pores.Pyrolysis shows that the Tmaxof samples D and E is lower than that of samples A,B and C,which indicates that the maturity of samples D and E is lower than that of samples A,B and C.In maturation,the solid kerogen turns into hydrocarbon and generates pores,which are mainly cylindershaped.Thus,a higher maturity of samples A,B and C is consistent with relatively developed OM-hosted pores in these samples.Among all these samples,sample B is characterized by the highest TOC content and the highest maturity(Tmax).Accordingly,we observed micropores,mesopores,even micro-fractures in OM.Also,sample B's clay content is relatively low(30.5%),and the clayhosted pores are not developed.The dominance of cylindershaped pores in sample B is due to the high content and maturity of OM.The clay content of sample C is higher than 50%.Also,we observed widely distributed clay-hosted pores in the pore space of sample C.Thus,though the OM-hosted pores are also developed in sample C,the slit-shaped pores(72.78%)dominate the pore space of sample C.The highest felspar content characterizes sample D-The total content of orthoclase and plagioclase reached 20.4%.It is almost three times the feldspars'content in samples A,B and C,and twice in sample E.The acidic fluids generated by the process of hydrocarbon generation will dissolve feldspar and generate dissolution pores.Although the TOC content of sample D is relatively high,its maturity is relatively low;simultaneously,its clay mineral content is relatively low,and the clay-hosted pores are not developed.The factors mentioned above result in sample D being dominated by noticeable dissolution pores(cylinder-shaped,72.78%).In contrast,sample E is characterized by the highest quartz content(59.3%).With the lowest TOC content(2.19%),lowest maturity,lowest clay content(20.2%),we did not observe any OMhosted or clay-hosted pores in sample E.Micro-fractures are primarily developed in the interface between quartz and other minerals.Thus,slit-shaped pores dominate the pore space of sample E.

Above all,it seems safe to conclude that Watershed Segmentation can provide a quantitative description of the pore geometry,and the result could realistically reflect the geometric composition of pore space.

5.2.The novel hybrid model provides PSD interpretation closer to the DFT model

Section 4.4,taking sample A as an example,demonstrated that PSDHybridis closer to PSDDFTthan PSDcyl.and PSDslit.So,does the hybrid model also valid for other samples?We created the hybrid model for the other four samples based on Watershed Segmentation on SEM images to verify it further.We compared the obtained PSDs from three thermodynamic models(PSDHybrid,PSDcyl.,PSDslit)with that from the DFT model(PSDDFT);the results are shown in Fig.11,and the overall quantitative analysis is shown in Table 5.

Qualitatively,we can see from Fig.11 that PSDHybridof all samples is closer to PSDDFTthan PSDcyl.and PSDslit.Quantitatively,the deviation of the hybrid model from the DFT model is the lowest,smaller than that of the cylindrical and slit model from 5.06% to 68.88%(Table 5).It indicates that the hybrid model can further improve the accuracy of PSD interpretation than monogeometric thermodynamic models.The hybrid model is valid for different samples with various pore spaces.

Fig.11.Comparison of different models for PSD calculation for the other four representative samples.

Table5 Deviation of three models to DFT models and the improvement of the hybrid model over monogeometric models.

5.3.The PSDHybrid help better understand the heterogeneity of the pore network

To clarify if the newly established hybrid model is valid or even better for heterogeneity interpretation,based on the FHH fractal model,we analysed all samples'nitrogen adsorption data,as shown in Fig.12.The results of piecewise linear fitting are satisfying,with fitting coefficients higher than 97%,demonstrating the FHH model's validity.

The D2of sample B(corresponding equivalent pore diameter is from 1 nm to 15 nm)is the highest,while the D2of sample A ranked second among all samples.As we know that OM-hosted pores are the primary space for adsorbed gas,which may account for as high as 80% of total gas in shale.OM's high capacity for shale gas adsorption is due to its highly heterogeneous pore surface,which leads to relatively high heterogeneity.The well-developed OMhosted pores in sample A＆B is the main reason for the high heterogeneity.Interestingly,samples E and D,whose pore network is dominated by micro-fractures and dissolution pores,have relatively high D2.The heterogeneity of sample C is the lowest.From the pore type's perspective,the prominent clay-hosted pores,undeveloped OM-hosted pores may be the reason.From the pore size perspective,the clay-hosted pores,whose size is larger than OMhosted pores,may contribute to a relatively low heterogeneity.

For quantitative comparison's sake,we calculated the PV utilizing monogeometric models(cylindrical and slit)and the hybrid model to compare the correlation between obtained PV with D2.As shown in Table 6,the obtained PV by the hybrid model of all the samples falls in between that by cylindrical and slit models.Moreover,the correlation relationship between D2and PV obtained by the hybrid model(Pearson's R=-0.6696)is higher than that by monogeometric thermodynamic models(Pearson's R=-0.6375 and-0.6380,respectively).Previous researches have demonstrated that the fractal dimension D2of shale negatively correlates with PV.A closer relationship between D2and PV demonstrated that the hybrid model could better present the relationship between the two.

Table6 Comparison of PV obtained by cylindrical,slit and hybrid model and their correlation with fractal dimension D2.

Fig.12.Fractal analysis of five representative samples for pore structure characterisation of continental shale(Two vertical green segments divide each graph into three sections,with macropores on the left,mesopores in the middle,and micropores on the right.K1 and K2 represent the slope of segment 1,and segment 2,D1 and D2 represent their fractal dimensions.There are three segments for sample B,while there is only one segment for sample C.When the fractal features in different pore sizes are similar,the calculated fractal dimension in different pore size ranges tend to be the same(Sample C);when there are two(Sample A,D,E)or three(Sample B)pore size ranges with distinct fractal features,there will be two/three intervals with distinct calculated fractal dimensions.).

5.4.Remained challenges

In this work,we proposed a novel hybrid model based on the combination of Watershed segmentation for geometric analysis ofpore space and monogeometric thermodynamic models for PSD calculation.We utilised the hybrid model for PSD analysis of Chang7 shale and clarified its validity and improvement for PSD interpretation.From PSD's perspective,it is appropriate to say that we made a step further to obtain a more realistic PSD that can more precisely reflect the pore geometry of shale pore space than monogeometric models.However,there are still some remained challenges.

The proposed hybrid model is developed based on monogeometric thermodynamic models.Thus,it has the shortcomings of conventional thermodynamic models--not suitable for analyzing micropores and extremely small mesopores.In this article,we set a limit of 5-100 nm.The obtained PSD is closer to the DFT model than that obtained by monogeometric thermodynamic models only in this pore size range.How to expand the practical application scale of the model is the focus of further research.

Watershed Segmentation by flooding helps us analyze whether the geometry of each pore is biased to cylinder or slit.It is a dichotomous scheme,but the geometric types of rock pore space are far more complex than the two.To further refine the research on pore morphology in future research is a difficult point that we need to tackle.

Also,it is worth noting that a good understanding of the pore network is the key to successfully utilizing the method described in this article.Since the ratios of slit-shaped pores and cylindershaped pores are extracted from SEM pictures,selecting SEM pictures representing the pore space's morphological characteristics is particularly important for obtaining correct results.If necessary,multiple SEM pictures can be preferably analyzed,and the average value shall be taken as the final ratio.This part of the work is currently more dependent on the researcher's understanding of the pore space.How to realize the automation and refinement of this part of work is the focus of the next step.

6.Conclusions

In the present work,by integrating two thermodynamic monogeometric models(cylinder-shaped and slit-shaped)and Watershed Segmentation by flooding,we proposed a hybrid model for interpreting shale PSD in the range of 5-100 nm by lowtemperature nitrogen adsorption data.The general process for generating a hybrid model is:Firstly,analyzing the nitrogen adsorption data utilizing two monogeometric models(cylindrical and slit)to generate PSDcyl.and PSDslit;Secondly,carrying out pore geometric segmentation using Watershed by flooding on typical SEM images to obtain the ratio of slit-shaped(?s)and cylindershaped pores(?c).Lastly,combining the results of the two to form a hybrid model based on Equation(15).Based on pyrolysis,XRD,FE-SEM observation,quantitative comparison with the results obtained by the DFT model,and fractal analysis,we discussed and verified the validity of the hybrid model.

The hybrid model proposed in this work could better reflect the pore space's real geometry and provide a more realistic PSD.Moreover,the hybrid model is valid for different samples with various dominating pore geometry.Compared with thermodynamic monogeometric models,PSDs obtained from the hybrid model are closer to that from the DFT model.Also,the fractal analysis suggested a closer relationship between D2and PV obtained by the hybrid model,demonstrating that the hybrid model could better reflect the heterogeneity of the pore network.

Acknowledgements

This research was financially supported by the National Key R＆D Program of China(Grant No.2017YFC0603106),the Youth Program of National Natural Science Foundation of China(Grant No.41802148);and the State Key Laboratory of Petroleum Resources and Prospecting(Grant No.2462017YJRC025,Grant No.PRP/indep-04-1611).We would also thank the PetroChina Changqing Oilfield Company for providing core samples and necessary data.