Mesoscale structure, mechanics, and transport properties of source rocks' organic pore networks.

Organic matter is responsible for the generation of hydrocarbons during the thermal maturation of source rock formation. This geochemical process engenders a network of organic hosted pores that governs the flow of hydrocarbons from the organic matter to fractures created during the stimulation of production wells. Therefore, it can be reasonably assumed that predictions of potentially recoverable confined hydrocarbons depend on the geometry of this pore network. Here, we analyze mesoscale structures of three organic porous networks at different thermal maturities. We use electron tomography with subnanometric resolution to characterize their morphology and topology. Our 3D reconstructions confirm the formation of nanopores and reveal increasingly tortuous and connected pore networks in the process of thermal maturation. We then turn the binarized reconstructions into lattice models including information from atomistic simulations to derive mechanical and confined fluid transport properties. Specifically, we highlight the influence of adsorbed fluids on the elastic response. The resulting elastic energy concentrations are localized at the vicinity of macropores at low maturity whereas these concentrations present more homogeneous distributions at higher thermal maturities, due to pores' topology. The lattice models finally allow us to capture the effect of sorption on diffusion mechanisms with a sole input of network geometry. Eventually, we corroborate the dominant impact of diffusion occurring within the connected nanopores, which constitute the limiting factor of confined hydrocarbon transport in source rocks.

The emergence of disruptive technology (1, 2) in the field of petroleum engineering has improved access to hydrocarbons from source rock formations. As a result, the historically declining production of natural gas has been on the rise. The stimulation technique of hydraulic fracturing aims to create a fracture network accessing the hydrocarbons and to force their migration to the crack’s surface (2, 3). However, as the bimodal production curves suggest, the large-scale deployment of hydraulic fracturing wells has led to invariably declining heterogeneous production rates (4–6). Early stages of these rates follow a square root of time drop, which in the long term becomes an exponential decay leading to limited production. The production of hydraulic fracturing wells depends on three factors: the capacity of creating fracture networks successfully reaching the confined hydrocarbons, the ability of maintaining these fractures opened, and the fluid transport properties of the materials storing hydrocarbons (4–7). At the molecular scale, hydrocarbons are generated within organic agglomerates—typically referred to as kerogens—composed of an amorphous porous carbon skeleton and dispersed in a fine-grained mineral matrix (7). The chemical composition of the carbon skeleton (carbon hybridization, functional groups) is known to evolve through thermal maturation induced during geological burial (8, 9). It is therefore imperative to establish a deep scientific understanding of the structural, mechanical, and confined fluid transport properties of these organic agglomerates with respect to their thermal maturation. This type of knowledge would allow predictions of the displacement of confined hydrocarbons from the organic matter to fracture networks and their contribution in the production of hydraulic fracturing wells.

Realistic atomistic-scale models of kerogens have begun being utilized to simulate hybridization of the carbon skeletons upon maturation (8). These chemical transitions strengthen the kerogen matrix and shift the rupture mechanism from a plastic to the brittle state. In addition, the structural arrangement of nanopores is responsible for their strong adsorption, the breakdown of continuum hydrodynamics, and the occurrence of interfacial wetting effects in agreement with the fast productivity declines of exploitation wells (10, 11). Such nanostructures demonstrate a self-diffusion governed transport regime (12). These recent breakthroughs, however, only account for the subnanopores network reached through atomistic-scale models (8). The presence of hierarchical structures of organic hosted pores evolving through thermal maturation, as probed in scanning electron microscopy (5, 13), small-angle neutron scattering (14–17), and adsorption (3, 14, 18), may affect the overall properties. For example, at an experimental macroscopic scale, the organic porosity has been found to be responsible for a significant reduction of the kerogen particle modulus with respect to thermal maturation (19). In addition, the analysis of gas adsorption isotherms demonstrates that thermal maturity of the organic matter affects the surface area, pore volume, and geometrical permeability (3, 14, 18, 20). However, any interpretation of adsorption data is restricted by the simplistic assumptions adopted in the models. The issue thus becomes the validity of results from both theoretical and experimental studies, which cannot be reliable without an accurate determination of geometrical arrangements of the organic pore networks at relevant scales (21).

