Simulation of Fluid Flow in Induced Fractures in Shale by the Lattice Boltzmann Method

Abstract

With increasing interest in unconventional resources, understanding the flow in fractures, the gathering system for fluid production in these reservoirs, becomes an essential building block for developing effective stimulation treatment designs. Accurate determination of stress-dependent permeability of fractures requires time-intensive physical experiments on fractured core samples. Unlike previous attempts to estimate permeability through experiments, we utilize 3D Lattice Boltzmann Method simulations for increased understanding of how rock properties and generated fracture geometries influence the flow. Here, both real induced shale rock fractures and synthetic fractures are studied. Digital representations are characterized for descriptive topological parameters, then duplicated, with the upper plane translated to yield an aperture and variable degree of throw. We present several results for steady LBM flow in characterized, unpropped fractures, demonstrating our methodology. Results with aperture variation in these complex, rough-walled geometries are described with a modification to the theoretical cubic law relation for flow in a smooth slit. Moreover, a series of simulations mimicking simple variation in proppant concentration, both in full and partial monolayers, are run to better understand their effects on the permeability of propped fractured systems.

1 Background

The structure of fractures and their complexity will generally depend on two main conditions: (1) the type of material undergoing the breakage, determined by its physical and chemical properties and (2) total mechanical stress acting on the material. The general definition of a fracture for the purpose of this work is any discontinuity in a rock volume created during rupturing of a rock mass, generating surface with annihilated cohesion.

In its most simple representation, we can imagine a fracture as flow in a smooth slit. Laminar flow between smooth parallel plates can be posed and solved analytically [1], yielding a cubic law relationship between permeability (k) or transmissivity (T) and geometrical parameters: aperture (h), cross sectional area (A), and width (W).

$$ T = kA = \frac{{Wh^{3} }}{12} $$

(1)

Flow in real fractures, however, is impacted by wall roughness, aperture, and shear displacement. According to Hakami [2], there are three main aspects to fluid flow prediction in single fractures (see Fig. 1): fluid properties, fracture void geometry, and imposed boundary conditions. We focus on the elements in boldface.

Fig. 1.

Figure modified from [2].

Main factors altering the flow in a single fracture. Context of this work is in bold italic.

Although once considered of no commercial hydrocarbon potential, with recent advancement in completion and stimulation practices in long horizontal wells, oil and gas production from shale formations is made economically feasible and is an increasingly important contributor in the fossil fuels portfolio with global potential [3]. We focus on the fracture void geometry of fractured shale and the subtopics of geometrical parameterization of the system and physical changes produced by the fracturing process, namely, normal and shear displacement.

1.1 Surface Characterization

Surface roughness can be regarded as any irregularity or deviation of the surface structure from the mean smooth plane. The larger these deviations, the rougher the surface, and in case they are relatively small, a surface is considered smooth. Fracture roughness is controlled by stress conditions which affect the crack propagation pattern, lithology of the rock matrix, including all types of heterogeneities, and finally, secondary physical and chemical processes, such as weathering, erosion, and mineral precipitation.

Standards from the field of tribology, studying effect of surface roughness on lubrication, wear and friction [4], and some amplitude parameters proposed by Stout and Blount [5] were considered. The most straight-forward parameter to describe roughness, representing the mean of absolute profile height deviations from the mean plane, is average roughness (S_a).

$$ S_{a} = \frac{1}{{N_{x} N_{y} }}\sum\nolimits_{j = 1}^{{N_{y} }} {\sum\nolimits_{i = 1}^{{N_{x} }} {\left| {z\left( {x_{i} ,y_{j} } \right) - m} \right|} } $$

(2)

where z is the vertical height at any point, N_x and N_y are total number of data points in x and y direction respectively, and m is the mean of all data points.

$$ m = \frac{1}{{N_{x} N_{y} }}\sum\nolimits_{j = 1}^{{N_{y} }} {\sum\nolimits_{i = 1}^{{N_{x} }} {z\left( {x_{i} ,y_{j} } \right).} } $$

(3)

However, the most widely used parameter, known in statistics as standard deviation, is root mean square (RMS) roughness (S_q) given by

$$ S_{q} = \sqrt {\frac{1}{{N_{x} N_{y} }}\sum\nolimits_{j = 1}^{{N_{y} }} {\sum\nolimits_{i = 1}^{{N_{x} }} {\left( {z\left( {x_{i} ,y_{j} } \right) - m} \right)^{2} } } } $$

(4)

We also have the difference between extrema, S_t, skewness, $S_{sk}$, and kurtosis, $S_{ku} .$

$$ S_{t} = \hbox{max} \left( {z\left( {x_{i} ,y_{j} } \right)} \right) - \hbox{min} \left( {z\left( {x_{i} ,y_{j} } \right)} \right) $$

(5)

