Prediction of decline in shale gas well production using stable carbon isotope technique

Abstract

Prediction of shale gas production is a challenging task because of the complex fracture-pore networks and gas flow mechanisms in shale reservoirs. Empirical methods, which are used in the industry to forecast the future production of shale gas, have not been assessed sufficiently to warrant high confidence in their results. Methane carbon isotopic signals have been used for producing gas wells, and are controlled by physical properties and physics-controlling production; they serve as a unique indicator of the gas production status. Here, a workable process, which is combined with a gas isotope interpretation tool (also known as a numerical simulator), has been implemented in Longmaxi shale gas wells to predict the production decline curves. The numerical simulator, which takes into account a convection-diffusion-adsorption model for the matrix and a convection model for fractures in ¹³CH₄ and ¹²CH₄ isotopologues, was used to stabilize the carbon isotope variation in the produced gas to elucidate gas recovery. Combined with the production rates of the four developing wells, the total reserves ranged from 1.72 × 10⁸ to 2.02 × 10⁸ m³, which were used to constrain the trend of two-segment production decline curves that exhibited a transition from a hyperbolic equation to an exponential one within 0.82–0.89 year. Two-segment production decline curves were used to forecast future production and estimate ultimate recovery.

Graphical abstract

Keywords

shale gas / production decline / Longmaxi formation / carbon isotope

1 Introduction

In the past several decades, great achievements have been made in the geological characterization of shale gas around the world. However, prediction of shale gas production is still a challenging task, partly because the existing gas transport models have not fully considered the complex mechanisms in the nanopore systems (Ilk et al., 2010; Lee and Sidle, 2010; Wang et al., 2015). Consequently, analytical, semi-analytical, and numerical solutions based on these gas transport models cannot accurately predict gas production (Moinfar et al., 2016; Mahmoud et al., 2018).

Regardless of the fluid mechanics in shale, the production rate is considered as a direct indicator of good performance and is crucial for forecasting future well production. On this basis, a decline curve analysis (DCA) that uses regression for historical production data was proposed to predict well production. The only input data required for this method are production data and basic reservoir information (Tan et al., 2018). The classical DCA model is based on the Arps equation (Arps, 1945), which is a hyperbolic function with three parameters (initial production rate, initial decline rate, and decline exponent) that has been used to simulate the decline in production rate. The Arps model was originally designed for boundary-dominated flow, which occurs in formations with medium and high permeability. However, the production performance of horizontal fractured wells in low-permeability reservoirs (e.g., shale gas or tight gas) is typically characterized by multiple flow regimes. Initially, a well produces a transient flow, and the production rate usually shows a rapid drop. At a later production time, the production rate drop slows down, indicating boundary-dominated flow. If the transient solution is represented by the hyperbolic Arps equation with a decline exponent greater than 1 (Currie et al., 2010; Lee and Sidle, 2010; Okouma et al., 2012), it would lead to an unbounded and unrealistic prediction of the estimated ultimate recovery (EUR).

To resolve this problem, a two-segment Arps equation was employed for different production time intervals (Seshadri and Mattar, 2010; SPEE, 2016). The hyperbolic Arps equation used in the early stage is replaced by the exponential Arps equation in the later stage when the decline rate reaches a predetermined value. Such a predetermined decline rate is empirical and has no physical basis (Meyet et al., 2013). Furthermore, several alternative models have been proposed for analyzing the decline curves of shale gas wells, such as MHD (Robertson, 1988), PLE (Ilk et al., 2008), SEDM (Valko, 2009), Duong (Duong, 2011), EEDCA (Zhang et al., 2016), and FDC (Zuo et al., 2016). Several comparative studies based on these DCA models show great differences when forecasting the production of the gas derived from Barnett and Eagle Ford shale formations (Kanfar and Wattenbarger, 2012; Mishra, 2012; Joshi and Lee, 2013). In general, there are approximately ±10% errors in the 30-year EUR forecasts based on various mathematical theories. Practically, it is difficult to determine the model that provides the best prediction because each result seems plausible under specific hypothetical conditions.

More recently, Gao et al. (2017) proposed a new method based on methane carbon isotope variation to predict the production decline of shale gas. The key step of this method is determining the critical point during shale gas production when piecewise functions are applied to forecast the production decline. Because of the occurrence of organic matter and nanoscale pore systems within shale matrix, ¹³CH₄ isotopologue exhibits a stronger bond with organic matter and nanoscale pores than that demonstrated by ¹²CH₄-isotopologue, resulting in a faster flow rate of the ¹²CH₄ isotopologue as compared to that of ¹³CH₄-isotopologue. This results in the formation of lighter and heavier methane carbon isotopes in the produced gas during the early and late production periods, respectively (Gao et al., 2017; Zhang et al., 2018; Li et al., 2020). A minima is observed in the profile of the methane carbon isotope of the produced shale gas. Gas recovery was determined by monitoring the presence of a minima, which was used to calculate the total reserve with regard to the cumulative production (Gao et al., 2017). In addition, the critical point was derived from the point where the integral of the decline curve becomes equal to the quotient of cumulative production and recovery, while the piecewise decline curve was determined to forecast the EUR.

In this study, we introduced a comprehensive isotope fractionation model for predicting the production of four horizontal shale gas wells located in the Luzhou deep shale gas field in China. Through the carbon isotope variations of the produced gas measured over the last nine months at the well site, total reserves and two-piece decline curves were determined from the simulation based on the isotope fractionation model. The tentative application of carbon isotopes in shale gas production forecasts indicates that isotope evolution can help predict future production, especially when the fractionation mechanisms are clear. This novel method may be a useful tool in the future.

2 Geological settings