Recent developments of bright-field electron tomography allow imaging pores of amorphous materials with sizes above 0.5 nm in three dimensions (21, 22). Such resolution overlaps with atomistic-scale models and permits a consistent multiscale approach. In this paper, we provide 3D reconstructions of the organic pore network of one oil-prone and two gas-prone source rocks with subnanometer resolution. We show that the parameters extracted from the analysis of the tomograms compare reasonably well to those measured with N2 adsorption isotherms. This approach therefore offers an insight on the evolution of the organic pore networks in terms of structural, mechanical, and fluid transport properties at different thermal maturities. We demonstrate that the structural evolution affects not only an average pore radius and geometric surface area, but more importantly the tortuosity and connectivity. Through imaging informed simulation procedure, we further reveal that these topological features impact both mechanical and transport properties. We show at the mesoscale a maturation-induced decrease of rigidity. This offset the strengthening of carbon matrix as noted through reactive potential atomistic-scale models (8). In addition, following the previous work in this field (23–25), we study the influence of pressurized fluid inside pore networks, which leads to different stress distributions that impact crack nucleation areas. Combined with the recent findings on hydrocarbon transport within atomistic simulations of nanostructures (10–12), our analysis allows us to capture the effect of sorption on diffusion mechanisms with just a single input of the porous network geometry. As a conclusion, we show how the diffusion, occurring in the connected nanopores, constitutes the limiting factor of confined hydrocarbon transport, despite the consideration of sorption. The multiscale simulation approach presented here hence provides a way to access the interplay between topology and physical properties of porous solids at the mesoscale.

Results

We first performed adsorption isotherms on three samples from economically valuable source rock formations containing type II kerogens: Marcellus (MAR), Haynesville (HAY), and Lower Eagle Ford (LEF). The organic geochemistry of these samples () plotted in the modified Van Krevelen diagram (26) illustrates their differences in terms of hydrocarbon generation potentials (). In brief, LEF is a source rock entering the oil window whereas HAY and MAR are source rocks in the dry gas window. The specific surface areas (A), pore volumes (V), and hydraulic radii (r) obtained from N2 adsorption analyses are summarized in . The specific surface area increases with thermal maturity, which leads to decrease in the average pore size (r = 4V/A) (). This aligns with previous studies performed on extended sets of shale gas plays (3, 14, 17). These results suggest that (i) the geothermal maturation (transition from oil-prone to gas-prone) triggers the formation of organic hosted nanopores (3, 13) and (ii) our source rocks samples correspond to the two extremes of this general tendency.

To further test this hypothesis, we performed direct observations and analysis of the pore networks. Among the organic agglomerates dispersed in source rocks, the migrated organic matter (solid bitumen and/or pyrobitumen) is recognized as the main contributor of organic hosted porosity (13). Focused ion beam (FIB) thin sections were extracted on those locations in each sample and their tridimensional structures were reconstructed through electron tomography (). Fig. 1 shows the tridimensional pore aperture maps of the electron tomograms of LEF, HAY, and MAR, respectively. Significantly different structures are readily observable, as evidenced by the decreasing ranges of aperture radii. The resulting pore-size distributions (PSDs), plotted in Fig. 1, confirm this qualitative description. LEF presents a convex PSD with pore radii in the range of 3–12.6 nm, while HAY and MAR expose concave distributions dominated by nanopores (0.5 < r < 2.0 nm). It is worth noting that our field of view precludes the observation of pore with r > 20 nm, which impacts the PSD. Nonetheless, these distributions agree with previous findings that present decrease of the average pore size () with respect to thermal maturity (3, 13, 14).

Fig. 1.