$$ S_{sk} = \frac{1}{{N_{x} N_{y} S_{q}^{3} }}\sum\nolimits_{j = 1}^{{N_{y} }} {\sum\nolimits_{i = 1}^{{N_{x} }} {\left( {z\left( {x_{i} ,y_{j} } \right) - m} \right)^{3} } } $$

(6)

$$ S_{ku} = \frac{1}{{N_{x} N_{y} S_{q}^{4} }}\sum\nolimits_{j = 1}^{{N_{y} }} {\sum\nolimits_{i = 1}^{{N_{x} }} {\left( {z\left( {x_{i} ,y_{j} } \right) - m} \right)^{4} } } $$

(7)

Other statistical characterizations of note are the joint roughness coefficient (JRC) [6], semi-variogram [7], and fractal dimension [8].

1.2 Fracture Characterization

While there are many reasonable ways to characterize roughness of a single surface, flow in fractures is between two similar surfaces. Many of the above concepts can be generalized to properties of the plane pair or the space they create, such as a semi-variogram on point aperture. The surfaces are similar due to the fracture generation process from an intact material. However, the fracturing process generates debris, giving dislocations in one or both surfaces, and any shear displacement creates aperture distributions with potential impingement “pillars” at points of contact and potential elimination of spatial correlation regarding flow, despite similarities only a short distance away. Since shale has significant clay and organic content, the resulting fracture is also sensitive to further imposed stress with plastic flow or creep, yielding apertures changing with time.

When considering flow in the fracture, we also have the notion of tortuosity (τ = L_a/L), as the arclengths for streamlines, L_a, are longer than the sample length, L. Yet another measure of tortuosity, T_s, was defined by Belem et al. [9] as a normal component of the average roughness ratio, R_s, which could also be extended to fractures as the average between top and bottom surfaces. The average roughness ratio compares the actual surface area to that of a nominal plane spanning the same region.

1.3 Flow Characterization

Early attempts were experimental in nature [10,11,12] and explored extensions of flow in rough conduits with a relative roughness parameter, ε/2h, where ε is absolute roughness height and 2h is equal to hydraulic diameter. Witherspoon et al. [13] incorporated a friction factor to model experiments on granite, basalt, and marble fractures. Other authors introduced a relative roughness ratio [14,15,16], the JRC coefficient [17], and a contact area fraction [18], c, to account for periodic collapse of the flow area by upper and lower surface contact. Many of these cubic law corrective models strongly depend upon the nature of the average aperture used.

2 Procedure

2.1 Experimental

We consider real rough-walled profiles with both uniform and variable aperture field distributions, including taking into account asperities at points of contact between the two surfaces. All of the models described herein are based on 3D surface profiles acquired from longitudinal fractures created by Brazilian tests performed on four 1 in. diameter, 1.5 in. long shale core samples that did not contain visible macroscopic fractures (see Fig. 2).

Fig. 2.

Fracture creation using the Brazilian test procedure.

The Brazilian test is commonly used for material tensile strength determination and can essentially be described as uniaxial normal stress compression process leading to failure.

Resulting fracture faces are digitized for surface topography using a commercial optical profilometer manufactured by Nanovea (see Fig. 3). The main characteristics of any profilometer are its maximum measurement range and resolution in x and y directions, determined by the stepping of the sample stage holder, and resolution, measurement range and accuracy in the z direction, dictated by the particular optical pen installed. Our pen had a measurement range of 27 mm with a vertical resolution of 600 nm and vertical accuracy of 3000 nm. We used steps of 100 µm in both x and y directions in the surface characterization of all samples. The data required preprocessing to convert the cloud point data to an STL mesh that was importable to Exa’s PowerDelta tool to ensure a dense compatible boundary suite that would serve as no-slip boundary conditions in PowerFlow, Exa’s Lattice Boltzmann Method flow simulator.

Fig. 3.

Surface profilometer used for data acquisition.

Figure 4 indicates the two types of simulations that dictate the level of geometrical characterization needed to correlate with flow behavior. We could entertain exclusively normal translation of the upper surface to achieve a nearly constant aperture, or we could allow normal and shear displacement to yield two unmatched surfaces and variable aperture.

Fig. 4.

Flow geometries without and with displacement in the direction of fracturing.

Due to the debris created in the fracturing process, however, we did not have ideal mating surfaces. Rather than work with both surfaces with dislocations and a possible debris field, we instead duplicated and translated a digital representation of the lower surface to create an upper surface and flow volume. This is illustrated in Fig. 5 where we had to construct the two surfaces, define the flow field, and encapsulate the open fracture into a regular solid geometry.

Fig. 5.

Dual surfaces to fracture to 3D volume transformation for simulation input.

