Layering in Shales Controls Microfracturing at the Onset of Primary Migration in Source Rocks

The process of primary migration, which controls the transfer of hydrocarbons from source to reservoir rocks, necessitates the existence of fluid pathways in low permeability sedimentary formations. Primary migration starts with the maturation of organic matter which produces fluids that increase the effective stress locally. The interactions between local fluid production, microfracturing, stress conditions, and transport remain difficult to apprehend in shale source rocks. Here, we analyze these interactions using a coupled hydro‐mechanical numerical model based on the discrete element method. The model is used to simulate the effects of fluid production emanating from kerogen patches contained within a shale rock alternating kerogen‐poor and kerogen‐rich layers. We identify two microfracturing mechanisms that control fluid migration: (i) propagation of hydraulically driven fractures induced by kerogen maturation in kerogen‐rich layers and (ii) compression‐induced fracturing in kerogen‐poor layers caused by fluid overpressurization of the surrounding kerogen‐rich layers. The relative importance of these two mechanisms is discussed considering different elastic properties contrasts between the shale layers, as well as various stress conditions encountered in sedimentary basins, from normal to reverse faulting regimes. The layering in shales causes local stress redistribution that controls the prevalence of each mechanism over the other and the onset of microfracturing during kerogen maturation. The approach is applied to the Draupne formation, a major source rock in the Norwegian continental shelf in the North Sea.