Morphological characterization of the organic mesoporous networks imaged in electron tomography. Aperture maps of the organic pore networks of (A) LEF (φmeso = 9.1%), (B) HAY (φmeso = 19.3%), and (C) MAR (φmeso = 16.9%). (D) Pore-size distributions from aperture maps calculation. (E) In-pore CLDs. (F) Specific surface areas (A) calculation from Eq. as a function of the cutoff distance (rc) of the in-pore C.L.D.

The in-pore chord length distribution (CLD) is recognized as an effective stereological tool for characterizing disordered porous media (27–29). The endpoints of frequency distribution segments appear at the pore interface, fp (r), and therefore were computed within the tomograms (Fig. 1). LEF showed a constant distribution below 1 nm, which decreases slightly in the range of 1–11 nm, followed by a sharp decay. In HAY and MAR, the distributions are constant at r < 1.5 nm and exhibit abrupt decays at larger values of r. The specific surface area, A, can be extracted from the first moment of the in-pore CLD using the following equation (28):where φ is the porosity of the tomogram, ρ is the density of the solid phase, and <ℓ> is the first moment of the normalized CLD. The porosity (φ) was measured by the ratio of pore volume with the total volume of the tomogram. The density of the amorphous carbon matrix (ρ) was arbitrarily estimated as ρ = 1.0 g/cm3 for LEF and ρ = 1.4 g/cm3 for HAY and MAR (8). If A is plotted as a function of the shortest chord length r (Fig. 1), then the geometrical surface area (S) can be defined by extrapolating r to zero (29) using Eq. :where the constant α equals 1.15, 0.55, and 0.95 nm−1 for LEF, HAY, and MAR, respectively. From this expression, we obtain S = 113 m2/g for LEF, 217 m2/g for HAY, and 338 m2/g for MAR, which is in a reasonable agreement with values measured in gas adsorption experiments (). Small discrepancies between specific surface areas obtained from adsorption isotherms and geometrically determined values were previously documented for mesoporous materials and evaluated at ∼20% (30), with S being greater than A. Despite the differences between two approaches, the results confirm the increase of accessible surface area as a function of thermal maturity and the predominance of subnanopores in MAR and HAY with respect to LEF. This decrease in average pore size can be reasonably related to the development of subnanopores within the organic agglomerate as it is converted to mobile alkanes such as methane (3).

One advantage of the 3D porous network reconstructions is the ability to directly measure their topological properties. The tortuosities calculated from random walkers traveling through pore networks were 1.9 ± 0.7, 3.7 ± 0.5, and 6.6 ± 0.7 for LEF, HAY, and MAR, respectively (). These mean values increase with thermal maturity, whereas their SDs are related to the variance in PSD (31). Other studies have previously demonstrated that the geometrical permeability, k, (k = φr/ξ) of organic porous networks increases with maturity (3) [where ξ is the Kozeny dimensionless constant determined from flow core plug experiments (32)]. The evolution of permeability was interpreted as the result of a significant decrease in tortuosity (reduction of ξ by a factor of 4–8). Contradictory to previous findings, we show that the tortuosity measured on the organic porous network of three source rocks increases with thermal maturity (). In parallel, the connectivity density (ζ), obtained with the ratio of the first Betti number and the total volume of the tomogram (33), increase with maturity (). This trend also rationalizes the differences observed between S and A at low maturity (LEF), where a significant part of the pore network is disconnected and thus inaccessible to N2 molecules. This indicates that the thermal evolution of the organic matter not only entails morphological changes of the porous network (decreasing of r and concomitant increasing of S) as shown previously (3, 14), but it also implies major topological changes toward increasingly connected networks. We emphasize that these topological features impact the porosity detected in transmission electron microscopy, including organic nanopores (r < 1 nm) that mainly connect the organic mesopores (r > 1 nm).

The next step is to utilize a numerical approach to study the mechanical behavior and the confined hydrocarbons transport at the mesoscale. Lattice models were created by converting voxels of the binarized tomograms to lattice nodes corresponding to void, if porosity is detected, and the organic carbon skeleton (8). The mechanical, transport, and thermodynamical properties of the organic carbon skeleton presented hereafter take into account results from molecular dynamics using realistic numerical models of kerogen (8).