In field stimulation operations, the fractures are created by exerting enough pressure transmitted by a fluid to crack the rock. Release of the pressure would collapse the fracture, which could, in some cases, heal. To maintain an open conduit, suspended sand or other proppant material is included in the pressurizing fluid. The proppant, often coated to produce a bond with the formation to avoid return of solids, maintains an aperture consistent with the sand grain diameter. We simulate the permeability of a proppant-ladden fracture as spherical elements placed between the fracture faces. Rather than simulate the delivery process, we investigate the impact of proppant density as deviations from monolayer coverage. That is, we start with a closed packing arrangement and create lower coverage with random removal of particles.

2.2 Computational

Method.

We employ a commercial, three dimensional, Lattice Boltzmann simulator (PowerFlow) provided by Exa Corporation. Application of Lattice Boltzmann Method (LBM) has several advantages compared to other numerical simulation schemes [19]:

• 1. Mesoscopic level of operations – simplified microscopic description allows more natural description of small scale effects and detailed modelling of highly complex geometric boundaries.

• 2. Simple and automatic volume discretization – simulation volumes can easily be meshed into a lattice of cubic voxels without any need for solid boundary adaptations.

• 3. Inherent parallelism – simulations can be run on multiple processors due to the localized nature of the performed operations.

Even though LBM is highly parallel, computational demand required for simulations to converge is quite high – larger domains of simulation can require somewhere from a couple of days up to a week of run time on 100 cores to obtain reasonable results.

Unlike molecular dynamics, instead of analyzing fluid as a collection of individual particles, LBM treats fluid volume as “collection of particles represented by a particle velocity distribution function at each grid point” [20]. Another critical part of the LBM development is introduction of a simplified collision operator introduced by Bhatnagar, Gross and Krook [21], which considers a stencil with 19 possible nodes for discrete particle velocities (see Fig. 6).

$$ u_{lat} = \frac{{\nu_{lat} Re}}{Resolution} $$

(8)

$$ u_{max} = \frac{{h^{2} }}{8\mu }\frac{dP}{dL} = \frac{3}{2}u_{lat} $$

(9)

$$ g = \frac{dP}{dL}/\rho_{lat} $$

(10)

where $u_{lat}$ in all of the equations represents average lattice velocity, $\nu_{lat}$ is kinematic lattice viscosity, $u_{max}$ is maximum lattice velocity, $\rho_{lat}$ is lattice density and g is gravity or acceleration applied to induce fluid flow. Notice that the main two parameters used to control other simulation parameters are Reynold’s Number, Re, and resolution, λ.

Fig. 6.

D3Q19 lattice arrangements for 3D problems [21].

$$ k_{{lat_{sim} }} = \frac{{\nu_{lat} u_{{lat_{sim} }} }}{g} $$

(11)

$$ k_{sim} = k_{{lat_{sim} }} \left( {\frac{{l_{char} }}{\lambda }} \right)^{2} $$

(12)

After reaching a steady state, we extract average simulated lattice velocity $u_{{lat_{sim} }}$, which we further use to calculate simulated permeability. Conversion from permeability in squared lattice length units to square meters is done using the voxel size which essentially is determined by the ratio of characteristic length $l_{char}$ and λ. Characteristic length in all of our simulations is chosen to be equal to the aperture height separating surfaces of a fracture.

An important task before the start of simulations with real rock surfaces is to perform verification exercises. This is achieved by running a benchmark case with parallel plate geometry and comparing simulation results with analytical solution given by Eq. 13.

$$ u\left( z \right) = 4u_{max} \frac{z}{h}\left( {1 - \frac{z}{h}} \right) $$

(13)

where z is elevation above the bottom plate surface in the range [0, h]. Comparison of the velocity profiles and computed permeability with varying aperture shows a good match (see Fig. 7) between simulation and theory.

Fig. 7.

Comparison of analytic and LBM simulation velocity profiles and permeability values for steady state flow between parallel plates.

Based on these results and on the measure of error between analytic and LBM permeability always lower than 0.7% for all the steady state Poiseuille flow cases in Fig. 7, we assume the correctness of the simulator.

Validation.

To determine proper values for Re and λ, we once again run a series of test simulations. First runs were made using parallel plate geometry with h = 0.8 mm. Error percent was calculated using the analytic cubic law solution.

Resolution.

Based on Fig. 8, we can easily conclude that the most appropriate Resolution to run the simulations is ten lattice units. If we go above this number, we increase accuracy, but the improvement will not be significant and will incur a higher cost in both discretization and simulation time.

Fig. 8.

Permeability change with Resolution variation for parallel plate geometry

Reynolds Number.

The next critical parameter affecting both accuracy and speed control is Re. Flow in smooth parallel plates indicated an abrupt change in flow behavior for Re ≥ 20, although errors compared to parallel plate analytical solution were less than 1%. Flow between two complementary rough surfaces separated by a 0.2 mm vertical aperture gave the results presented in Fig. 9. In the absence of reference solution, we assumed the result from the simulation with the smallest Re number to be the most correct one and a basis for computed error. Based on the rough fracture simulations, we concluded that the maximum Re number value we should utilize to run our simulations without sacrificing accuracy is ten.