Introduction
In several geological systems, fluid may be expelled from low permeability rocks. This process happens during dehydration of sediments and serpentine rocks in subduction zones, dehydration of clay minerals in sedimentary basins during diagenesis, or hydrocarbon expulsion from source rocks, a process called primary migration. Shale rocks deposited in sedimentary basins can act as source rocks of hydrocarbon-bearing reservoirs, as well as cap rocks for the same reservoirs, or for underground storage, for example, Carbon Capture Storage. Shales prevent fluids to escape due to their low permeability and to a capillary sealing mechanism controlled by small pores (Horsrud et al., 1998). Shale rocks, when acting as source rocks, are of primary interest for hydrocarbon exploration. However, the apparent low permeability of shales questions the mechanisms that lead to hydrocarbon expulsion. Some shales collected in outcrops or cores samples contain microfractures and microporosity at small scales that can increase the overall permeability of the rock (Gale et al., 2014;Kalani et al., 2015;Ougier-Simonin et al., 2016). However the mechanisms at the origin of these microfractures remain unclear.
Kerogen maturation into hydrocarbons leads to overpressure at the microscale that can induce fracture propagation and therefore fluid migration. This mechanism has been proposed to explain shale microfracturing (Fan et al., 2012;Jin et al., 2010;Pelet & Tissot, 1971;Vernik, 1994) and has been numerically investigated using a coupled hydro-mechanical model by Teixeira et al. (2017). These authors showed that kerogen maturation induces the formation of hydraulically driven microfractures that tend to propagate toward the largest principal stress direction within an initially homogeneous and isotropic material. Under normal faulting conditions encountered in many sedimentary basins, the largest principal stress direction being vertical, fractures would predominantly open along this direction, thus facilitating vertical fluid migration. Left: Schematic shale rock structure before maturation. Right: Microfracturing process during kerogen maturation. The increase of fluid pressure in kerogen patches produces microfractures through two mechanisms. Mechanism 1: nucleation, propagation, and rotation of hydraulically driven microfractures emanating from kerogen patches (in black) in layers with high kerogen content (dark red). Mechanism 2: compression induced microfractures (in white) in layers with low kerogen content (light red). Both mechanisms increase the permeability of the medium and may lead to the formation of connected fluid pathways, allowing fluid expulsion from a tight rock.
Nonetheless, shale rocks are complex materials due to their formation process that involves the successive deposition of millimeter to centimeter thick layers with variable mineral and organic contents. These rocks are inherently anisotropic materials (Valès et al., 2004) that may show fractal distribution of layer thickness (Li et al., 2017). Such structural layering can lead to local stress inversion between successive sedimentary units within sedimentary basins (Gunzburger & Magnenet, 2014). Depending on the deposition environment and on the regional tectonics, local stress conditions can vary from normal faulting stress regime when the vertical stress v is larger than the horizontal stress h , to reverse faulting stress regime when the horizontal stress is larger than the vertical stress (Zoback, 2007, Chapter 9;Zhang & and Zhang, 2017). Obviously, the interactions between the far-field stress conditions and the mechanical properties of the materials play a major role on the deformation processes taking place in shale formations.
We develop here a numerical study to investigate the influence of shale layering on the microfracturing processes occurring during kerogen maturation and the effect of these processes on permeability. Two microfracturing mechanisms are proposed and discussed as candidates to increase the permeability ( Figure 1). The first mechanism is the propagation of microfractures emanating from kerogen patches due to the local increase of fluid pressure generated by kerogen maturation. The second mechanism is the ability of this overpressure to damage the surrounding layers that do not contain kerogen. The first mechanism has been discussed in previous studies (e.g., Jin et al., 2010, and references therein), while, to our knowledge, the second mechanism has not been proposed yet to explain primary migration in source rocks.
The main goal of the present study is thus to quantify the importance of these two different microfracturing mechanisms induced by kerogen maturation in layered rocks. In particular, the study focuses on the effect of the mechanical contrast (Young's modulus) between layers on the formation of microfractures considering three different mechanical loading conditions: a hydrostatic case ( v = h ), a normal faulting case (i.e., extension, v > h ), and a reverse faulting case (i.e., compression, v < h ). Finally, we apply our findings to explain primary migration in the Draupne shale formation, a major source rock in the Norwegian continental shelf.

The Draupne Shale Formation, Norwegian Continental Shelf
The Upper Jurassic (Kimmeridgian) Draupne shale formation is found throughout the Norwegian Continental Shelf (NCS) and is the main source, carrier, or seal rock in several oil and gas fields. It is a world-class source rock in the Central and Northern North Sea (Figure 2a). This formation stretches over a broad range of depths (500-6,500 m), thicknesses (0-550 m) (Figures 2c and 2d), and maturity. It provides an ideal candidate to study microfracturing due to maturation of organic matter and expulsion of hydrocarbon. The Draupne formation has been substantially indurated in the North Sea by burial and later uplifted with the maximum uplift documented towards the Norwegian coast. In the Southern North Sea, the equivalents to the Draupne formation are the Mandal and Tau formations, which vary in thickness from 1 to 170 m (Mandal formation) and 2-118 m (Tau formation), respectively. In the Norwegian Sea and the Barents Sea, the equivalents to the Draupne formation are the Spekk and Hekkingen formations, which vary in thickness from 1 to 150 m (Spekk formation) and 2 to 190 m (Hekkingen formation), respectively (NPD, 2020).
The Draupne formation is made of dark, gray-brown to black, usually non-calcareous, carbonaceous, occasionally fissile shale (Figure 2b) deposited in a marine environment with restricted bottom circulation and often with anaerobic conditions. Minor limestone streaks and concretions occur throughout the formation. It has low elastic wave velocity (V p = [2, 800 − 3, 600] m . s −1 , V s = [1, 000 − 1, 400] m . s −1 ), low density (2.24 g·cm −3 ), and high resistivity due to high organic content (TOC in the range 5-12 wt%) (Mondol, 2019). Several Draupne shale samples have been collected at depth in wells. For example, the well 16/8-3S was drilled in the Ling Depression, Central North Sea, where the thickness of the formation is 86 m. The core is naturally split into subsections along with bedding parallel fractures. This shale has the following composition: total organic content (6.8%), clays (51.1%), quartz (22.4%), microcline (17.8%), pyrite (7.7%), dolomite (1.2%), and calcite (0.7%). It is layered at the millimeter scale (Figure 2e), with alternations of kerogen-rich and kerogen-poor layers. In the following, we consider this shale representative of the Draupne formation to perform numerical simulations of kerogen maturation-induced microfracturing.

Numerical Modeling of Shale Maturation
To investigate microfracturing induced by kerogen maturation in shale rocks, we use the coupled hydro-mechanical model proposed by Papachristos et al. (2017). The model couples the discrete element method (DEM) and the pore-scale finite volume method (FVM) implemented in the open source platform YADE DEM (Kozicki & Donzé, 2008, 2009Smilauer et al., 2015). The same approach has been previously used to investigate kerogen maturation in homogeneous isotropic rock samples by Teixeira et al. (2017).
The rock material is modeled as an assembly of bonded rigid particles whose respective displacement is ruled by Newton's second law of motion ( Figure 3). The bonds binding each pair of interacting particles obey an elastic-brittle behavior in directions both normal and tangential to the contact plane. Under stress or strain loading, these bonds can break by either tensile or shear failure mechanism, following a modified Mohr-Coulomb criterion with a tension cut-off (Scholtès & Donzé, 2013).
When an interparticle bond breaks, a microfracture is explicitly formed in the material. Each microfracture has a surface area proportional to the radii of the interacting particles formerly bonded (i.e., ( R 1 +R 2 2 ) 2 , with R 1 and R 2 the radii of the particles), and its normal vector is collinear to the vector joining the two particles centroids. When several microfractures form within the solid, they can coalesce to form fractures. The aperture of the microfractures can increase or decrease during the deformation of the solid and is calculated from the normal relative displacement between the particles.
Besides, the poral space is discretized through a Delaunay triangulation that allows simulating a compressible fluid flowing between the solid particles . Because the present study focuses on the response of very low permeability rocks to kerogen maturation, we assume that flow only occurs within microfractures that preexisted in the rock or were induced by the loading and the maturation.
The flow along each microfracture is modeled as a flow between two parallel plates with constant hydraulic aperture and permeability described by the cubic law (Witherspoon et al., 1980). The hydraulic aperture corresponds to the aperture of the microfracture measured from the particle displacement. The method is hydro-mechanically coupled because (i) local fluid overpressure induces forces on the solid particles that eventually generate hydraulically driven microfractures as a result of interparticle bond breakage; and (ii) a deformation of the solid skeleton induces local fluid overpressures and associated fluid flow within the microfracture network.

Numerical Simulation Setup 2.3.1. Model Geometry
The simulated rock specimen is made of 50,000 spherical particles packed in a volume with dimensions 14 × 10 × 10cm 3 . The medium is divided into seven horizontal layers, alternating material with high organic content (#1, red) and material with low organic content (#2, blue), as shown in Figure 4. The layers with high organic content contain kerogen patches, and the layers with low organic content do not. The interfaces between the layers have the same cohesion as the layers.
The kerogen patches within high organic content layers are modeled as penny-shaped microfractures with a radius of 1 cm. Each microfracture extends over 6.5 discrete elements in the bedding direction of the shale material to reproduce the bedding-parallel preferential orientation of kerogen patches observed using synchrotron X-ray tomography imaging (Kobchenko et al., 2011;Panahi et al., 2013Panahi et al., , 2018. The number of elements has been chosen as a compromise between a fine discretization of the layers, with 6.5 discrete elements across the layer thickness, and a reasonable simulation time.

Stress Boundary Conditions
In the following, we use the engineering stress convention where compressive stresses are negative and extensional stresses are positive. External (boundary) stresses are applied on the numerical specimens through frictionless rigid walls whose movements are servo-controlled during the simulations. When the specimens contain several layers presenting mismatching mechanical (elastic) properties, the walls are subdivided into smaller units that are all independently servo-controlled in order to avoid stress concentrations at the layers boundaries. For all the configurations tested here, the loading conditions are chosen such that the response of the simulated medium remains in a quasi-static regime. Therefore, the results do not depend on the loading rate.
Kerogen patches are inserted in the numerical medium at the initial stage of the simulations as closed pre-existing penny-shaped microfractures ( Figure 4). These patches correspond to interparticle interfaces where the cohesive strength is set to zero. These interfaces behave in an elastic-frictional manner following the smooth-contact approach proposed in (Scholtès & Donzé, 2012). They all share the same orientation, parallel to the bedding, and are hydraulically connected one with another.
The particle assembly is first compressed to the desired state of stress by moving the surrounding walls at a constant strain rate without considering any hydro-mechanical coupling. Once the predefined state of stress and static equilibrium are reached, the stresses at the boundaries are kept constant by adjusting the position of the surrounding walls.
Three sets of numerical experiments were designed to simulate three different stress states: hydrostatic con- Our simulations are performed for stress states relevant to shale maturation processes and particularly to the Draupne formation in the North Sea where maturation occurred at a depth close to 2,570 m ( v = −55MPa, h = −36MPa, and pore pressure P P = 25MPa) (Zadeh et al., 2017). The reverse faulting and hydrostatic stress conditions are defined such that the first invariant of the stress tensor (i.e., the pressure, I 1 = trace(̄)) is the same for all simulations.

Kerogen Maturation, Overpressure, and Fluid Flow
Black shales deposited in basins commonly contain a significant volume fraction of kerogen that transforms to hydrocarbons as temperature and pressure increase during burial (Luo & Vasseur, 1996). We simulate kerogen maturation by injecting a fluid at a constant flow rate in the kerogen patches, as proposed by Teixeira et al. (2017). Injecting fluid in the patch mimics the pressure buildup that occurs during the maturation of organic matter without considering the full chemo-thermo-hydro-mechanical processes at stake.
The volume expansion associated with conversion of kerogen into hydrocarbon results in a significant increase in fluid pressure, which leads to microfracturing around the kerogen patches . Here, we injected a Newtonian fluid with a bulk modulus of 2.210 9 Pa and a viscosity of 10 3 Pa·s. The injection generates overpressure in the kerogen patches that eventually leads to the formation of microfractures in its vicinity as a result of interparticle bond breakage. The fluid can then flow within these newly created microfractures to propagate further into the rock matrix, creating preferential fluid pathways.
The elastic waves induced by bond breakages are damped to avoid dynamic effects, such that simulations are quasi-static. The injection rate and the simulation time are not directly related to the geological time under which maturation takes place. Hence, the rate of fluid production in the kerogen patches defines the time scale, and therefore, no absolute time is given when presenting the results of the numerical simulations in the following.

Material Properties
Each layer is modeled as an isotropic and impermeable rock material. The mechanical properties are chosen from values measured at depth in the Draupne shale formation (Mondol, 2019). The Young's modulus and the Poisson's ratio are different for each layer whereas the uniaxial tensile strength is identical. In addition, because the migration of fluid is driven by the propagation of overpressurized kerogen-filled penny-shaped patches, we estimated the mode I fracture toughness, K IC , of the numerical specimen by running a dedicated numerical test where a single penny-shaped microfracture was placed in a homogeneous sample subjected to a tensile uniaxial external loading. We calculated a value K IC = 0.82MPa·m −1/2 . The mode I fracture toughness of Draupne shale is not known. However, our calculated value is comparable to the K IC measured by Chandler et al. (2016) on Mancos shale, which varies in the range [0.21, 0.75]MPa·m −1/2 depending on the orientation of the bedding relative to the loading direction. The properties of the high and low organic content layers are summarized in Table 1.

State of Stress in a Layered Elastic Material
The state of stress in a stratified material formed by two series of isotropic layers with different elastic properties can be obtained from Hooke's law: The local fluid production in kerogen patches may fracture the rock. Because the largest principal stress direction controls the direction of hydrofracture propagation (Teixeira et al., 2017), the ratio hi ∕ v between the horizontal stress and the vertical stress in the layer, #i, plays a crucial role. If hi ∕ v > 1, the largest principal stress direction is horizontal (i.e., parallel to bedding), and if hi ∕ v < 1, the principal stress direction is vertical (i.e., perpendicular to bedding). The expression of this ratio is given in equations (2) and (3) for a rock made up of two different materials, where E i and nu i are the elastic parameters of material #i, P = h ∕ v is the ratio between the horizontal stress and the vertical stress imposed on sample boundaries, and q1 is the volume proportion of material #1. A derivation of those expressions is given in Appendix A.
Figure 5 presents a comparison between the local stress distributions calculated from the analytical solutions (equations (2) and (3)) and the local stress distributions within the layered DEM model shown in Figure 4 for hydrostatic loading conditions ( v = h = −17.3MPa). Figure 5 shows the stress distribution within the rock medium under mechanical equilibrium, before the fluid was injected into the kerogen patches. At  (2) and (3), and numerical local stresses predicted by the DEM model presented in Figure 4 (q1 = 3∕7), considering a hydrostatic state of stress (P = 1 with v = h = −17.3MPa). The numerical local stresses are computed by averaging the stresses on particles over a layer of thickness x along the vertical axis. These results are obtained for two sets of high organic content and low organic content rock parameters: (a) Material #1, E 1 = 6.8GPa, this stage of the simulation, no overpressure due to kerogen maturated is generated, and the response of the numerical model is mostly elastic as no or very few microfractures are induced by the mechanical loading.
The model is in good agreement with the analytical solutions for both the local vertical stress v and the local horizontal stress h . The results also show that, depending on the mechanical contrast between the layers, the orientation of the largest principal stress can be either vertical or horizontal in the different shale layers, as observed in some log measurements in the Paris basin (Gunzburger & Magnenet, 2014).

Results for a Nonlayered Rock
The way microfractures propagate directly depends on the mechanical properties of the rock and the state of stress. Here, we study the effect of each property on maturation-induced microfracturing in a nonlayered medium without considering the structural effects caused by the layering. The case of a layered rock is studied in section 4.

Hydraulically Driven Microfractures in Layers With High Organic Content
We simulate here the maturation of a sample containing a unique kerogen patch subjected to different far field stress regimes ( Figure 6). The four stress conditions are defined keeping the same value of the first invariant of the stress tensor I 1 , and different differential stresses equal to 0,9,11.5,and 19GPa. Under hydrostatic loading (zero differential stress), the hydraulically induced fractures are preferentially oriented in the x0z plane, parallel to the bedding, as shown on the fluid pressure map as well as on the stereographic projection of the microfractures orientation distribution. In this case, the direction of microfracture propagation is driven by the orientation of the initial kerogen patch.
For the simulation with the highest differential stress (19GPa), the hydraulically induced fractures are preferentially oriented in the yOz vertical plane, perpendicular to the bedding, as confirmed by the stereographic projection of the microfracture planes that shows a maximum in the x direction. This result confirms that fluid-induced microfractures propagate and rotate toward the direction of the largest principal stress direction, along the y-axis in this case.
Both the stereographic plots and the fluid pressure maps show that an increase of the differential stress in the vertical direction enhances the number of vertical fractures, an effect that can be interpreted by a faster reorientation of the fractures and the creation of multiple vertical fluid pathways, as shown by the differences between the −25MPa and −30 MPa cases on Figure 6. To conclude, for a given set of material properties, an increase of differential stress enhances microfracture propagation toward the principal stress direction.

Overpressure Due to Kerogen Maturation in Layers With High Organic Content
The transformation of kerogen into hydrocarbon changes the stress state by two processes. First, dilation of the kerogen patches produces a local overpressurization in the surrounding solid. Second, the propagation of hydrocarbon-filled microfractures modifies the local stress field in the vicinity of the kerogen patches. These two effects combined together increase the differential stress in the kerogen-poor layers located below and above the kerogen-rich layers. This process leads to the formation of compression induced microfractures in the kerogen-poor layers (Mechanism 2 in Figure 1).
In order to assess the effects of material properties and external loading conditions on this mechanism, we performed simulations with a single kerogen patch contained in a homogeneous layer whose Young's modulus was varied between 6.8 and 9.9GPa (Table 1). Figure 7 displays the maximum fluid pressure recorded at the center of the kerogen patch for three different stress conditions: hydrostatic, normal faulting, and reverse faulting stress regimes. While the flow rate is imposed at a source point located at the center of the kerogen patch, the pressure is homogeneously distributed inside the patch before fracturing. For the range of Young's moduli tested, the maximum pressure is independent of the Young's modulus (Figure 7). The maximum pressure depends on the external stress conditions and is correlated to the minimum stress value.
To conclude, the maximum fluid pressure reached within the kerogen patches is independent of the Young's modulus in the range of parameters studied here and strongly depends on the applied external stress.

Compression-Induced Microfractures in Layers With Low Organic Content
In order to assess the dependency of damage on stress increase and material properties (Young's modulus), we performed triaxial compression test simulations on a low organic content rock (material #2), without kerogen. The initial stress conditions correspond to the normal faulting stress with v = −30 MPa and h = −11 MPa. The vertical stress was increased by increments of 1 MPa up to −40 MPa. Figure 8 presents the number of microfractures generated in the rock as a function of the vertical stress considering five simulations run with different Young's moduli ranging from 6.8 GPa to 9.9 GPa (Table 1). In the range of values investigated here, the Young's modulus has a negligible effect on the damage produced within the sample. Damage results from the increase of the vertical stress. To conclude, microfracturing in low organic content layers is controlled by the vertical stress in this layer.

Summary of the Results for a Nonlayered Rock
In this section, we summarize how two microfracturing mechanisms take place in a nonlayered rock during kerogen maturation (see Figure 1). The results show that:   (2) and (3). h = v = −17.3MPa. The total volume of injected fluid is the same for each case: 0.2mm 3 .
• (i) The propagation of fluid induced microfractures is directly influenced by the far field stress: high deviatoric stress tends to favor the propagation of hydraulically driven fractures in the direction of the largest principal stress. • (ii) The overpressure needed to propagate hydraulically driven microfractures from a kerogen patch is independent of the Young's modulus of the solid and strongly depends on the applied external stress and the stress perpendicular to the organic patch. • (iii) The number of microfractures in a layer with no organic content increases with the first stress invariant I 1 # 2 and is independent of the Young's modulus in the range of parameters tested here.

Results for a Layered Rock
We investigate three different stress conditions: hydrostatic, normal faulting, and reverse faulting regimes. For each scenario, 20 simulations were performed to consider every possible combinations of the rock material properties given in Table 1.
In each simulation, the number of microfractures induced in the layer with low organic content is computed after the same amount of fluid was injected in the kerogen patches. Results are displayed in Figures 9,11, and 13. Each figure contains two plots where the vertical and horizontal axes correspond to E 1 and E 2 , the Young's moduli of the layers with high organic content #1, and low organic content #2, respectively. The background color maps represent the ratio h1 ∕ v in the layer with high organic content (left panel, equation (2)), and the ratio h2 ∕ v in the layer with low organic content (right panel, equation (3)) because these parameters discriminate which mechanism induces the fracturing (see section 3.4). Figure 9 shows the results of simulations performed under hydrostatic stress condition ( v = h = −17.3MPa). Two different areas in the (E 1 , E 2 ) plots can be distinguished.

Hydrostatic Stress State:
First, the red shaded area in the left panel corresponds to a local stress condition where the largest principal stress is vertical in the layers with high organic content. In this configuration, because of the deviatoric stress in the high organic content layer, the fluid-induced microfractures emanating from the kerogen patches are oriented toward the vertical direction and propagate into the layers with low organic content, as shown in Figure 10a.
Second, the blue shaded area in the left panel corresponds to a local stress condition where the largest principal stress is horizontal in the layers with high organic content (#1, red layers in Figure 4). The microfractures appearing within the layers with low organic content are induced by the compression of the layer with low Scenario with E 1 = 9.9GPa in high organic content layers, and E 2 = 11.5GPa in low organic content layers. (c) Scenario with E 1 = 9.9GPa in high organic content layers, and E 2 = 8.0GPa in low organic content layers.
organic content (#2, blue layers in Figure 4) caused by the over-pressurization of the high organic content layers. As shown in Figures 10b and 10c, the hydraulically driven microfractures propagate horizontally inside the high organic content layers because the largest principal stress is horizontal. The amount of damage in the layers with low organic content is correlated to the initial confinement of the layer with low organic content (I 1,#2 ). This effect is observed on Figures 10b and 10c where the case with higher I 1,#2 corresponds to Figure 10b, which shows more fractures in the low organic content layers.
For a given value of E 1 (e.g., 10GPa), the number of microfractures increases when E 2 and I 1,#2 increase and when I 1,#1 decreases. A decrease of I 1,#1 results in a decrease of the overpressure within the kerogen patches, which should lead to less damage in the layers with low organic content. However, this effect is counterbalanced by the increase of the confinement in the layers with low organic content I 1,#2 . As a consequence, the increase of confinement in the layers with low organic content I 1,#2 has a stronger effect on the development of microfractures in those layers.
To conclude, under hydrostatic stress condition, the mechanism and degree of microfracturing in the layers with low organic content depend on the contrast of elastic properties between the layers. On the one hand, for an elastic contrast leading to a local largest principal stress oriented vertically in the layers with high organic content (red area in Figure 9, left), microfracturing of the layers with low organic content results from the reorientation of the fluid driven microfractures emanating from the kerogen patches in the direction perpendicular to the bedding. On the other hand, for an elastic contrast leading to the local largest principal stress oriented horizontally in the layers with high organic content (blue area in Figure 9, left), microfracturing of the layers with low organic content is caused by its compression due to the overpressurization of the two layers with high organic content located above and below. The parameter that has the strongest impact on this damaging process is the initial confinement of the layer with low organic content (I 1,#2 ). Figure 11 shows the results of simulations performed under normal faulting stress conditions, with stress values similar to that measured at depth in the Draupne shale formation ( v = −30MPa, h = −11MPa) (Zadeh et al., 2017).

Normal Faulting Stress State: v = −30MPa, h = −11MPa
Under such strong deviatoric state of stress (19MPa), the largest local principal stress in both layer is vertical ( Figure 11). As shown in Figure 11, the stress distribution over the range of elastic parameters investigated is similar for all cases (i.e., h1 ∕ v ∈ [0.45, 0.5], h2 ∕ v ∈ [0.26, 0.3]). The number of microfractures induced within the low organic content layers after the injection of 0.20 mm 3 of fluid in the kerogen patches is significantly higher compared to the hydrostatic case ( Figure 9). As shown in Figure 12, kerogen maturation Figure 11. Normal faulting stress conditions. Number of broken bonds (i.e., microfractures) in the layers with low organic content (#2, blue layers in Figure 4). Each circle corresponds to a simulation. The background colors show the value of h1 ∕ v in the layers with high organic content (left) and the value of h2 ∕ v in the layers with low organic content (right). These values are computed using equations (2) and (3). h = −11MPa, v = −30MPa. The total volume of injected fluid is the same for each case: 0.20mm 3 .
causes the propagation of fluid driven microfractures toward the vertical direction that corresponds here to the direction of both the local and global largest principal stresses. Hence, the main mechanism responsible for microfracturing in the layers with low organic content in this configuration is the reorientation of the fluid driven microfractures emanating from the kerogen patches. Figure 13 shows the results of simulations performed under reverse faulting stress conditions ( v = −10MPa, h = −21MPa). In this case, the largest local principal stress is horizontal in both the low and high organic content layers ( hi ∕ v > 1). Kerogen maturation induces the propagation of fluid-driven microfractures along the horizontal direction in all configurations, as shown in Figure 14.   (2) and (3). h = −21MPa, v = −10MPa. The total injected fluid volume is the same for each case: 0.2mm 3 .

Reverse Faulting Stress State: v = −10MPa, h = −21MPa
Under this reverse faulting regime, the number of microfractures in the low organic content layers increases with the ratio h2 ∕ v and therefore with the initial confinement of the low organic content layers (I 1,#2 ) as shown in Figure 13.
On the contrary, the initial confinement of the high organic content layers (I 1,#1 ), which is link with an increase of of h1 ∕ v is not correlated with the number of fracture ( Figure 13). Therefore, the difference of Figure 14. Reverse faulting stress conditions. Three-dimensional view of the microfractures induced by kerogen maturation for two different sets of Young's moduli in the layers with high organic content (red) and low organic content (blue). The yellow cells correspond to microfractures filled with fluid and the dark blue cells to microfractures that do not contain fluid. The boundary stresses are h = −21MPa and v = −10MPa, and the total injected fluid volume is 0.2mm 3 . (a) Scenario with E 1 = 6.8GPa in the high organic content layers, and E 2 = 8.0GPa in the low organic content layers. (b) Scenario with E 1 = 9.9GPa in the high organic content layers, and E 2 = 14.9GPa in the low organic content layers.
pressure increase link with the confinement increase of the layer does not have a significant impact on the fracture mechanism of the low organic content layers.
To conclude, the result indicates that the initial confinement of the low organic content layer I 1,#2 control the fracturation within this layer. These conclusions are similar to those obtained for the hydrostatic loading conditions (section 4.1).

Mechanics of Layered Shales
Materials that contain weak and strong layers have been widely studied to characterize the variety of fracturing mechanisms they experience (Bourne, 2003;Boersma et al., 2020;Douma et al., 2019;Guo et al., 2017;Teufel & Clark, 1984). These studies showed that a change in fracture behavior in vertically compressed layered samples can be attributed to the nature of the horizontal stress distribution within each layer that depends on the elastic contrast between the layers. Mode I fractures develop in layers under extensive horizontal stress while mode II fractures develop in layers under horizontal compressive stress. The strength of the interface between the layers controls the transmission of stress across them (Bourne, 2003;Guo et al., 2017). Considering a cohesive interface between layers, we developed the analytical expressions that give the stress distribution in a layered material as a function of the Young's modulus ratio P E , the macroscopic stress applied P , the proportion of material #1, and both Poisson's ratio 1 and 2 (equations (A7) and (A8)).
These equations allow defining the stress regime diagrams displayed in Figures 15 and 16. These diagrams show the stress ratio in the layers #1 and #2 as a function of the applied stress and the Young's modulus ratio. These figures are calculated for constant Poisson's ratio ( 1 = 0.34, 2 = 0.24) and different proportions of material #1 (0.9 in Figure 15 and 3∕7 in Figure 16). The proportion of each material strongly controls the distribution of stress across the layers. Considering the case with a high proportion of high organic content material (90%, q 1 = 0.9) (Figure 15), a wide range of parameters (E 2 ∕E 1 > 2 and P < 0.4) produces horizontal extensive stresses in the layers with low organic content ( h2 ∕ v < 0). This configuration is likely to be reached in shales of the Draupne formation (E 2 ∕E 1 ∼ 2, P ∼ 1∕3) as shown by the gray area on Figure 15 that represent Draupne range of parameters. It can also be applied to other shales formations where the variability of the Young's modulus between layers can be larger than a factor two (Sone & Zoback, 2013). Extensive horizontal stress conditions are favored in rocks with a large proportion of layers with high organic content, which can lead to mode I microfracturing in layers with low organic content, creating pathways for fluid migration. The approach developed in the present study helps discriminating when this mechanism occurs in shales.

Microfracturing and Primary Migration in Shales
Fracture propagation in shales has been extensively studied in the context of hydrofracturing for unconventional hydrocarbon production. Studies have mainly focused on the role of heterogeneities, such as pre-existing fracture planes (Zhang et al., 2019) and layer interfaces on the propagation of hydraulically driven fractures (Haddad & Sepehrnoori, 2015;Guo et al., 2017).
Fluid-driven microfracturing induced by kerogen maturation has been proposed as a possible mechanism to increase shale permeability and enhance fluid migration from source to reservoir rocks during primary migration (Pelet & Tissot, 1971;Teixeira et al., 2017;Vernik, 1994). For a homogeneous, elastic and isotropic shale, Teixeira et al. (2017) have shown that, under strong deviatoric stress, fluid-induced microfractures emanating from kerogen patches propagate toward the largest principal stress direction, allowing the fluid to migrate toward the reservoir.
Here, we specifically study the effect of structural layering in shales on the development of microfractures induced by kerogen maturation. These layers can present different elastic properties and a variety of thicknesses (Li et al., 2017). The main results are summarized in Figure 16 that displays the different stress regimes admissible in the layers with low or high organic content. The blue color shades correspond to macroscopic reverse faulting stress regimes and the red color shades correspond to normal faulting regimes. Two different domains can be defined: (1) v > h1 , where the rotation of the hydraulically driven microfractures emanating from kerogen patches is the main fracturing mechanism in layers with low organic content (Mechanism 1 in Figure 1), and (2) v < h1 , where hydraulically driven microfractures emanating from kerogen patches propagate along the bedding, within the layers with high organic content. In such case, the mechanism creating microfractures in the layer with low organic content is the compression of this layer  (2) and (3), for given Poisson's ratios ( 1 = 0.34, 2 = 0.24) and a given proportion of high organic content rock (q 1 = 0.9).
due to the overpressurization of the kerogen patches present in the organic-rich layers located below and above (Mechanism 2 in Figure 1). For given external stress conditions, microfracturing within the layers with low organic content is mainly controlled by the initial confinement I 1,#2 . If the initial vertical stress v remains constant, the increase of I 1,#2 is related to the increase of h2 ∕ v .
These results applied on Draupne shales show that, given the high deviatoric stress in this formation, the microfracturing of the low organic content layers is caused by the rotation of the hydraulically driven microfractures toward the vertical direction (Mechanism 1 in Figure 1). In other shales, other mechanisms could explain microfracturing in the low organic content layer, depending on the external stress conditions ( Figure 16).
Laboratory studies on immature shale samples have been performed and fracturing induced by kerogen maturation was visualized using either scanning electron microscopy (Allan et al., 2014;Ma et al., 2017) or time-lapse synchrotron X-ray microtomography (Kobchenko et al., 2011;Panahi et al., 2013). All these experiments were carried out without external stress applied and showed kerogen maturation-induced fracture propagation along the bedding plane, which corresponds to the preferential direction of the penny-shaped kerogen patches. Panahi et al. (2019) and Teixeira et al. (2017) presented X-ray microtomography images of a sample where a small differential stress was applied during maturation and highlighted the formation of vertical fractures, as confirmed by our simulations. Microfractures developing in low organic content layers connected with high organic content layers through channels, thus creating a three-dimensional fluid pathways for primary migration. The simulations of the present study provide additional explanation of these experimental results.  (2) and (3), for a given Poisson's ratio ( 1 = 0.34, 2 = 0.24) and a given proportion of high organic content rock (q 1 = 3∕7), corresponding to the simulations.
Given the Young's modulus E i , the Poisson's ratio i , the proportion of each kind of layer q 1 in a shale rock, and the external stress conditon P , equations (2) and (3) can be used to characterize which fracturing mechanism is most likely to occur within the low organic content layers.

Conclusion
The influence of structural layering in shales on the microfracturing process during kerogen maturation is investigated for different stress conditions using hydro-mechanically coupled numerical simulations. The nucleation and growth of microfractures emanating from kerogen patches increases the permeability and favors primary migration from the source rock to the reservoir. Figures 15 and 16 summarize the following results: 1. Hydraulically driven microfractures form due to local fluid overpressure induced by kerogen maturation. They emanate from the kerogen patches and may propagate and rotate in the direction of the local largest principal stress. Under normal faulting stress conditions, this direction is vertical and flow pathways connected in three dimensions may form (Mechanism 1, black microfractures in Figure 1). This process is directly controlled by the state of stress in the layers with high organic content. 2. The local stress variations due to the contrast of elastic parameters between kerogen-rich and kerogen-poor layers can induce mode I fracture in shales (Bourne, 2003;Boersma et al., 2020). Figure 15 shows that this mechanism can occur in the Draupne shale formation. This mechanism is favored when the proportion of low organic content layer is low. Equations (2) and 3 can be used to decipher if this mechanism is likely to occur. 3. If the local horizontal stress is larger than the vertical stress in the layers with high organic content, hydraulically driven microfractures caused by kerogen maturation propagate within these layers. Then, damage in the layers with low organic content is caused by the compression of those layers induced by the pressure buildup of the surrounding high organic content layers (Mechanism 2, white microfractures in Figure 1). The main parameter controlling this fracturing process is the initial confinement of the layers with low organic content (I 2,#2 ), which results from the elastic contrast between the layers. 4. Stress redistribution due to the layering controls microfracturing of the low organic content layers. Under high deviatoric stress, such as in the Draupne formation, rotation of the microfractures emanating from the kerogen patches is the main microfracturing mechanism (Mechanism 1 in Figure 1). Under reverse faulting stress conditions, microfracturing by compaction of the low organic content layers can be the dominant mechanism (Mechanism 2 in Figure 1).
The present study proposes a tool to characterize which microfracturing mechanism likely occurs in a layered shale rock for a given stress configuration based on equations (2) and (3). Further studies could be carried out in the future to better constrain these mechanisms with regard to the complexity of sedimentary basins by taking into account true triaxial stress conditions ( v ≠ h1 ≠ h2 ), a more complete description of the anisotropy of mechanical properties, and geochemical models of kerogen maturation.

Appendix A: State of Stress in a Layered Linear Elastic Material
A layered rock made of two linear elastic isotropic materials is considered. Both materials are characterized by a Young's modulus E i , a Poisson's ratio i , and a volume fraction q i (with q 1 + q 2 = 1), as shown in Figure  A1. The applied loading is assumed to be transverse isotropic, with ( v , h ), respectively, perpendicular and parallel to the layers ( Figure A1). In this configuration, the vertical stress, along the z-axis, is the same in both layers and equal to v .
Applying Hooke's law on an elastic isotropic material (equation (1)), we can project the strain tensor for both materials on the ⃗ x axis. The strain along the ⃗ x axis is geometrically constrained to be the same for both materials and is noted t .
and h1 and h2 are related to v through: q 1 h1 + (1 − q 1 ) h2 = h (A3) Figure A1. Schematic representation of a layered rock with imposed elastic parameters and loading conditions.

10.1029/2020JB019444
Solving the set of equation (A4) and introducing the parameter P = h v , it is possible to calculate analytically hi v from six known parameters: hi v = g(E i , i , P , q 1 ) (equations (A5) and A6)).
Equations (A5) and (A6) can be written by using the Young's modulus ratio P E = E 1 ∕E 2 . Five independent parameters are needed to calculate the ratio hi ∕ v , the Young's modulus ratio, the applied stress ratio P , the proportion q 1 of material #1, and the Poisson's ratio of both layers 1 and 2 . , macroscopic first stress invariant I 1,#1 local first stress invariant in layers with high organic content I 1,#2 local first stress invariant in layers with low organic content V f injected fluid volume