We used a fast Fourier transfrom (FFT)-based numerical method (34–36) to compute the elastic response of the lattice models from LEF and MAR (). We arbitrarily assigned to the organic carbon skeleton a local Young’s modulus of 2 GPa for LEF and 15 GPa for MAR with a Poisson’s ratio of 0.25 for both, in agreement with published data from molecular simulations on oil-prone (37) and gas-prone kerogens (8). Assuming the same Poisson’s ratio at the box scale, the global Young’s moduli resulting from the simulations are 1.4 GPa and 6.9 GPa for LEF and MAR, respectively. Considering the linear elasticity, the ratio between macro- and micro-Young’s moduli (i.e., the decrease of rigidity) constitutes the main relevant quantity. It comes to 70% for LEF and down to 46% for MAR. Since the porosity alone cannot explain this difference, we argue that the mesopore’s topology is responsible for such decreases. The presence of a homogeneous and highly disordered porous network (Fig. 1) therefore explains the counterintuitive link between the rigid organic carbon skeleton (8) of the gas-prone sample (MAR) and its flexible behavior at the upper scale. On the contrary, the oil-prone sample (LEF) presents a lower decrease in rigidity at the macroscopic scale due to heterogeneous distribution of pores, including isolated mesopores (Fig. 1). This decrease of rigidity with thermal maturity is confirmed by experimental indentations (19).

In Fig. 2, we illustrate how the elastic energy is stored inside two mesostructures under agent (a), corresponding to an imposed averaged strain = 0.01, and agent (b), including an adsorbed fluid characterized by a fluid pressure P. The total elastic energy, reduced per volume unit (106 J/m3), is 0.82 for LEF-(a), 0.9 for LEF-(b), 4.13 for MAR-(a), and 4.5 for MAR-(b). For (a) and (b) we observe a similar total stored energy for a given averaged strain. However, agent (a) leads to Gaussian energy distributions, whereas agent (b) leads to exponential energy distributions (Fig. 2). These two types of distributions have already been established in porous networks where it was shown that exponential distributions govern the crack initiation zones (38). This demonstrates that the pressurized fluid, not studied by classical tensile tests, is critical in predicting the mechanical damage of organic matter contained in source rocks. In fact, hydraulic fracture propagation is directly influenced by stress distribution as crack nucleation appears more likely in high elastic energy areas. It corresponds to the vicinity of the largest pores for LEF-(b), leaving the rest of the matter undamaged. Contrary, crack nucleation can be more homogeneous at higher thermal maturity [see Fig. 2, MAR-(b)].

Fig. 2.

Elastic simulations of the two porous mesostructures originating from LEF and MAR under two types of agent: (a) imposed strain boundary, and (b) adsorbed fluid with free strain boundary. The averaged strain of the simulation box is = 0.01 for all cases. The size is 700 × 700 pixels with a resolution of 0.42 nm for LEF and 0.35 nm for MAR. Mesopore phase is in white. (Top) Elastic energy density fields (106 J/m3). (Bottom) Distribution of the elastic energy density (106 J/m3) for each simulation.

We applied the continuous time random walk (CTRW) (39–42) algorithm to study the confined hydrocarbons transport properties over the geometry given by the mesoscale lattice models. The numerical scheme used here to simulate transport consists of constructing random walks, where walkers can move isotropically on the lattice. Because two different diffusion timescales are considered, one for each phase (i.e., the nanoporous and the mesoporous phases), we also assign random transit times at each step of the random walk from exponential distributions characterized by an average transit time that depends on the phase where the walker is located. To account for sorption between nano- and mesoporous phases, we use partial boundary conditions at the interface (see for more details on the homogenization procedure). Due to limited size of the mesopores and their connectivity through the nanoporous phase, where transport has been shown to be purely diffusive (10, 12, 42), we neglect advection in the system. Moreover, the amorphous nature of the nanoporous phase would prevent hydrodynamic slip at the interface. As a result, we treat both phases as homogeneous, which means that we attribute a single fluid concentration ρ and diffusion coefficient D per phase. Fig. 3 shows the homogenized diffusion coefficient, , adimensionned by as function of the diffusion ratio in the MAR case for various concentration ratios . We observe that the upscaled diffusion coefficient increases with the diffusion ratio. We use an effective medium theory (EMT) (43, 44) to rationalize the dependence on aswhere is the ratio of the volume fractions of the mesoporous phase (φ) and the nanoporous phase , and is an obstruction factor that captures the impact of geometry and sorption at the interface on the homogenized diffusion coefficient. Fig. 3 reports an excellent agreement between numerical results and the EMT where is taken as for each curve. We notice that the influence of sorption depends nonmonotonically on the concentration ratio that controls the sorption effect.

