Solvable model of gas production decline from hydrofractured networks

We address questions that arose from studying gas and oil production from hydrofractured wells. Does past production predict the future? This depends on deducing from production as much as possible about the plausible geometries of the fracture network. We address the problem through a solvable model and use kinetic Monte Carlo and Green's function techniques to solve it. We have three main findings. First, at sufficiently long times, the production from all compact fracture networks is described by a universal function with two scaling parameters, one of which is the diffusivity of unbroken rock $\alpha$ and the second of which is a parameter ${\mathcal{V}}_{\mathrm{ext}}$ with units of volume. Second, for fracture networks where the power-law distribution of fracture spacings falls below a critical value (and this appears to be the case in practice), early-time production always scales as one over the square root of time. Third, the diffusivity $\alpha$ that sets the scale for late-time production is inherently difficult to estimate from production data, but the methods here provide the best hope of obtaining it and thus can determine the physics that will govern the decline of unconventional gas and oil wells.

Figure 1

The geometry of the parallel plate model [6] consists of a horizontal well with a uniform array of planar hydrofractures. The cuboid volume that contains the well and the hydrofractures is known as the stimulated reservoir volume. Gas transport in this model involves diffusion of gas pressure in the interior of the stimulated volume and normal to the hydrofractures. Gas transport from the exterior is neglected.

Figure 2

Geometry of our Monte Carlo model, with a stimulated reservoir volume (fractures in blue) embedded in a fracture-free matrix (white). The precise structure of the fracture network comes from a geologically constrained percolation model described in detail in [17, 18]. In the model, fracture length results from a power-law probability distribution characterized by the exponent $e$. The density of the fracture network is controlled by the filling fraction $p$, defined as the ratio of the broken lattice bonds to the total number of bonds.

Figure 3

Validation of our Brownian motion code. To construct the scaled rate vs time plot, we fit the parallel plate solution (29) to the probability density function for arrival time, which we get from the simulation. The fitting procedure yields a timescale $\tau$ and a mass scale $\mathcal{M}$, which we then use to scale time and rate: scaled time $\tilde{t}$ is defined by $t/\tau$; scaled rate, by $\dot{Q}/(\mathcal{M}/\tau ).$ (a) Validation geometry, a fractured (blue) $5\times 5$ lattice. (b) Simulation (jagged6 line) matches the theory from the parallel plate model (smooth line).

Figure 4

Production decline for the fracture network in Fig. 2 (jagged line). The rate initially declines according to the solution for the parallel plate geometry, Eq. (29) (smooth curving line), but then transitions to a power-law decline (straight line).

Figure 5

Geometry and decline rate for structure with $e=2$ and $p\approx 40\%$. (a) Fracture network with $e=2$ and $p\approx 40\%$. (b) Production rate (jagged line) initially declines according to the solution for the parallel plate geometry (smooth curve), Eq. (29), but then transitions to a power-law decline (straight line).

Figure 6

Illustration of the interior of a fractured network. The absorbers are the black squares, and these are indicated by the condition ${\theta }_j=0$. The space $L$ is the space of all sites interior to the absorbers but not absorbers themselves: light blue (or gray), indexed by $l$.

Figure 7

Illustration of the exterior of a fractured network. Space $M$ is the set of all absorbers in contact with the exterior: black squares, indexed by $m$. Space $M^{\prime }$ is the set of all nearest neighbors of absorbers that are not themselves absorbers and in contact with the exterior: orange (or gray), indexed by $m^{\prime }$. Space $N$ is the union of $M$ and $M^{\prime }$: black and orange (or gray) squares, indexed by $\nu$. Space $O$ is the set of absorbers in contact with the exterior and all exterior sites: black squares plus everything outside them.

Figure 8

Spectrum and production for a single absorber. (a) Continuous spectrum on the negative real axis for single absorber. (b) Gas production, single absorber.

Figure 9

(a) Square with interior of size $64\times 64$. (b) Interior production from the square, comparing exact computation ${\dot{Q}}_{\mathrm{interior}}$, Eq. (15), with approximate analytical expression ${\mathcal{Q}}_{\mathrm{int}}$, Eq. (29), scaled as described in the text.