The study area of the Luzhou Shale Gas Field, which has been the most profitable deep shale gas field in China, is located in the south of the Sichuan Basin (Fig. 1(a)). Black organic-rich shale with visible graptolite fossils in the Upper Ordovician Wufeng Formation and Lower Silurian Longmaxi Formation was deposited under euxinic to anoxic conditions in a deep-shelf depositional environment (Yang et al., 2017); it was proved that it can be considered as the main production zone. High-maturity Longmaxi shale exhibits low porosity and low permeability; however, the gas content can be higher than 7 m³/t rock (Ma et al., 2020). In particular, shale formations from the Wufeng Formation. to the first member of the Longmaxi Formation. (Long-1) exhibit a total gas content of up to 15% (Fig. 1(b)) and are stably distributed laterally.

In this study, four horizontal wells (A, B, C, and D) from the Luzhou Deep Shale Gas Field were drilled at the bottom of the Longmaxi Formation. The Long-1 profiles of wells with different dates of production were investigated. The production in Well A was initiated on November 15, 2018. Wells B, C, and D were from the same platform, and their production was initiated on May 28, June 6, and June 20, 2019, respectively. The current production rates for wells A, B, C, and D are 5.6 × 10⁴, 7.8 × 10⁴, 9.2 × 10⁴, and 8.3 × 10⁴ m³/day, respectively.

3 Data and method

3.1 Real-time data collection for methane carbon isotope ratios

For the purpose of real-time measurement of methane carbon isotope ratios (δ¹³C₁) of the produced gas samples, two sets of gas chromatographic-isotope ratio infrared spectrometer (GC-IRIS; imported from Power, Environmental and Energy Research Institute or self-developed with similar performance) were deployed in the well field from September 2018 to May 2019. During this period, gas samples from wells A, B, C, and D were collected using IsoTubes almost every week, and methane carbon isotope ratios were measured via the imported infrared spectrometers for wells A, B, and C. Self-developed infrared spectrometer was applied to well D.

GC-IRIS integrated gas chromatography and infrared spectrometry were used for performing carbon isotope measurements (Fig. 2(a)). The laser measurement system is advantageous in the well field because mass spectrometry would otherwise require a significantly high vacuum level, clean power supplies, and stable environment temperature for detection.

GC was used to separate the natural gas components into methane, ethane, and propane (Fig. 2(b)). A combustion oven transformed all hydrocarbon gases into carbon dioxide. Direct absorption of the CO₂ obeys the Beer-Lambert rule, and the absorption spectrum is precisely measured to determine the isotope ratio (Fig. 2(c)). The analytical error was approximately 0.4‰, and the typical measurement time for each natural gas sample was approximately 3 min.

3.2 Production data collection in study wells

The production rate data (× 10⁴ m³/day) of well A were recorded from the date when it was first put into production (October 31, 2019). The production data of wells B, C, and D were recorded from their own dates when they were put into production (October 24, 2019). Raw daily production data were normalized to each day (24 h). Some days with no data due to workover of the well or pressure build-up were abandoned in this study. The cumulative production for any interval can be calculated from the production rate data. The hyperbolic Arps equation was used to fit some of the processed production data. The parameters for the hyperbolic decline function were determined using the least-squares regression method.

3.3 Simulation and prediction of isotopic signals

The dual-porosity model consisted of fractures and matrices. The pressure gradient toward the fracture and well causes the methane in the matrix to flow out. With pressure depletion, methane that is adsorbed onto the organic matter tends to desorbs into the matrix pore space, where it can flow as free gas into the fracture and the wellbore. For the shale matrix, a modified convection-diffusion-adsorption model (CDAM) of the ¹³CH₄ and ¹²CH₄ isotopologues, which is based on Equation (2) in Cao et al. (2019), has been employed by considering real gas equation of state (EOS) and competitive adsorption:

(1)\{\begin{array}{l}V{\displaystyle \frac{\partial }{\partial t}}\left[{\displaystyle \frac{p_{\mathrm{1}}{M}_{\mathrm{1}}}{zRT}}{\phi }_{\mathrm{m}}+(\mathrm{1}-{\phi }_{\mathrm{m}}){\rho }_{\mathrm{r}}{\rho }_{\mathrm{s}\mathrm{t}\mathrm{d},\mathrm{1}}{V}_{\mathrm{L}}{\displaystyle \frac{b_{\mathrm{1}}{p}_{\mathrm{1}}}{\mathrm{1}+b_{\mathrm{1}}{p}_{\mathrm{1}}+b_{\mathrm{2}}{p}_{\mathrm{2}}}}\right]=-{\displaystyle \underset{V}{\iiint }\nabla N_{\mathrm{1}}\mathrm{d}V}\hfill \\ \\ \\ V{\displaystyle \frac{\partial }{\partial t}}\left[{\displaystyle \frac{p_{\mathrm{2}}{M}_{\mathrm{2}}}{zRT}}{\phi }_{\mathrm{m}}+(\mathrm{1}-{\phi }_{\mathrm{m}}){\rho }_{\mathrm{r}}{\rho }_{\mathrm{s}\mathrm{t}\mathrm{d},\mathrm{2}}{V}_{\mathrm{L}}{\displaystyle \frac{b_{\mathrm{2}}{p}_{\mathrm{2}}}{\mathrm{1}+b_{\mathrm{1}}{p}_{\mathrm{1}}+b_{\mathrm{2}}{p}_{\mathrm{2}}}}\right]=-{\displaystyle \underset{V}{\iiint }\nabla N_{\mathrm{2}}dV}\hfill \end{array},

where V is the equivalent volume of the fractured matrix; R is the gas constant; T is the absolute temperature; φ_m is the porosity of the matrix; ρ_r is the density of shale, which is equal to 2600 kg/m³ (Yu et al., 2021); z is the gas compressibility factor; V_L is the ultimate adsorbed amount of shale; ρ_std,1 and ρ_std,2 are the densities of ¹³CH₄ and ¹²CH₄ under standard conditions, respectively; and b₁ and b₂ are the Langmuir coefficients of ¹³CH₄ and ¹²CH₄, respectively.

The Langmuir coefficient is a function of the adsorption thermodynamic parameters (Xia and Tang, 2012):

(2)b=\exp \left({\displaystyle \frac{q}{RT}}+{\displaystyle \frac{\mathrm{\Delta }s}{R}}\right),

where q is the heat of adsorption and Δs is the standard entropy of adsorption. For ¹²CH₄, the values of q = 16 kJ/mol and Δs = −73 J/(mol_˙K) are used according to the methane adsorption isotherms of Longmaxi shale (Tian et al., 2016; Tang et al., 2017). The ratio (b₁/b₂) of the Langmuir coefficients of ¹³CH₄ and ¹²CH₄ is the fractionation factor between the adsorbate and bulk phases. The fractionation factor is dependent on temperature (van Hook, 1967):

(3)\alpha \mathrm{=}{\displaystyle \frac{b_{\mathrm{1}}}{b_{\mathrm{2}}}}=\exp \left({\displaystyle \frac{A}{T^{\mathrm{2}}}}-{\displaystyle \frac{B}{T}}\right),

where A = 95.17 is for the isotope fractionation of ¹³CH₄ and ¹²CH₄ (Xia and Tang, 2012), and B = 0.125 is calculated from the experimental data by Bruner et al. (1966).

For a fracture network, considering the communication of free gas between the fracture and the matrix, the governing flow equation for gas migration can be written as a continuity equation (Cao et al., 2019):

(4)\{\begin{array}{l}{\displaystyle \frac{\partial }{\partial t}}\left({\displaystyle \frac{p_{\mathrm{1}}{M}_{\mathrm{1}}}{zRT}}{\phi }_{\mathrm{f}}V\right)=-{\displaystyle \underset{V}{\iiint }\nabla U_{\mathrm{1}}\mathrm{d}V}+F_{\mathrm{1}}\hfill \\ {\displaystyle \frac{\partial }{\partial t}}\left({\displaystyle \frac{p_{\mathrm{2}}{M}_{\mathrm{2}}}{zRT}}{\phi }_{\mathrm{f}}V\right)=-{\displaystyle \underset{V}{\iiint }\nabla U_{\mathrm{2}}\mathrm{d}V}+F_{\mathrm{2}}\hfill \end{array},

where Ф_f is the porosity of the fracture; U₁ and U₂ are the mass fluxes of ¹³CH₄ and ¹²CH₄ in the fracture network, respectively; and F₁ and F₂ are the amounts of ¹³CH₄ and ¹²CH₄ from the matrix, respectively.

Owing to the differences in the sizes of pores in the matrix and fractures, the flow regimes vary. The significant slippage effect in the nanopores of the shale matrix results in a conversion from continuum flow to transition flow (Javadpour et al., 2007; Freeman et al., 2011), rendering the total mass flux through the matrix a combination of Knudsen diffusion and continuum flow (Cao et al., 2019).

(5)N_i=-{\displaystyle \frac{M_i\nabla p_i}{RT}}\left(D_{K,i}+{\displaystyle \frac{p{r}^{\mathrm{2}}}{\mathrm{8}\mu }}\right),

where N_i is the mass flux of methane in the matrix (i = 1 for ¹³CH₄ and i = 2 for ¹²CH₄), M_i is the molecular weight of ¹³CH₄ and ¹²CH₄, r is the nanopore radius, μ is the gas viscosity, p is the total pressure, and D_K,i is the Knudsen diffusion coefficient of ¹³CH₄ and ¹²CH₄.

(6)D_{K,i}={\displaystyle \frac{\mathrm{2}r}{\mathrm{3}}}\sqrt{{\displaystyle \frac{\mathrm{8}RT}{\mathrm{\pi }{M}_i}}}.

For fracture, only continuum flow is considered, and the mass flux is:

(7)U_i=-{\displaystyle \frac{M_i\nabla p_i}{RT}}\left({\displaystyle \frac{p{h}^{\mathrm{2}}}{\mathrm{12}\mu }}\right),

where U_i is the mass flux of methane in fractures (i = 1 for ¹³CH₄ and i = 2 for ¹²CH₄) and h is the width of the fracture.

3.4 Stable carbon isotope technique

According to Gao et al. (2017), the key to using methane carbon isotopes for constraining the production decline is to record the minimum δ¹³C₁ value of the produced gas and the corresponding time. Because of the uncertainty of the time when the minimum δ¹³C₁ value is recorded, which is caused by the complicated reservoir performance and reservoir behavior of the horizontal multi-stage fractured well, it is costly and unrealistic to monitor δ¹³C₁ of the produced gas until its minimum value is detected.

In this study, by considering Knudsen diffusion in nanopores and adsorption of organic matter, a dual porosity model was employed to simulate the carbon isotope variation of the produced gas. The δ¹³C₁ value of the produced gas at time t was calculated from the production rates of ¹³CH₄ and ¹²CH₄ (Cao et al., 2019).

(8)\delta {}^{\mathrm{13}}C{}_{\mathrm{1}}\left(t\right)=\left({\displaystyle \frac{q_{\mathrm{1}}\left(t\right)/q_{\mathrm{2}}\left(t\right)}{{\left({}^{\mathrm{13}}C/{}^{\mathrm{12}}C\right)}_{std}}}-\mathrm{1}\right)\times \mathrm{1000},

where q₁(t) and q₂(t) are the production rates of ¹³CH₄ and ¹²CH₄ at time t, respectively.

In this way, by studying a period when the produced gas exhibits δ¹³C1 variation (for example, three months or half a year), the δ¹³C1 variation during production is determined by matching the δ¹³C1 variation of the produced gas under given pressure conditions and parameters.

The gas recovery at time t is calculated as:

(9)R={\displaystyle \frac{Q_{\mathrm{1}}\left(t\right)+Q_{\mathrm{2}}\left(t\right)}{GIP}}\times \mathrm{100},

where Q₁(t) and Q₂(t) are the cumulative productions of ¹³CH₄ and ¹²CH₄ at time t, respectively. GIP is the initial gas content determined from the model at time zero.

The gas recovery corresponding to the latest monitored δ¹³C₁ is obtained from the relationship between the gas recovery and carbon isotopes from the model simulation. The cumulative production is the sum of all daily production data from the beginning to the day of the newest δ¹³C₁. The total reserves can be calculated from the cumulative production divided by the gas recovery. The time of the early hyperbolic decline transition to the later exponential decline is determined using the critical point where the total reserves derived from the integration of decline functions are equal to the value obtained using the isotope method.

4 Results and discussion

4.1 Production data and carbon isotope variations

Figure 3 presents the production curves and carbon isotope variations of the four wells. These four wells were produced in the output-restricted mode. The production histories were segmented by the replacement of nozzles of different sizes and the change in production policies, which were very different from those of typical shale gas wells in North America (Baihly et al., 2010). Each segment declined rapidly in the early stages, which was followed by a decreasing production decline; however, different decline rates were observed for each segment. The boosters were applied to enhance production at 500 days for well A, 470 days for well B, 455 days for well C, and 440 days for well D. Until October 31, 2019, the cumulative productions of wells A, B, C, and D were 3575 × 10⁴, 3645 × 10⁴, 3982 × 10⁴, and 3068 × 10⁴ m³, respectively.

The carbon isotope variations of the produced gas for these four wells declined gradually from –28‰ to –31‰ from September 2018 to May 2019 (Fig. 3). This is different from the reported carbon isotope variations of the produced gas in the Weiyuan and Changning shale gas fields in the Sichuan Basin. During shale gas production in the Weiyuan and Changning gas fields from 2012 to 2016, δ¹³C₁ displayed a slight increase of 1.3‰ (Zhang et al., 2018). In addition, this result is quite different from the carbon isotope variations of desorbed gas from shale core samples in well sites, which displayed a significant increase (Qin et al., 2017; Yang et al., 2017; Li et al., 2018). According to Li et al. (2020), various isotope variations are present during different shale gas production periods. The studied wells are in stage II, with adsorbed gas increasing in the produced gas according to the isotope variations (Fig. 3).

4.2 Mechanisms of isotope fractionation during shale gas production

Xia and Tang (2012) suggested that isotopic fractionation in adsorption/desorption systems is controlled by mass transport, and not by the partitioning between the adsorbate and free gas phases. Significant gas depletion in the fracture network during shale gas production caused the gas to flow along the pressure gradient from the matrix to the fracture, leaving behind heavier methane carbon isotopes because of the different effective diffusion coefficients of ¹³CH₄ and ¹²CH₄. The difference between the isotopic compositions of the instantaneous gas and the original gas without adsorption can be quantified using an empirical equation (Xia and Tang, 2012):

(10){\delta }_{\mathrm{i}\mathrm{n}\mathrm{s}}-{\delta }_{\mathrm{0}}=\mathrm{1000}\ln (D_{{\mathrm{13}}_{{\mathrm{CH}}_{\mathrm{4}}}}/D_{{\mathrm{12}}_{{\mathrm{CH}}_{\mathrm{4}}}})\mathrm{·(}\mathrm{1.21}+\ln (\mathrm{1}-R_{\mathrm{g}})),

where R_g is the gas recovery, D_{{\mathrm{13}}_{{\mathrm{CH}}_{\mathrm{4}}}}/D_{{\mathrm{12}}_{{\mathrm{CH}}_{\mathrm{4}}}} is the ratio of the effective diffusion coefficients of ¹³CH₄ and ¹²CH_4, and δ_ins and δ₀ are the instantaneous carbon isotope values of the produced gas and the initial value in the gas reservoir, respectively.

This equation suggests a trend of constant increase in instantaneous carbon isotope values of produced gas, which possibly can be used to simulate the carbon isotope variations of the produced gas in the Weiyuan and Changning shale gas fields (Zhang et al., 2018); however, this is not applicable to four wells from the Luzhou Shale Gas Field.

The ratio of effective diffusion coefficients between ¹³CH₄ and ¹²CH₄ is closely related to the total organic content (TOC) of the rock (Zhang and Krooss, 2001; Li et al., 2020). The higher the TOC value, the smaller is the D_{{{\mathrm{}}^{\mathrm{13}}}_{{\mathrm{CH}}_{\mathrm{4}}}}/D_{{{\mathrm{}}^{\mathrm{12}}}_{{\mathrm{CH}}_{\mathrm{4}}}} ratio. This is related to the strong adsorption of organic matter, as higher TOC values correlate with stronger adsorption of organic matter (Gasparik et al., 2014; Tian et al., 2016), which further results in a greater difference between the effective diffusion coefficients of ¹³CH₄ and ¹²CH₄. Moreover, Li et al. (2021) indicated that the Knudson number should be an intrinsic parameter controlling the D_{{{\mathrm{}}^{\mathrm{13}}}_{{\mathrm{CH}}_{\mathrm{4}}}}/D_{{{\mathrm{}}^{\mathrm{12}}}_{{\mathrm{CH}}_{\mathrm{4}}}} ratio, which is related to the apparent correlation between the TOC and D_{{{\mathrm{}}^{\mathrm{13}}}_{{\mathrm{CH}}_{\mathrm{4}}}}/D_{{{\mathrm{}}^{\mathrm{12}}}_{{\mathrm{CH}}_{\mathrm{4}}}} ratio. Fundamentally, in fluid dynamics, gas flow within the nanopore system in shale is discontinuous (Lilley and Sader, 2008; Dongari et al., 2009), which requires consideration of both molecule-molecule and molecule-wall interactions that always use slip models (Beskok and Karniadakis, 1999; Zhang et al., 2012; Taghavinejad et al., 2020) to describe the apparent permeabilities of fractured porous media (Freeman et al., 2011; Wang et al., 2020). The apparent permeability is closely related to the pore size and pressure; correspondingly, the effective diffusion coefficient varies during the depletion of pressure in a shale gas reservoir (Javadpour et al., 2007; Ziarani and Aguilera, 2012). Consequently, instantaneous carbon isotope fractionation during shale gas production is a complex process, and the trend can be easily omitted owing to the limited field data points. It is reasonable to assume that such a declining trend was not observed based on only six isotope values over 3.5 years in the Weiyuan and Changning shale gas fields (Zhang et al., 2018).

Recently, gas transportation simulation experiments of Longmaxi shales indicate four stages of methane carbon isotope variations during the pressure depletion from 10 MPa to 0.1 MPa (Li et al., 2020). The advection of free gas at high pressure, resulting in an almost constant carbon isotope value, is the main feature of the first stage. In the second stage, adsorbed ¹²CH₄ is preferentially desorbed into free gas as the pressure decreases, leading to isotopically lighter gas flowing out of the shale. In the third stage, the gas is mainly produced from adsorbed gas, and owing to the fractionation mechanism of gas adsorption/desorption, it quickly becomes isotopically heavier because the desorbed gas becomes heavier. Minimum methane carbon isotope values are observed between the second and third stages, which also appeared during the production of a Barnett shale gas well (Gao et al., 2017). Compared with the above results, the methane carbon isotope variations of the four wells in Fig. 3 are distributed in stage II, with a gradual decrease in δ¹³C₁ values. Li et al. (2020) claimed that the amount of adsorbed gas increases in stage II, which contradicts the estimation of adsorbed gas in place. If the produced gas becomes 100% adsorbed gas in stage III, the gas recovery at the point of the minimum isotope value is greater than 50% (Li et al., 2020); this is because the adsorbed gas accounts for approximately 20% to 40% of the reserves in the Longmaxi shale (Zou et al., 2021). This presents a large difference due to the recovery factors being less than 5% at the minimum point (Gao et al., 2017) and it disagrees with the US shale gas production that exhibits recovery factors ranging from 20% to 30% (EIA, 2015).

When the gas recovery and the ratio of adsorbed to free gas in shale gas wells are considered, it can be concluded that only a small portion of free gas is produced. The difference between the effective diffusion coefficients or apparent permeabilities of ¹³CH₄ and ¹²CH₄ within long tortuous nano-pathways in the shale matrix leads to carbon isotope fractionation during the gas production process. The desorbed ¹²CH₄ molecules have a higher flow velocity than ¹³CH₄ molecules in free gas (Zhang and Krooss, 2001), resulting in a decrease in δ¹³C₁ values. When isotopically heavier free gas exceeds isotopically lighter desorbed gas with gas production, isotope rollover occurs. The reversal is related to the nanopore structure and the adsorbed gas ratio; more complex nanopore structures and greater adsorbed gas ratios result in stronger isotope reversal.

The carbon isotope variations in shale gas during production are related to the well production capacity. Once the pore pressure is reduced, isotopically lighter gas desorbs and becomes free gas, which alters the δ¹³C₁ value of the producing gas. A lower decrease in carbon isotope values corresponds to a greater production capacity. In light of this statement, wells B and C exhibit a greater production capacity than that of wells A and D (Fig. 3), which is in line with the actual cumulative production.

4.3 Decline curve analysis for early production

The production histories were quite erratic (Fig. 3), which impractically made the use of DCA. However, the cumulative production-time relationship under hyperbolic decline was regular and could be used in the DCA. Discrete least squares approximation was applied to determine the initial production rate q_i, initial decline rate D_i, and decline exponent b.

The cumulative productions before May 30, 2019, which corresponded to the latest δ¹³C₁ values of the produced gas, were used for the fitting calculation. Figure 4 presents the results from the cumulative curves for these four wells using the hyperbolic equation (Table 1). The actual production curves were undetermined owing to erratic production histories. However, the modeling results agree well with the hyperbolic decline. The cumulative productions obtained by matching the hyperbolic decline were very close to the actual cumulative productions, with errors of less than 4%. Table 1 lists the fitted parameters of hyperbolic decline. Noticeably, all the decline exponents are greater than 1, which leads to an unbounded and unrealistic prediction of the total reserve if the Arps hyperbolic decline equation is used.

4.4 Total reserve from carbon isotope simulation

The ¹³CH₄ and ¹²CH₄ molecules have different flow velocities within the nanopore system of the shale matrix, owing to the Knudsen diffusion and adsorption effect related to their respective molecular weights. Significant carbon isotope fractionation occurs during different stages of gas production (Gao et al., 2017; Li et al., 2020). To simulate the carbon isotope ratio variations during gas production from a fractured shale gas reservoir, a bicomponent gas CDAM (Cao et al., 2019) was proposed to calculate the ¹³CH₄ and ¹²CH₄ methane released from the shale matrix (Knudsen diffusion in nanopores and adsorption of organic matter coexist). Then, this bicomponent gas convection model was employed to calculate the gas production and carbon isotope variation from the fracture network in the shale gas reservoir.

For the four studied wells, the initial carbon isotopic signature was set to −28‰ by extrapolating backward to zero time (Fig. 3). The parameters of the formation pressure and temperature are provided according to our knowledge of the Longmaxi shale gas. The boundary pressure was set as the average value based on the production pipeline pressure. The parameters of pore size, porosity, and adsorption capacity of the Longmaxi shale are average values based on sample measurements. The size of the equivalent volume is empirically given. An iterative algorithm, which is also known as linear quadratic estimation, has been used to simulate the isotope variations. Meanwhile, a damped least-squares method is used to determine the optimal results by fine-tuning the facture parameters in Table 2.

The relationship between carbon isotope values (Eq. (8)) and gas recovery (Eq. (9)) for each well was used to determine the recovery at the latest point of the actual carbon isotope value of the produced gas. For well A, the gas recovery on May 30, 2019 was approximately 15.6% based on the carbon isotope simulation results. At that time, the cumulative production from the production rate data was 2985 × 10⁴ m³. Thus, the total reserve for well A was estimated to be 19134 × 10⁴ m³. The gas recoveries on May 30, 2019 were approximately 15.3%, 14.8%, and 14.4% for wells B, C, and D, respectively, which translate to cumulative productions of 2787 × 10⁴, 2985 × 10⁴, and 2475 × 10⁴ m³, respectively. The total reserves for these three wells were thus estimated to be 18216 × 10⁴, 20169 × 10⁴, and 17187 × 10⁴ m³, respectively.

4.5 Prediction of production decline

Prediction from the hyperbolic decline equation yielded a greater decline exponent than 3 (Table 1), leading to unrealistic total reserves. At a certain time point (t_c), early hyperbolic decline would decrease exponentially, thereby representing later stages. Determination of t_c for a given well required that the total reserve calculated from the integration of the production decline curve was equal to the value derived from the isotope fractionation model. The total reserve from the integration of the production decline curve was divided into two parts: one from the hyperbolic decline curve (0–t_c) and the other from the exponential decline curve (t_c-infinite).

(11)\{\begin{array}{l}{\displaystyle \underset{\mathrm{0}}{\overset{t_{\mathrm{c}}}{\int }}{\displaystyle \frac{q_{\mathrm{i}}}{{\left(\mathrm{1}+b{D}_{\mathrm{i}}t\right)}^{\mathrm{1}/b}}}\mathrm{d}t={\displaystyle \frac{q_{\mathrm{i}}}{\left(\mathrm{1}-b\right)D_{\mathrm{i}}}}\left[\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{{\left(\mathrm{1}+b{D}_{\mathrm{i}}{t}_{\mathrm{c}}\right)}^{\mathrm{1}/b-\mathrm{1}}}}\right]}\hfill \\ {\displaystyle \underset{t_{\mathrm{c}}}{\overset{+\infty }{\int }}{q}_{\mathrm{i}}{e}^{-Dt}\mathrm{d}t={\displaystyle \frac{q_{\mathrm{i}}}{D}}{e}^{-D{t}_{\mathrm{c}}}}\hfill \end{array}.

The convergence (t_c) of the hyperbolic decline and exponential decline was observed on days 518, 324, 315, and 301 for wells A, B, C, and D, respectively (Fig. 5). A hyperbolic decline was applied before the convergence, and an exponential decline was used afterwards. Based on the hybrid decline curve, future production for a given well can be forecasted, and the EUR can be determined during a limited time.

The hyperbolic decline overestimates the EUR for these four shale gas wells (Fig. 5). The duration of the hyperbolic decline was short compared to the long life of the shale gas well, which implied that the fracture-dominated flow only existed in a short period of time when gas flow penetrated into the matrix and the area of production continued to increase. When the matrix-dominated flow prevailed, the production rate decreased significantly because the matrix permeability was several orders of magnitude smaller than the fracture permeability. However, the decline rate is very small and the production rate decreases very slowly, which guarantees a long life for shale gas production.

Based on the decline curves in Fig. 5, the predicted cumulative production until October 31, 2019 for well A was 3639 × 10⁴ m³, which was 1.8% greater than the actual cumulative production of 3575 × 10⁴ m³. For wells B, C, and D, the predicted cumulative production until October 24, 2019 was 3545 × 10⁴, 3823 × 10⁴, and 3241 × 10⁴ m³, respectively; these values are similar to the actual cumulative productions for each well.

5 Conclusions

In this study, we have further improved a novel tool that was proposed by Gao et al. (2017), with the aim of interpreting gas stable carbon isotope values for predicting the production decline curves of shale gas wells. A large number of sequential real-time methane carbon isotope values of the produced gas is a prerequisite for this method. Combined with the variation in the carbon isotope values of the produced gas, the evolution of the methane carbon isotope during production could be simulated by using a numerical simulator composed of a CDA model, while taking into account Knudsen diffusion and competitive adsorption in the matrix as well as the convection model for fractures in the ¹³CH₄ and ¹²CH₄ isotopologues; this is aimed at determining the relationship between isotope variation and gas recovery. The two-segment production decline curve of a transition from a hyperbolic to an exponential equation is determined by constraining the total reserves from the integration of the production curve and the isotope fractionation model.

The field study in the Longmaxi shale gas wells sets a good example for the future application of gas isotope tools in the determination of production curves. The total reserves for the four developing wells range from 1.72 × 10⁸ to 2.02 × 10⁸ m³, as derived from the simulation results of the isotope fractionation model. By 0.82–0.89 year, the production rates begin to decline after the exponential equation was used for the EUR forecast. This is a trial that combines the isotope variation of the produced gas to improve the confidence in production prediction. Our ability to model controls on the produced gas isotopes is limited by the incomplete reservoir data and the difficulty in validating models; this is because of the limited duration for which the isotopic data are collected from the wells. The two-segment Arps model may be inappropriate for use in extremely low-permeability reservoirs. Recent alternative models (e.g., SEDM and PLE) for shale gas wells may be more appropriate for these applications. As the physics of the production mechanism becomes better known, both applications of isotope variation and prediction models can be further developed to provide increasingly more reliable forecasts of the future performance of shale gas reservoirs.

References

[1] Arps J J (1945). Analysis of decline curves. Transactions of the AIME, 160(01): 228–247

[2] Baihly J D, Altman R M, Malpani R, Luo F (2010). Shale gas production decline trend comparison over time and basins. In: SPE Annual Technical Conference and Exhibition, SPE-135555-MS

[3] Beskok A, Karniadakis G E (1999). A model for flows in channels, pipes, and ducts at micro and nano scales. Microscale Thermophysical Engineering, 3(1): 43–77

[4] Bruner F, Cartoni G P, Liberti A (1966). Gas chromatography of isotopic molecules on open tubular columns. Anal Chem, 38(2): 298–303

[5] Cao G H, Zhang H X, Jiang W B, Wu S, Zhu D, Lin M (2019). A new gas-content-evaluation method for organic-rich shale using the fractionation of carbon isotopes of methane. SPE J, 24(06): 2574–2589

[6] Dongari N, Sambasivam R, Durst F (2009). Extended Navier-Stokes equations and treatments of micro-channel gas flows. J Fluid Sci Tech, 4(2): 454–467

[7] Currie S M, Ilk D, Blasingame T A (2010). Continuous estimation of ultimate recovery. In: SPE Unconventional Gas Conference, SPE-132352-MS

[8] Duong A N (2011). Rate-decline analysis for fracture-dominated shale reservoirs. SPE Reservoir Eval Eng, 14(03): 377–387

[9] Freeman C M, Moridis G J, Blasingame T A (2011). A numerical study of microscale flow behavior in tight gas and shale gas reservoir systems. Transp Porous Media, 90(1): 253–268

[10] Gao L, Wu S, Deev A, Olson R, Mosca F, Zhang S, Ni Y, Qu Q, LaFollette R, Chen G, Tang Y (2017). The gas isotope interpretation tool: a novel method to better predict production decline. AAPG Bull, 101(08): 1263–1275

[11] Gasparik M, Bertier P, Gensterblum Y, Ghanizadeh A, Krooss B M, Littke R (2014). Geological controls on the methane storage capacity in organic-rich shales. Int J Coal Geol, 123: 34–51

[12] Ilk D, Anderson D M, Stotts W J, Mattar L, Blasingame T A A (2010). Production-date analysis—challenges, pitfalls diagostics. SPE Reservoir Eval Eng, 13(3): 538–552

[13] Ilk D, Rushing J A, Perego A D, Blasingame T A (2008). Exponential vs. hyperbolic decline in tight gas sands: understanding the origin and implications for reserve estimates using Arps’ decline curves. In: SPE Annual Technical Conference and Exhibition, SPE-116731-MS

[14] Javadpour F, Fisher D, Unsworth M (2007). Nanoscale gas flow in shale gas sediments. J Can Pet Technol, 46(10): 55–61

[15] Joshi K, Lee W J (2013). Comparison of various deterministic forecasting techniques in shale gas reservoirs. In: SPE Hydraulic Fracturing Technology Conference, SPE-163870-MS

[16] Kanfar M S, Wattenbarger R A (2012). Comparison of empirical decline curve methods for shale wells. In: SPE Canadian Unconventional Resources Conferences, SPE-162648-MS

[17] Lee W J, Sidle R (2010). Gas-reserves estimation in resource plays. SPE Econ Manage, 2(2): 86–91

[18] Li K, Meng Z Y, Ji J, Zheng X W, Zhang Q, Zhou W (2018). Characteristics and influencing factors of desorption gas in Wufeng-Longmaxi formations in Fuling area, Sichuan Basin. Petrol Geo Exper, 40(1): 90–96 (in Chinese)

[19] Li W B, Lu S F, Li J Q, Zhang P F, Wang S Y, Feng W J, Wei Y B (2020). Carbon isotope fractionation during shale gas transport: mechanism, characterization and significance. Sci China Earth Sci, 63: 674–689

[20] Li W B, Lu S F, Li J Q, Wei Y B, Feng W J, Zhang P F, Song Z J (2021). Geochemical modeling of carbon isotope fractionation during methane transport in tight sedimentary rocks. Chem Geol, 566: 120033

[21] Lilley C R, Sader J E (2008). Velocity profile in the Knudsen layer according to the Boltzman equation. Proc R Soc Lond A Math Phy Sci, 464(2096): 2015–2035

[22] Ma X H, Xie J, Yong R, Zhu Y Q (2020). Geological characteristics and high production control factors of shale gas reservoirs in Silurian Longmaxi Formation, southern Sichuan Basin, SW China. Pet Explor Dev, 47(5): 901–915

[23] Mahmoud O, Ibrahim M, Pieprzica C, Larsen S (2018). EUR prediction for unconventional reservoirs: state of the art and field case. In: SPE Trinidad and Tobago Section Energy Resources Conference, SPE-191160-MS

[24] Meyet M, Dutta R, Burns C (2013). Comparison of decline curve analysis methods with analytical models in unconventional plays. In: SPE Annual Technical Conference and Exhibition, SPE-166365-MS

[25] Mishra S (2012). A new approach to reserves estimation in shale gas reservoirs using multiple decline curve analysis models. In: SPE Eastern Regional Meeting, SPE-161092-MS

[26] Moinfar A, Erdle J C, Patel K (2016). Comparison of numerical vs analytical models for EUR calculation and optimization in unconventional reservoirs. In: SPE Low Perm Symposium, SPE-180209-MS

[27] Okouma V, Symmons D, Hosseinpour-Zonoozi N, Ilk D, Blasingame T A (2012). Practical considerations for decline curve analysis in unconventional reservoirs—application of recently developed time-rate relations. In: SPE Hydrocarbon, Economics, and Evaluation Symposium, SPE-162910-MS

[28] Qin H, Fan X J, Liu M, Hao J Y, Liang B (2017). Carbon isotope reversal of desorbed gas in Longmaxi shale of Jiaoshiba area, Sichuan Basin. Petrol Res, 2(2): 169–177

[29] Robertson S (1988). Generalized Hyperbolic Equation. Society of Petroleum Engineers, SPE-18731-MS

[30] Seshadri J N, Mattar L (2010). Comparison of power law and modified hyperbolic decline methods. In: SPE Canadian Unconventional Resources and International Petroleum Conference, SPE-137320-MS

[31] Society of Petroleum Evaluation Engineers (SPEE) (2016). Monograph 4: Estimating ultimate recovery of developed wells in unconventional reservoirs

[32] Taghavinejad A, Sharifi M, Heidaryan E, Liu K, Ostadhassan M (2020). Flow modeling in shale gas reservoirs: a comprehensive review. J Nat Gas Sci Eng, 83: 103535

[33] Tan L, Zuo L H, Wang B B (2018). Methods of decline curve analysis for shale gas reservoirs. Energies, 11(3): 552–569

[34] Tang X, Ripepi N, Stadie N P, Yu L J (2017). Thermodynamic analysis of high pressure methane adsorption in Longmaxi Shale. Fuel, 193: 411–418

[35] Tian H, Li T F, Zhang T W, Xiao X M (2016). Characterization of methane adsorption on overmature Lower Silurian–Upper Ordovician shales in Sichuan Basin, southwest China: experimental results and geological implications. Int J Coal Geol, 156: 36–49

[36] U.S. Energy Information Administration (EIA) (2015). Technically Recoverable Shale Oil and Shale Gas Resources: Brazil.

[37] Valko P P (2009). Assigning value to stimulation in the Barnett shale: a simultaneous analysis of 7000 plus production histories and well completion records. In: SPE Hydraulic Fracturing Technology Conference, SPE-119369-MS

[38] van Hook W A (1967). Isotope effects on vaporization from the adsorbed state: methane system. J Phys Chem, 71(10): 3270–3275

[39] Wang L, Torres A, Xiang L, Fei X, Naido A, Wu W (2015). A technical review on shale gas production and unconventional reservoirs modeling. Nat Resour, 6(03): 141–151

[40] Wang H, Chen L, Qu Z, Yin Y, Kang Q, Yu B, Tao W Q (2020). Modeling of multi-scale transport phenomena in shale gas production—a critical review. Appl Energy, 262: 114575

[41] Xia X Y, Tang Y C (2012). Isotope fractionation of methane during natural gas flow with coupled diffusion and adsorption/desorption. Geochim Cosmochim Acta, 77: 489–503

[42] Yang Z H, Wei Z H, He W B, Fan M, Yu L J, Xu E S, Qian M H (2017). Characteristics and significance of on-site gas desorption from Wufeng-Longmaxi shales in southeastern Sichuan Basin. Nat Gas Geosci, 28(1): 156–163 in Chinese)

[43] Yu S Y, Wang X X, Li S H, Liu Y M, Xiao L M, Liu X, Zhao W, Zhang J E (2021). Impact of geological factors on marine shale gas enrichment and reserve estimation: a case study of Jiaoshiba area in Fuling gas field. Geofluids, 2021: 6637360

[44] Zhang H, Rietz D, Cagle A, Cocco M, Lee J (2016). Extended exponential decline curve analysis. J Nat Gas Sci Eng, 36: 402–413

[45] Zhang M J, Tang Q Y, Cao C H, Lv Z, Zhang T, Zhang D, Li Z, Du L (2018). Molecular and carbon isotopic variation in 3.5 years shale gas production from Longmaxi Formation in Sichuan Basin, China. Mar Pet Geol, 89: 27–37

[46] Zhang T W, Krooss B M (2001). Experimental investigation of the carbon isotope fractionation of methane during gas migration by diffusion through sedimentary rocks at elevated temperature and pressure. Geochim Cosmochim Acta, 65(16): 2723–2742

[47] Zhang W M, Meng G, Wei X (2012). A review on slip models for gas microflows. Microfluid Nanofluidics, 13(6): 845–882

[48] Ziarani A S, Aguilera R (2012). Knudsen’s permeability correction for tight porous media. Transp Porous Media, 91(1): 239–260

[49] Zou C, Zhao Q, Cong L Z, Wang H Y, Shi Z S, Wu J, Pang S Q (2021). Development progress, potential and prospect of shale gas in China. Natural Gas Indust, 41(1): 1–14 in Chinese)

[50] Zuo L H, Yu W, Wu K (2016). A fractional decline curve analysis model for shale gas reservoirs. Int J Coal Geol, 163: 140–148

RIGHTS & PERMISSIONS

Higher Education Press