Fig. 3.

Homogenized fluid transport properties of the mesoscale lattice models. (A) Homogenized diffusion coefficient of the MAR reconstruction scaled by as function of the diffusion contrast for different concentration ratios indicated by colors. The corresponding curves stand for the effective medium theory (Eq. ) where the value of is when (i.e., ). (B) Concentration-dependent obstruction factor as function of the probability for a random walker to cross the interface from the nano- to the mesoporous phase (red: LEF, green: HAY, blue: MAR). The corresponding curves display the analytical model (Eq. ) where is simply taken as . (C) Homogenized diffusion coefficient of methane for LEF (red) and MAR (blue) as function of the bulk fluid pressure. The dashed curves correspond to the diffusion in the nanoporous phase only, while the solid lines represent the homogenized model with mesoporosity. Standard geological conditions (T = 400 K and lithostatic pressure of 25 MPa) have been considered.

Fig. 3 displays, for LEF, HAY, and MAR, an evolution of the obstruction factor (γ2) as a function of, , the probability of a random walker crossing an interface from nano- to mesoporous phase. When P = 0 (or = 0) there is no fluid in the mesopores; is a constant that corresponds to the obstruction factor of the nanoporous phase, while the mesopores are considered as an impermeable phase. In practice, this can be computed by a simple random-walk simulation with bounce-back rules at the interface. If the concentration in mesopores increases and the probability p stays below 0.5, the diffusion increases due to an increased probability of the fluid entering mesopores. In this regime the fluid can use the mesopores as shortcuts, thus lowering the influence of the tortuosity of the nanoporous phase on the overall diffusion. When p gets closer to 0.5, the efficiency of the shortcut effect reaches a maximum because the decreasing probability of exiting from a mesopore balances the advantage of entering it. At , there is no influence of the interface because the probability of crossing the interface in both directions is the same, thus resulting in . Beyond this point (), mesopores start to act as traps because the probability of entering a mesopore becomes larger than the probability of escaping it. The homogenized diffusion coefficient thus diminishes and ultimately goes to zero in the limit of perfect trapping ( or ). An analytical model () for the concentration-dependent obstruction factor becomesFig. 3 shows an excellent agreement between the proposed model and numerical results. For three very different geometries, the model (Eq. ) captures the effect of sorption on diffusion through p with only one parameter , which depends on the geometry of the system.