Figure 10

Cumulative production for flow to network shown in Fig. 9. From top to bottom, total, exterior, and interior production.

Figure 11

Network of absorbers used as example.

Figure 12

Continuous spectrum from network of Fig. 11. (a) Continuous spectrum along negative real axis for network in Fig. 11. (b) Real part of continuous spectrum along contour in complex plane in Fig. 23, as function of real part of contour variable.

Figure 13

(a) Gas production rate for network in Fig. 11. (b) Cumulative gas production from interior of network in Fig. 11.

Figure 14

Lattice realizations of the fracture network of Fig. 5 at three separate resolutions. (a) Fracture network of Fig. 5 on lattice whose spacing equals minimum distance between fractures. (b) Fracture network of Fig. 5 on lattice whose spacing equals half the minimum distance between fractures (c) Fracture network of Fig. 5 on lattice whose spacing equals one sixth the minimum distance between fractures.

Figure 15

Comparison of kinetic Monte Carlo and lattice Green's function result for fracture network of Fig. 5. (a) Convergence of Green's function solution as lattice of Fig. 5 is scaled by factors $s$ (from lowest curve to highest curve) of 1, 2, 4, and 6, and thus approaching a continuum limit where the early-time behavior is a $t^{-1/2}$ power law. The scaled time is time divided by the scaling factor squared, leading the curves to lie on top of one another. (b) Comparison of Green's function solution scaled by factor of 6 with kinetic Monte Carlo solution in Fig. 5.

Figure 16

Late-time and very-late-time production functions ${\dot{Q}}_{\mathrm{lt}}\left(t\right)$ and ${\dot{Q}}_{\mathrm{vlt}}$ for $\sigma ={\sigma }_0.$ The late-time expression produces a value of $4\pi$ for $t\ll 1$, but this is not meaningful; only $t\gg 1$ provides information about production rates.

Figure 17

Grid geometry with $N=12$.

Figure 18

Plots of cumulative production when the probability of interference time $\tau$ goes as ${\tau }^{\beta }$ up to a maximum interference time ${\tau }_m$. Curves are ordered with $\beta =0.9$ on top and $\beta =-0.7$ on bottom. (a) As $\beta$ drops below 0.5 the early-time asymptotic behavior of ${\mathcal{Q}}_{\mathrm{int}}(t,\beta )$ stops depending upon $\beta$. (b) However, the detailed shape of the cumulative production curve in the vicinity of $t\approx {\tau }_m$ does depend upon $\beta ,$ and this does affect attempts to provide good phenomenological fits to exact solutions of model systems.

Figure 19

Multiscale rectangular network used as test case.

Figure 20

Fits of the model in Eq. (58) to the exact production rate and cumulative production for the structure in Fig. 19. (a) Production rate, comparing the exact result ${\dot{Q}}_{\mathrm{Total}}$ to the parameterized ${\dot{\mathcal{Q}}}_{\mathrm{grid}\beta }$. (b) Cumulative production, comparing the exact result $Q_{\mathrm{Total}}$ to the parameterized ${\mathcal{Q}}_{\mathrm{grid}\beta }$.

Figure 21

Slope of log-log production plot, comparing models with parameterized fits. (a) Slope of log-log cumulative production including fractal exponent $\beta$ from Eq. (58). (b) Slope of log-log cumulative production without fractal exponent $\beta$ in Eq. (49).

Figure 22

Comparison of fits obtained to production data as maximum available time increases. The bands around each line show one-$\sigma$ uncertainties obtained from the statistical routine conf_interval in the package lmfit [22]. The uncertainties are almost always underestimates. (a) Estimates of model parameters, including long-time behavior, making use of fractal exponent $\beta$. (b) Estimates of model parameters, including long-time behavior, without making use of fractal exponent $\beta$.

Figure 23

Integration contour used for Eq. (D5). The values used in computation were ${\lambda }_1=-8+{10}^{-5}, \alpha =-0.1,{\lambda }_2=-{10}^{-5}.$