Fig. 9.

Permeability change with Re number variation for rough mated surfaces (h = 0.2 mm).

3 Results

Summary statistics on four fractured shale samples are given in Table 1. LBM simulation results are provided in Table 2.

Table 1. Summary statistics of four core samples

Table 2. LBM simulation results for mated surfaces

3.1 Aperture

Based on a series of simulations with apertures varying form 0.02 mm up to 0.4 mm, we propose the correlation

$$ k = \frac{{h^{2} }}{{12\tau_{x} \tau_{y} T_{s}^{2} }} $$

(14)

which is compared with the simulation and smooth slit cubic law results in Fig. 10. Since these simulations only translate the duplicate fracture face normal to the direction of flow, the aperture is everywhere constant. The correlation includes measures of tortuosity in both x and y directions, since flow is actually 3D in nature, and a third measure of tortuosity, T_s, based upon area roughness ratio.

Fig. 10.

Comparison of simulation results with cubic law and proposed equation for surface A1.

3.2 Shear

When considering lateral or shear translation, aperture is no longer constant, as illustrated in Fig. 11, and Eq. 14 can no longer adequately capture the added complexity.

Fig. 11.

Typical aperture distribution created with vertical displacement plus lateral shear.

We propose modifications to arrive at a new correlation shown in Eq. 8.

$$ k = \frac{{h_{G}^{2} }}{{12\tau_{x} \tau_{y} T_{s}^{2} \left[ {1 + \frac{{S_{q} }}{{2h_{G} }}} \right]}} \left( {\frac{1 - c}{1 + c}} \right) $$

(15)

where the RMS roughness, S_q, and contact area fraction, c, are also introduced. A variety of mean apertures were tested: arithmetic mean aperture, geometric average aperture, arithmetic mean excluding zero values (NNZ), and the mechanical aperture (distance between mean planes) with the geometric mean aperture, h_G, providing the closest match in Eq. 15 to the simulated permeability, as shown in Fig. 12.

Fig. 12.

Results generated using Eq. 4 for different average apertures vs simulation results (surface A1). Legend: Geom – $h_{G}$, Average – $h_{A}$, Average_NNZ – $h_{\varvec{A}}$ excluding zero apertures, Mech - $h_{m}$.

3.3 Proppant Density

When using a relatively large dispersed solid to bridge surfaces and maintain an aperture, it is believed that the impact of the surface asperities and roughness play a much less dominant role. The roughness and shape distribution of the proppant could possibly become the focal point affecting resulting flow behavior, but these elements are controllable. The preliminary work in this area, therefore, uses smooth parallel plates with variable concentration of proppant to examine first order effects. Starting with closed packing and monolayer filling, lower concentrations (f < 1) were obtained by random removal of spheres. The observed relation for concentration dependent permeability is displayed in Fig. 13 and captured by the curve fit

$$ k = 7.78 \times 10^{4} \,e^{ - 4.14f} $$

(16)

with a relatively high correlation coefficient (R² = 0.905). The scatter seen in this figure at low values of concentration is due to multiple placement patterns at each of these concentrations. While field efforts target placement of increasing proppant volume [22], LBM simulations clearly demonstrate the value in only partial coverage, provided proppant is delivered throughout the rock failure zone.

Fig. 13.

Permeability and proppant concentration relationship.

4 Conclusions

A procedure was assembled for systematic characterization of induced fractures, with regard to impact on fluid flow, and their re-assembly to create digital fractures of uniform aperture and those with aperture distributions created through shear displacement of one face. In both cases, the simple flow in a smooth slit model required empirical modification using phenomenological parameters related to either frictional drag or pathlength extension to successfully represent the pressure drop relationship to flow rate. Elemental work with propped fractures supplements flow characterization with size and concentration of artificially indicated the value in partial fracture filling in sustainable fracture conductivity (Fig. 14).

Fig. 14.

Velocity streamlines in full monolayer (top) and partial monolayer (bottom). Uniform flow through all channels in full monolayer while in partial monolayer flow mostly occurs through the biggest channel.

5 Future Work

It is highly desirable to extend this work to multiphase flow due to the known bimodal wettability in such rocks with organic and inorganic porosity. The rock matrix was treated as impermeable, though real systems actually feed the fracture. Additionally, the rock matrix can be treated as a system under stress and undergoing plastic flow, resulting in partial embedment and loss of permeability. The impact of fracture roughness in proppant studies was ignored and should be investigated, as well as characterization of shape, distribution, and surface roughness of proppant. A test matrix that includes all the major commercial shale plays would be advantageous, as the surface properties examined should be functions of brittleness, related to mineralogy and organic content.