We then compare, in Fig. 3, the diffusion coefficient in a single nanoporous phase with the one accounting for the distribution of mesopores given by the tomograms as a function of pressure P of bulk methane in mesopores. For the diffusion coefficient in the nanoporous phase, we use a free volume theory. The methane concentration in the nanoporosity is given by a Langmuir adsorption isotherm. Transport and adsorption behaviors in subnanopores for MAR were obtained by molecular simulation studies using the atomistic-scale model of Marcellus (8, 45). As a first approximation, we treat methane in mesopores as a bulk fluid. Bulk properties of methane, diffusion coefficient, and concentration as a function of P are taken from the National Institute of Standards and Technology database. For the MAR case, diffusion at low pressure is dominated by the nanoporous phase because of the strong heat of adsorption of nanopores due to confinement effects. Then, instead of decreasing exponentially with pressure, the diffusion model that accounts for mesopores reaches a minimum around 15 MPa and increases at higher pressure to become almost constant when P approaches 100 MPa (Fig. 3). Consequently, the recovery of confined hydrocarbons is impacted at the mesoscopic level by nonmonotonic evolution of transport properties in time. We also report the same quantities in the LEF case but with diffusion in the nanoporous phase described by a more general free-volume model that accounts for adsorption-induced swelling (), which leads to an exponential increase as a function of pressure (46, 47). Due to swelling, the adsorption isotherm evolves linearly with P (until P > 50 MPa) (48), thus lowering the impact of sorption on transport at the mesoscale. As a result, mesopores do not induce a qualitative change on the transport behavior at the mesoscale. Quantitatively, in both cases, the application of the CRTW algorithm showed that the transport of confined hydrocarbons at the mesoscale is mainly dictated by the nanoporosity. If mesopores in the gas-prone (rigid) case can induce nonmonotonical trend for the diffusion coefficient, they do not induce that effect for the oil-prone case (flexible nanoporosity). At the engineering scale, these trends complement the effect of fractures on wells’ productivity. As already pointed out in previous studies (49–51), the pressure dependence of fractures’ permeability can lead to a decrease of productivity when pressure decreases. The diffusion mechanisms evidenced here would either accentuate the productivity decrease for oil-prone organic matter, or partially mitigate it for gas-prone organic matter.

Discussion

Bright-field electron tomography was applied on source rocks’ organic matter, which allowed 3D reconstructions of the organic porous networks with a subnanometer resolution. Structural analysis of the tomograms demonstrated that thermal maturation entails increase of the accessible surface area, tortuosity, and connectivity of the nanoporous networks. Lattice models were then produced using the coupling between (i) geometrical arrangements of the porous networks at the mesoscale with (ii) nanoscale features deduced from a molecular dynamics study (8). These models allowed studying numerically an elastic response and the confined hydrocarbons transport properties at the mesoscale. Our results highlight the effect of topological changes on the organic matter rigidity. Moreover, we showed that pressurized fluid inside pore networks implies different potential crack nucleation areas, localized around the largest pores for the oil-prone case (LEF), and more homogeneously distributed at higher thermal maturity (MAR). The dominant impact of diffusion mechanisms occurring within the nanoporosity on confined hydrocarbon transport, previously sensed through atomistic studies, was also confirmed at the mesoscale. The constraint from imaging allowed a relatively simple computation of the only two parameters required by the fluid transport model: the volume fraction of the mesoporosity and the obstruction factor of the nanoporous phase. These results have also led to some questions, specifically whether the evidenced mesoscale effects may be generalized to a larger ensemble of source rocks, and whether one can identify the key geochemical parameters of the thermal history triggering the evolution. To conclude, the herein constructed analytical model can be readily used for predictions of “in-place” confined hydrocarbon recovery in multiscale approaches accounting for macropores, organic/inorganic interfaces, and fracture networks (52–54). This further improves the link between nanoscale constitutive models of porous materials (8, 10–12) and their engineering applications (4–6).

Materials and Methods

This study was conducted on samples from economically valuable source rock formations (MAR, HAY, and LEF). The thermal maturity was evaluated by RockEval pyrolysis and standard organic petrology technique. Low-pressure N2 isotherms were performed at 77 K using a Micromeritics 3Flex. Thin sections of organic inclusions were extracted in an FEI Helios NanoLab 660 FIB. Electron tomography acquisitions were carried out with a transmission electron microscope (JEOL JEM 2011) fitted with a LaB6 electron gun, under a 200-kV accelerating voltage, and assisted by the Digital Micrograph software (GATAN). The structural characterization was performed using iMorph (imorph.sourceforge.net) and imageJ. Lattice models were then created by converting voxels of the tomograms to lattice nodes. An FFT-based numerical model was used to compute the elastic response under two different agents: (i) an imposed averaged strain, and (ii) an adsorbed fluid inside pores. The CTRW algorithm was finally used to study the transport properties over the geometry of the lattice models (see for complementary details on ).