Introductions

Natural gas resource is an important energy and strategic material, which is very important to ensure the sustainable and healthy development of economy and society. With the sustained and rapid development of China's national economy, the natural gas industry has achieved leapfrog development. Compared with 2000, the total consumption of natural gas in 2015 reached 193.175 billion cubic meters, an increase of nearly seven times. According to the National Energy Administration, the proportion of natural gas in primary energy consumption is expected to increase to 35.7% by 2022 (Dong et al. 2017; Lin et al. 2021). Therefore, it is necessary to study the distribution of natural gas resources to reasonably plan the consumption pattern of regional natural gas resources. Among them, to strengthen the study of natural gas reserves, growth law is conducive to the reasonable formulation of future natural gas resources development plan. Furthermore, it can effectively promote the development of national industry (Hu et al. 2015; Jia et al. 2014).

At present, more than 80% of the total proven natural gas in China is distributed in the five basins of Sichuan, Ordos, Qaidam, Tarim and Heyingyiqiong (Zheng et al. 2018; Gong et al. 2019). After more than 60 years of exploration and development, Sichuan Basin has made great progress and formed a certain scale of reserves, which is attributed to the progress in the prediction of natural gas reserves in Sichuan Basin (Lu et al. 2017). Therefore, strengthening the study on the growth trend of natural gas reserves in Sichuan Basin can promote the further exploration and development of natural gas in Sichuan Basin in the future, which is conducive to the development of natural gas industry in China.

However, the distribution of oil and gas resources in Sichuan Basin is relatively dispersed, which increases the difficulty of natural gas reserves prediction (Liu et al. 2015). Breakthroughs in all fields of natural gas exploration and development will lead to a rapid increase in new proven reserves in the short term. However, the tendency of natural gas reserves to decline after rapid growth is difficult to predict. Therefore, the prediction of natural gas reserves in Sichuan Basin needs to be carried out deeply, so as to realize deeper exploration and development of natural gas resources (Davis and Harbaugh 1983).

Based on statistical research on the actual dynamic data of a large number of oil and gas field developments, many scholars have derived and proposed a series of natural gas reserves prediction models, such as the life cycle method (Ramezanianpour and Sivakumar 2017), the gray system method (Karahan and Ayvaz., 2006) and the neural network method (Hubbert 1949), among others. Most oil and gas field reserves and production growth trends conform to the peak model theory. In 1956, Hubbert predicted the future US oil supply by fitting the historical oil production curve. The prediction is based on two assumptions: First, that production will eventually decline exponentially and, second, that the area enclosed by the curve must be equal to the ultimate recoverable oil resources of the USA. Using the industry’s universally predicted US ultimate recoverable resources (150–200 billion barrels), Hubbert predicted that US oil production will reach its peak during 1965–1971 (Kaufmann and Cleveland 2001; Abdideh and Fatha Ba 2013).

At present, the research study on the natural gas reserves growth trend prediction is relatively extensive. Liu H (2017) used the Poisson model to predict the proven reserves of natural gas in the Sichuan Basin in the next 30 years. The prediction results show that the growth trend of proved natural gas reserves presents a significant cyclical characteristic, with a peak year roughly every 6 years. The final proved reserves of natural gas in Sichuan Basin will reach about 25,000 × 108m3 in 2025 (Liu et al. 2017). Yang Y (2021) applied the improved Weng's model and Weibull model in the study on natural gas production prediction in Sichuan Basin. The prediction results show that the annual gas output predicted by the two models in 2030 is \(\left( {786 - 856} \right) \times 10^{8} {\text{m}}^{3}\), which conform the development situation of natural gas in Sichuan Basin (Yang et al. 2021). Zhi D (2021) analyzed the potential and exploration status of natural gas resources in Junggar Basin. The HCZ model is used to predict the growth trend of proven natural gas reserves throughout the life cycle. The results show that the peak of proved reserves will occur in 2027, the annual proved reserves will reach 268.49 × 108m3, and the corresponding proved rate will be 25.8% in 2027. The natural gas reserves prediction results provide a theoretical basis for natural gas exploration and development (Zhi et al. 2021). Peng W (2019) analyzed the exploration and development history of the Ordos Basin. They combined the multi-peak Gauss model with the multi-peak characteristics of oil reserves and production growth in the Ordos Basin to fit the historical curve of oil reserves and production in Ordos Basin and ultimately predict the growth trend of reserves and production in the basin from 2006 to 2030 (Peng et al. 2019).

Obviously, the research work of natural gas reserves prediction is mostly carried out by using a single model, which can easily cause large prediction error. Due to the complex geology in the Sichuan Basin, the discovery of natural gas reserves has obvious fluctuations. Therefore, multiple methods can be combined together prediction methods to predict gas reserves. In this manuscript, 13 multi-cycle peaks are screened out from the growth curve of natural gas reserves through the original multi-cycle peak judgment standard. The history of natural gas exploration and development is combined with the trend of reserves change. And the historical change rule of natural gas reserves is comprehensively analyzed. Metabolism modified GM(1,3) model is used to predict the parameters of multi-cycle model, which improves the accuracy of parameter prediction results. The parameter prediction results are substituted into the multi-cycle Hubbert and Gauss models to predict the natural gas reserves. The correlation coefficient is used to screen the reserves prediction results with higher accuracy. The research results in the manuscript provide a theoretical basis for the exploration and development of natural gas in Sichuan Basin.

Natural gas reserves prediction theory

Overview of natural gas resources in Sichuan Basin

Sichuan Basin is rich in conventional gas, shale gas, tight gas and other oil and gas resources. It is one of the basins with the greatest potential for oil and gas reserves growth in China. Sichuan Basin belongs to a second-order tectonic unit on the “Yangtze quasi-platform,” which is a large sedimentary basin with stable subsidence inside the platform. Longitudinally, Sichuan Basin is characterized by multilayer system and multi-cycle (Gao et al. 2017).

Since 1956, the natural gas exploration and development process in Sichuan Basin has experienced more than 60 years and is currently in a period of rapid reserves growth (Guo et al. 2016). The exploration and development process is mainly divided into the following four stages.

  1. (1)

    The first stage: natural gas industry initial stage (1956–1977): the exploration of surface structure and fractured gas reservoirs was the main focus. In 1950 and 1960, two large-scale oil production projects are carried out, and no high-yielding large oil fields are found. The oil and gas exploration strategy was changed from “oil-dominated” to “oil–gas-dominated” (Lei 2010).

  2. (2)

    The second stage: steady natural gas growth stage (1978–2004): the discovery of the Xiangguosi carboniferous pore gas reservoir in 1977 is a sign of a major change in the exploration of Sichuan Basin. The exploration idea of targeting large- and medium-sized gas fields with fractured and pore reservoirs as the object is clarified. The carboniferous system is mainly explored, and the second and Triassic systems are also explored. A number of large- and medium-sized gas fields represented by the Wuwuti-carboniferous system gas field and the Feixianguan formation gas field in Dukouhe are discovered (Xian et al. 2017).

  3. (3)

    The third stage: rapid natural gas growth stage (2005–2010): since 2005, under the guidance of the lithologic reservoirs theory, great efforts have been maken to integrate exploration and development by targeting large- and medium-sized gas fields and combining Marine and land facies. The Xujiahe Formation and Pertriassic reef-shoal in the basin have achieved fruitful results in natural gas exploration. Sichuan Basin has entered the peak period of natural gas storage and production growth (Zhang et al. 2007).

  4. (4)

    The fourth stage: natural gas development and expansion stage (2011—present): taking strategic discovery as the goal, risk exploration is carried out in Sichuan Basin. It has achieved a breakthrough in Sinian exploration and a historic new discovery of Longwangmiao Formation. The ancient theory of carbonate reservoir formation with rifting as the core has been formed. This theory has contributed to the rapid exploration, efficient construction and steady expansion of Anyue extra-large gas fields in natural gas exploration and development (Chen et al. 2014).

Throughout the exploration and development of Sichuan Basin for more than 60 years, it is obvious that the annual proven reserves of natural gas vary with many peaks. As shown in Fig. 1, the annual proved reserves of the Sichuan Basin show multiple ups and downs over time, with obvious characteristics of multi-cycle. From the beginning of exploration and development in Sichuan Basin to 2018, the accumulative proved geological reserves of natural gas in Sichuan Basin are 28,747.55 × 108m3. And the unproved geological resources of natural gas are 249,336.03 × 108m3, with a proven rate of only 10.34%. The annual average proved reserves began to increase rapidly in 2011. At present, the basin is in the early stage of exploration and development with rapidly increasing reserves.

Fig. 1
figure 1

The growth trend of natural gas reserves in Sichuan Basin from 1956 to 2018

Full size image

Multi-cycle prediction model of natural gas reserves

At present, the prediction methods of oil and gas reserves mainly include analogy method, volume method, material balance method, production decline method, extrapolation prediction method, gray system prediction method (Cavallo 2002). Q82 Among them, based on the past reserves data, there are two basic methods for the prediction of future new reserves: extrapolation prediction method and gray system prediction method (Seifabad and Ehteshami 2013).

The extrapolation prediction method transforms and deduces the mathematical formula of the model according to its properties. The predicted reserves is obtained from known model parameters. The actual data of natural gas exploration and development area are analyzed. Through strict modeling, the extrapolation prediction method is used to predict the future development trend and evolution process of natural gas exploration and development areas (Zhang and Rong 2018). Accurate model parameters are the basis of accurate prediction by this method. Therefore, it is necessary to find other methods to predict the model parameters separately.

Gray system prediction method can predict the system containing both known information and unknown information. Gray prediction can effectively find some laws between data and predict the future trend of data change. However, gray prediction takes the predicted data as the prediction sample for a new round of prediction, which can easily cause the prediction results to deviate from the original law (Luo et al. 2020; Sheng et al. 2020). Therefore, it is necessary to improve the gray prediction method effectively to ensure the stability of the prediction results.

Common extrapolation prediction methods include logistic model, Weng cycle model, Hubbert model, Gauss model, Weibull model and HCZ model (Wang et al. 2015; Maggio and Cacciola 2009). The newly proven reserves of natural gas resources in Sichuan Basin show multi-peak with obvious multi-cycle state over time (Fig. 1). Therefore, the multi-cycle Hubbert model and multi-cycle Gauss model are both suitable for the study of the reserves growth rule of this basin.

Hubbert model of natural gas reserves prediction

Hubbert model is one of the multi-cycle reserves prediction models. The change process of the curve is from a gradual increase at the beginning, then to a stable period at the peak and finally to a rapid decline to complete resource consumption (Mahetaji et al. 2020). According to the single-cycle Hubbert model, the relationship between annual newly added proved reserves and time is as follows:

$$ N = \frac{{2N_{m} }}{{1 + \cosh \left[ {b\left( {t - t_{m} } \right)} \right]}} $$
(1)

where N is new proved reserves per year, (\(10^{8} m^{3} /a\)); \(N_{m}\) is the peak of annual newly added proved reserves, (\(10^{8} m^{3} /a\)); t is the reserves submission time; \(t_{m}\) is the peak occurrence time of annual newly added proved reserves; and b is the model parameter, which can represent the ratio of height to width of peak to some extent; it is not exactly the same as the slope.

Assume that there is an important prediction theory innovation or technological progress in the natural gas exploration and development work. Natural gas reserves will again show the cyclical change of growth-level decline. If the number of technological progress is more, there may be multiple cycle changes, that is, multi-peak phenomenon. The multi-cycle Hubbert model can be expressed as:

$$ N = \mathop \sum \limits_{{i = 1}}^{k} \frac{{2N_{{{\text{mi}}}} }}{{1 + \cosh \left[ {b_{i} \left( {t - t_{{{\text{mi}}}} } \right)} \right]}} $$
(2)

where k is the total number of cycles; i is the number of cycles, represents the parameter subscript of some cycle peak. For example, when \(i = 3\), \(N_{{{\text{mi}}}}\) represents the peak of annual newly added proved reserves in the third cycle.

Gauss model of natural gas reserves prediction

Hubbert model and Gauss model have been widely used in the research of gas storage and production prediction. The curve shape of Gauss model is similar to that of Hubbert curve, both of which are symmetrical models. However, the peak time of Gauss model is relatively late, and the curve tends to be flat. The curve increases gently at first, then rises rapidly, then increases slowly and then reaches a peak and begins to decline in a symmetrical form. This model is more suitable for predicting the growth law of natural gas reserves with relatively gentle changes (Mfoubat and Zaky 2019; H Refsnes., Diaz M., Stanko M. 2019).

The principle of Gauss model is the same as that of Hubbert model, which is also based on the growth curve obtained by life cycle. According to the single-cycle Gauss model, the relationship between annual newly added proved reserves and time is as follows:

$$ N = N_{m} \cdot e^{{ - \left( {t - t_{m} } \right)^{2} /2s^{2} }} $$
(3)

where s is the model parameter, which can represent the fluctuation of peak to some extent, and it is not exactly the same thing as the standard deviation.

The multi-cycle Gauss model can be expressed as:

$$ N = \mathop \sum \limits_{{i = 1}}^{k} N_{{{\text{mi}}}} \cdot e^{{ - \left( {t - t_{{{\text{mi}}}} } \right)^{2} /2s_{i} ^{2} }} $$
(4)

where k is the total number of cycles and i is the number of cycles, which represents the parameter subscript of some cycle peak. For example, when \(i = 3\), \(N_{{{\text{mi}}}}\) represents the peak of annual newly added proved reserves in the third cycle.

Multi-cycle parameter prediction model based on GM(1,n) gray theory

Hubbert and Gauss models can be used to analyze the multi-cycle growth law of natural gas reserves, on the premise that the values of the multi-cycle model parameters must be known (Arnaud and Reinout 2019; Ebrahimi and Cheshme 2015). Therefore, it is necessary to select effective mathematical methods to carry out the prediction research of multi-cycle parameters. The GM(1,n) gray prediction model can be applied to the numerical prediction.

The GM(1,n) model means that for n variables \(x_{1}\), \(x_{2}\)…\(x_{n}\) is a gray model based on first-order differential equation. Assuming that

$$ \left\{ {\begin{array}{*{20}c} {x_{i} ^{{\left( 0 \right)}} = \left[ {x_{i} ^{{\left( 0 \right)}} \left( 1 \right),\;x_{i} ^{{\left( 0 \right)}} \left( 2 \right) \cdots x_{i} ^{{\left( 0 \right)}} \left( p \right)} \right]} \\ {i = 1,\;2 \ldots n} \\ \end{array} } \right. $$
(5)

Add up \(x_{i} ^{{\left( 0 \right)}}\) to generate a new sequence

$$ \left\{ {\begin{array}{*{20}c} {x_{i} ^{{\left( 1 \right)}} = \left[ {x_{i} ^{{\left( 1 \right)}} \left( 1 \right),\;x_{i} ^{{\left( 1 \right)}} \left( 2 \right) \ldots x_{i} ^{{\left( 1 \right)}} \left( p \right)} \right]} \\ {x_{i} ^{{\left( 1 \right)}} \left( k \right) = \mathop \sum \limits_{{m = 1}}^{k} x_{i} ^{{\left( 0 \right)}} \left( m \right)} \\ \end{array} } \right. $$
(6)

where \(k = 1,\;2 \ldots p\;;\;i = 1,\;2 \ldots n\).Sequence x i (1) meet the following first-order linear differential equation

$$ \frac{{dx_{1} ^{{\left( 1 \right)}} }}{{dt}} + ax_{1} ^{{\left( 1 \right)}} = b_{1} x_{2} ^{{\left( 1 \right)}} + b_{2} x_{3} ^{{\left( 1 \right)}} + \cdots b_{{n - 1}} x_{n} ^{{\left( 1 \right)}} $$
(7)

The discrete format of the above equation can be rewritten as

$$ {\text{x}}_{1} ^{{\left( 0 \right)}} \left( {{\text{k}} + 1} \right) + \frac{{\text{a}}}{2}\left[ {{\text{x}}_{1} ^{{\left( 1 \right)}} \left( {\text{k}} \right) + {\text{x}}_{1} ^{{\left( 1 \right)}} \left( {{\text{k}} + 1} \right)} \right] = \mathop \sum \limits_{{{\text{i}} = 1}}^{{{\text{n}} - 1}} {\text{b}}_{{\text{i}}} {\text{x}}_{{{\text{i}} + 1}} ^{{\left( 1 \right)}} \left( {{\text{k}} + 1} \right) $$
(8)

Assuming that

$$ Y = \left[ {x_{1} ^{{\left( 0 \right)}} \left( 2 \right),\;x_{1} ^{{\left( 0 \right)}} \left( 3 \right) \ldots x_{1} ^{{\left( 0 \right)}} (p)} \right] $$
(9)
$$ \beta = \left[ {a,\;b_{1} ,\;b_{2} \ldots b_{{n - 1}} } \right] $$
(10)
$$ B = \left[ {\begin{array}{*{20}c} { - 0.5\left[ {x_{1} ^{{\left( 1 \right)}} \left( 1 \right) + x_{1} ^{{\left( 1 \right)}} \left( 2 \right)} \right]} & {x_{2} ^{{\left( 1 \right)}} \left( 2 \right)} & {\begin{array}{*{20}c} \cdots & {x_{n} ^{{\left( 1 \right)}} \left( 2 \right)} \\ \end{array} } \\ { - 0.5\left[ {x_{1} ^{{\left( 1 \right)}} \left( 2 \right) + x_{1} ^{{\left( 1 \right)}} \left( 3 \right)} \right]} & {x_{2} ^{{\left( 1 \right)}} \left( 3 \right)} & {\begin{array}{*{20}c} \cdots & {x_{n} ^{{\left( 1 \right)}} \left( 3 \right)} \\ \end{array} } \\ {\begin{array}{*{20}c} \cdots \\ { - 0.5\left[ {x_{1} ^{{\left( 1 \right)}} \left( {p - 1} \right) + x_{1} ^{{\left( 1 \right)}} \left( p \right)} \right]} \\ \end{array} } & {\begin{array}{*{20}c} \cdots \\ {x_{2} ^{{\left( 1 \right)}} \left( p \right)} \\ \end{array} } & {\begin{array}{*{20}c} {\begin{array}{*{20}c} \cdots \\ \cdots \\ \end{array} } & {\begin{array}{*{20}c} \ldots \\ {x_{n} ^{{\left( 1 \right)}} \left( p \right)} \\ \end{array} } \\ \end{array} } \\ \end{array} } \right] $$
(11)

Therefore, Eq. (8) can be denoted as

$$ Y = B \cdot \beta $$
(12)

In the above equation, Y and B are known quantities and β is undetermined parameter. The parameters can be obtained by least square method. Namely

$$ \beta = \left( {B^{T} \cdot B} \right)^{{ - 1}} B^{T} \cdot Y $$
(13)

Substituting the obtained values back in Eq. (7), the analytical equation can be obtained as follows.

$$ x_{1} ^{1} \left( {k + 1} \right) = e^{{ - ak}} \left[ {x_{1} ^{{\left( 0 \right)}} \left( 1 \right) - \frac{1}{a}\mathop \sum \limits_{{i = 2}}^{n} b_{{i - 1}} x_{i} ^{{\left( 1 \right)}} \left( k \right)} \right] + \mathop \sum \limits_{{i = 2}}^{n} \frac{{b_{{i - 1}} }}{a}x_{i} ^{{\left( 1 \right)}} \left( k \right) $$
(14)

Natural gas reserves prediction in Sichuan Basin

The research methods explained in 2.2 Multi-cycle prediction model of natural gas reserves are used to carry out research on natural gas reserves prediction. The research process is shown in Fig. 2. The main research process is divided into 4 steps.

Fig. 2
figure 2

Flow chart of studies on natural gas reserves prediction

Full size image

Step 1: multi-cycle peak recognition (3.1.1 Multi-cycle peak recognition).

In the natural gas reserves growth curve with large fluctuation, multi-cycle peaks are screened out. There are three judgment processes: 1. peak attribute judgment, 2. symmetry judgment of left and right peaks, 3. half-peak elimination judgment. The screening process of multi-cycle peak is beneficial to analyzing the growth law of natural gas reserves.

Step 2: analysis of the multi-cycle peak change law (3.1.2 Analysis of the change law of multi-cyclic peak).

In order to highlight the change rule of reserves, the reserves baseline is drawn, which is to connect the peaks and valleys of multi-cycle peaks. The reserves growth stage is divided into four stages of natural gas exploration and development. Combined with the history of exploration and development and the change trend of reserves baseline, the growth law of reserves is analyzed. The comprehensive analysis of natural gas reserves change law (3.1.3 Comprehensive evaluation) provides a theoretical basis for the establishment of multi-cycle reserves prediction model.

Step 3: multi-cycle model parameters prediction (3.2.2 Parameters prediction of multi-cycle model).

The existing GM(1,3) gray prediction model is improved. Metabolism method is used to replace the most advanced data in the training sample with the predicted new data. The improved method is applied to multi-cycle model parameter prediction. The accuracy of parameter prediction is improved, which provides parameter support for multi-cycle natural gas reserves prediction.

Step 4: natural gas reserves prediction (3.2.3 Prediction of natural gas reserves based on multi-cycle Hubbert and Gauss models).

The predicted multi-cycle parameters are substituted in the improved multi-cycle Hubbert and Gauss models (3.2.1 Improvement of multi-cycle prediction model). The accuracy of reserves prediction results of different models is tested. The prediction curve with large correlation coefficient is selected as a result of reserves prediction. According to the result data of reserves prediction, the natural gas reserves prediction model of Sichuan Basin is comprehensively evaluated (3.2.4 Comprehensive evaluation).

Analysis of natural gas reserves growth law in Sichuan Basin

As shown in Fig. 1, the reserves curve shows obvious multi-cycle changes on the basis of the overall growth, and the cumulative proved reserves are approximately exponential growth. In order to analyze the change rule of Sichuan Basin new proved reserves, it is necessary to study the multi-cycle change of reserves in this Sichuan Basin. Furthermore, it can provide a theoretical basis for multi-cycle parameter prediction.

Multi-cycle peak recognition

The variation of natural gas reserves curve is complex with many numerical fluctuations. Therefore, it is impossible to treat all the natural gas reserves numerical increases and decreases as separate cyclical peaks. Find out all the extreme value points in the reserves value, and there are 18 maxima and 19 minima, which constitute 18 multi-cycle peaks (as shown in Fig. 3). Before determining the validity of these peaks, call them temporary peaks.

Fig. 3
figure 3

Analysis of multi-cycle changes of natural gas reserves

Full size image

The known natural gas reserves data in Fig. 3 are provided by the Exploration and Development Research Institute of Petro China Southwest Oil and Gas Field Company, who counts the annual proven reserves of natural gas in each natural gas area in the Sichuan Basin and sort out the reserves growth curve.

Now it is necessary to determine the validity of these 18 temporary peaks one by one. In order to describe the properties of the peak, some new basic concepts need to be defined before judging the validity. Each temporary peak has a peak point and two peak and valley points. And the adjacent peaks share a peak and valley point, which are, respectively, the rear valley point of the preceding peak and the front valley point of the following peak.

Suppose there is a temporary peak, and the ordinate difference between the peak top and the peak valley is \(D_{1}\), and the abscissa difference is \(D_{2}\). The peak slope \(L = \left| {D_{1} /D_{2} } \right|\). Each peak has a front peak slope \(L_{f}\) and a back peak slope \(L_{b}\). The larger one is the maximum peak slope \(L_{{\max }}\) and the other is the minimum peak slope \(L_{{\min }}\). The slope ratio is defined as \(L_{r} = L_{{\max }} /L_{{\min }}\). The absolute value difference between peak and valley is peak-to-peak value \({\text{PP}}\). Each peak has a front peak-to-peak value \({\text{PP}}_{f}\) and a back peak-to-peak value \({\text{PP}}_{b}\). The larger one is the maximum peak-to-peak value \({\text{PP}}_{{\max }}\), and the other one is the minimum peak-to-peak value \({\text{PP}}_{{\min }}\). The peak-to-peak ratio is defined as \({\text{PP}}_{r} = {\text{PP}}_{{\max }} /{\text{PP}}_{{\min }}\). Table 1 shows the calculated parameters of each temporary peak.

Table 1 Temporary peak parameter
Full size table

According to the above newly defined parameters, the threshold value of multi-cycle parameters is reasonably determined. And the effective judgment is carried out. The criteria for determining the multi-cycle peaks are divided into three steps.

Step 1: peak attribute judgment.

The peak should have the property of large fluctuation. The maximum slope \(L_{{\max }}\) and the maximum peak-to-peak value \({\text{PP}}_{{\max }}\) are used. Judging criteria: the maximum slope \(L_{{\max }} > 20\) and the maximum peak-to-peak value \({\text{PP}}_{{\max }} > 40\). When the two criteria are both met, the basic attribute of peak is obtained, and the temporary peak becomes “half peak” (the peak with half attribute of multi-cycle peak, that is, the undetermined peak). Then the second step is entered. Otherwise, it is not classified as the property of multi-cycle peak.

As shown in Fig. 4, the temporary peak 2 in the first stage, \({\text{PP}}_{{\max }} = 9.03 < 40\) and \(L_{{\max }} = 9.03 < 20\), did not meet the peak attribute requirements. It can also be seen from Fig. 6 that the changes at peak 2 are relatively gentle with minor fluctuation, which does not conform to the characteristic of large peak fluctuation. The other temporary peaks meet the judgment requirement and then carry on the next step judgment.

Fig. 4
figure 4

Peak property judgment

Full size image

Step 2: symmetry judgment of left and right peaks.

The left and right peaks of the cycle peak both should fluctuate up and down. The absolute slope ratio \(L_{r} = L_{{\max }} /L_{{\min }}\) and the absolute peak-to-peak ratio \({\text{PP}}_{r} = {\text{PP}}_{{\max }} /{\text{PP}}_{{\min _{r} }}\) are use. Judging criteria: absolute slope ratio \(L_{r} < 5\) and absolute peak-to-peak ratio \({\text{PP}}_{r} < 4\). When the two criteria are both met, the “half peak” can be classified as a cycle peak. Otherwise, the “half peak” is still a “half peak” rather than a cycle peak because only one side has a large fluctuation. The remaining “half peak” in this judgment step enters the third judgment step.

As shown in Fig. 5, the two parameters of “half peak” 3, 11 and 12 do not meet the judgment standard. And the absolute peak-to-peak ratio \({\text{PP}}_{r}\) of “half peak” 9, 15 and 16 does not meet the judgment standard. So these “half peaks” are still treated as “half peaks” and enters the third judgment step. The other “half peaks” meet the judgment requirements and can be classified as cycle peaks.

Fig. 5
figure 5

Symmetry judgment of left and right peaks

Full size image

Step 3: “half peaks” elimination judgment.

If there are several “half peaks” adjacent to each other and they can be connected from head to tail, it is considered that several “half peaks” coincide to form a complete peak. Otherwise, non-adjacent or can not be connected “half peaks” do not meet the coincidence standard and should be eliminated from the multi-cycle peak sequence.

In the 6 “half peaks” entering step 3, temporary peaks 3 and 9 are not adjacent to other half peaks. So they are eliminated from the multi-cycle peak sequence. “Half peak” 11 is adjacent to 12 and “half peak” 15 is adjacent to 16, forming two groups of coincident peaks. In order to facilitate the description of half peak coincidence criteria, the two adjacent peaks whose middle sequence is near the front and the back are called left peak and right peak, respectively. And their front and rear peaks are called left front peak, left posterior peak, right anterior peak and right posterior peak, respectively.

Coincidence standard: head to tail connection. Both the slope ratio of left posterior peak and the right anterior peak is within 0.5–2. And both the peak-to-peak ratio of left posterior peak and the right anterior peak is within 0.5–2. Through calculation, it is found that these two pairs of adjacent peaks accord with the standard and can be classified as two multi-cycle peaks, respectively. As shown in Fig. 6, the new cycle peak obtained by this step is the splice of two “half peaks” conforming to coincidence standards. Therefore, the top of the peak is flat and the peak width is large. So the combined peak type obtained by this step is trapezoidal peak. The historical reserves growth peaks in this basin include 4 kinds of peaks: ordinary peak, forward extended peak, trailing peak and trapezoidal peak.

Fig. 6
figure 6

Multi-cycle peak type

Full size image

The end of peak 18 is in the process of change, but it meets the judgment criteria and is regarded as a multi-cycle peak. After three steps of multi-cycle peak judgment, 18 temporary peaks are screened out, 3 unqualified peaks are eliminated, and 4 half peaks are combined into two new peaks. Finally, 13 multi-cycle peaks are obtained, as shown in Fig. 7.

Fig. 7
figure 7

Periodic analysis of annual newly added proved reserves

Full size image

Analysis of the change law of multi-cyclic peak

The change of natural gas reserves is characterized by multi-cycle, which is not conducive to directly analyzes the fluctuation of natural gas reserves over the years. First, draw the baseline below each peak. Connect the left and right valley points of each peak in the natural gas reserves curve, and the curve drawn is regarded as the baseline. The variation trend of reserves curve and baseline in the four stages from 1956 to 2018 is analyzed one by one. As shown in Fig. 7, in the first two stages of exploration and development, the reserves and the baseline show a trend of steady growth. In the latter two stages, the reserves curve tends to have obvious exponential growth and decline trend. And the change degree of the baseline also begins to increase dramatically.

The first stage: natural gas industry initial stage (1956–1977).

The natural gas reserves curve change is relatively stable in 1960–1964 and 1966–1971, so the curve over the period can be considered as the baseline. In half the time of this stage, the change of reserves showed a steady trend. No large reserves of natural gas reservoirs were found in the oil and gas exploration before 1960 in this period, and due to the discovery of a large number of limestone fracture and cave gas reservoirs in September 1964. Therefore, the annual proved reserves in the second year increased from 43.17 × 108m3 in the previous year to 425.14 × 108m3 in the previous year. However, as the major gas reservoirs discovered in this stage were small gas fields, the reserves quickly decreased to 45.81 × 108m3 in the second year, followed by a slow reserves growth process and two minor fluctuations in 1971–1976.

The second stage: steady natural gas growth stage (1978–2004).

As the natural gas exploration and development objects in this stage are mainly large- and medium-sized gas fields, the natural gas reserves begin to rise steadily. Due to the different discovery time of gas fields, the reserves increase steadily on the basis of the multi-cycle change pattern. Reserves baselines for this phase are also steadily rising, suggesting a breakthrough in gas exploration compared with the previous phase.

The third stage: rapid natural gas growth stage (2005–2010).

The annual average proved reserves showed a multi-cycle trend on the basis of rapid increase. Compared with the previous stage, the proved reserves in this stage increased significantly, and the baseline showed a trend of rising and then falling. This is because the exploration and development work in this stage was carried out at the same time, and the proved reserves were accelerated. No large gas fields were found in the later stage of this stage, so the reserves tended to decline after rising rapidly.

The fourth stage: natural gas development and expansion stage (2011–2018).

Due to the progress of exploration and development concepts and technologies, a large number of large natural gas reservoirs have been discovered in this stage. And the annual proved reserves are in a state of rapid increase. As the exploration rate of natural gas resources is only 10.34% at present, a large number of natural gas deposits are likely to be discovered in the future. The rapid increase of natural gas reserves in this stage will continue for a period of time.

During these four natural gas exploration stages, the curve of newly added reserves changes in 13 cycles. Since the fourth stage (i.e., the natural gas development and expansion stage) is not yet over, the 13th cycle (2016–2018) has experienced a change process of first increasing and then decreasing. Although the decreasing process may not be over yet, the 13th cycle can still be regarded as a complete cycle process.

Comprehensive evaluation

Analyze the change rule of the data of natural gas reserves, and combine with the progress of exploration and development in Sichuan Basin. It is found that the natural gas reserves growth is mainly affected by the number and scale of newly discovered gas fields. Factors influencing the natural gas exploration results mainly include technological progress, exploration area, market demand, oil and gas prices. Among them, the most influential factor is the development of oil and gas exploration concept, which mainly affects the growth of oil and gas reserves. And, there are also two criteria that define the four exploration and development stages: the progress of exploration concepts and the growth of reserves. Advances in exploration concepts are largely determined by the extent to which natural gas fields are explored. Therefore, the reserves growth are mainly determined by the natural gas exploration and development process.

By the end of 2018, the proven rate of natural gas in Sichuan Basin was only 10.34%, which means that the Sichuan Basin is still in the early stage of exploration and development with rapid increase of natural gas reserves. Moreover, since 2011, natural gas reserves have shown an obvious quasi-exponential growth trend. Therefore, the trend of reserves growth in a short period in the future can be predicted largely based on the trend of reserves data.

The analysis of gas reserves prediction result in Sichuan Basin

Improvement of multi-cycle prediction model

After obtaining all the unimodal parameters in the existing prediction model, the numerical value of each unimodal is added up in the whole time domain, which has great limitations in the multi-cycle prediction. As described in 3.1.1 Multi-cycle peak recognition, the historical reserves growth peaks in this basin include four peaks: ordinary peak, forward extended peak, trailing peak and trapezoidal peak. However, the existing multi-cycle prediction model is only applicable to the prediction of left–right symmetric ordinary peaks.

As shown in Eqs. (1) and (3), in unimodal prediction, the peaks obtained by the two models are symmetric about \(t = t_{m}\). It is also easy to maintain this feature in multi-peak prediction. When the distance between adjacent peaks is small and the peak change is slow, the prediction accuracy between adjacent peaks is easily affected by each other. The natural gas exploration in Sichuan Basin has made great progress, and the natural gas reserves have increased rapidly over the years. Therefore, it is not suitable to use the existing method to predict directly.

Therefore, the prediction method is improved from accumulation in time to splicing in space. That is, it no longer accumulates all multi-cycle peaks according to the existing method; instead, each independent peak is calculated separately. For example, the time span of the first peak is 1957–1960, and its calculation formula will not affect the calculation result of the second peak in 1964–1966. The specific steps are as follows:

Step 1: use the single-cycle model to analyze the change of each cycle peak.

During the single-cycle peak analysis, the model time domain length is only the time span of the peak, but not the whole time domain.

Step 2: use the GM(1,n) model predict multi-cycle peak parameters.

In view of the phenomenon about peak asymmetry, unimodal prediction around the peak should be predicted, respectively. Respectively, calculate left peak slope parameters \(b_{L}\), left peak standard deviation parameters \(s_{L}\), right peaks slope parameters \(b_{R}\) and right peak standard deviation parameters \(s_{R}\). The \(N_{m}\) and \(t_{m}\) are the left and right peak sharing parameter (the parameter meaning has been explained in 3.1.1 Multi-cycle peak recognition), which are also the sharing parameters of Hubbert and Gauss model.

Step 3: take the prediction value of multi-cycle in each period as the new prediction result without summation. Then test the validity of the model.

Use the improved Hubbert and Gauss models as described above to predict the growth trend of natural gas reserves. Compare the prediction results and select the prediction result which is more consistent with the original reserves data.

Parameters prediction of multi-cycle model

In this manuscript, the peak time \(t_{m}\), peak \(N_{m}\), model parameters \(b_{L}\), \(b_{R}\) and standard deviation parameters \(s_{L}\) and \(s_{R}\) are extracted from the reserves curve. Among them, \(b_{L}\) and \(b_{R}\) are parameters of multi-cycle Hubbert model, \(s_{L}\) and \(s_{R}\) are parameters of multi-cycle Gauss model, \(t_{m}\) and \(N_{m}\) are common parameters of both multi-cycle models. In the gray prediction GM(1,n) method, each multi-cycle model has four parameters of multi-cycle variables.

In the multi-cycle Hubbert model, there are four reserves parameters (\(N_{m}\), \(~t_{m}\), \(~b_{L}\), \(~b_{R}\)). For predicting each set of data, there are three groups of independent variables and a group of dependent variables. Therefore, in GM(1,n) method, n = 3. Through the gray prediction GM(1,3) method, the values of the remaining 3 parameters are used as the original data to predict the change rule of the 4th parameter in each model. The method is then used to predict the variation rules of the other three groups of parameters. Finally, use the newly predicted parameter values as the original data for the next round of prediction. The prediction principle of multi-cycle Gauss model is the same.

The above methods can easily cause the prediction results to deviate from the original law. Therefore, the prediction method is effectively improved by adopting the metabolic method. When a new set of parameter values is obtained, the most front set of data are removed from the training sample. The parameter values obtained by the new prediction are taken as the new training sample data, and a new round of prediction is made. This is repeated. In each new round of prediction, the data obtained from the new prediction are replaced by the front-end data in the training sample. In this way, the convergence of the predicted results can be guaranteed.

In 3.1.2 Analysis of the change law of multi-cyclic peak, the number of multi-peak cycles of the confirmed reserves change curve is 12.5, which is temporarily regarded as 13 cycles in the reserves prediction. To verify the effectiveness of the method, the data from 1956 to 2005 are used as the training sequence. And the reserves data from 2006 to 2018 are used as the test sequence. There are 9 cycles of natural gas reserves growth from 1956 to 2005. Table 2 shows the multi-peak parameters of the reserves curve during that period.

Table 2 Multi-cycle parameters of natural gas reserves
Full size table

Figure 8a and b is the prediction results of multi-cycle parameters of Hubbert model and Gauss model, respectively. The curve before 2018 is the known parameters shown in Table 2, and the time period curve after 2018 is the prediction trend. The peak time \(t_{m}\) is the abscissa, the model parameters b and s are the left ordinate, and the peak \(N_{m}\) is the right ordinate. It can be seen that the variation law of \(t_{m}\) − \(N_{m}\) conforms to the fluctuation characteristics of the original curve. And the predicted results of the two models \(N_{m}\) are numerically close. Parameter b and standard deviation parameter s, which represent the slope, vary greatly and also show the feature of multi-cycle. Among them, the left peak parameters (\(b_{L}\), \(s_{L}\)) and the corresponding right peak parameters (\(b_{R}\), \(s_{R}\)) have low coincidence rates, but the local maximum occurs at the same time.

Fig.8
figure 8

Multi-cycle parameter prediction results

Full size image

Prediction of natural gas reserves based on multi-cycle Hubbert and Gauss models

The multi-cycle parameters predicted in Fig. 8 are substituted into the multi-cycle Hubbert and Gauss calculation formulas. Then get the prediction results of multi-cycle Hubbert and Gauss model for annual newly added proved reserves (Fig. 9). The prediction period ends in 2030, the training period is from 1956 to 2005, and the test period is from 2006 to 2030.

Fig. 9
figure 9

Natural gas reserves prediction results of multi-cycle model

Full size image

The reserves prediction curve predicted by the two multi-cycle models is compared with the actual annual proved reserves curve. It is found that the fitting degree of Hubbert model prediction curve and actual reserves curve is higher. The Gauss model prediction curve changes more rapidly than the actual reserves curve. Therefore, Hubbert model is more suitable than Gauss model for natural gas reserves prediction research in Sichuan Basin.

After the existing peak of \(122992.3 \times 10^{8} {\text{m}}^{3}\) in 2014, the reserves prediction curve of Hubbert model decreased slowly and then increased sharply. Then the reserves curve reach to one new peak of \(4718.86 \times 10^{8} {\text{m}}^{3}\) in 2022 and another new peak of \(5383.64 \times 10^{8} {\text{m}}^{3}\) in 2027. The prediction results are in line with the trend that the natural gas reserves in Sichuan Basin will increase rapidly and continuously in the near future.

The ratio of the accumulative proved reserves to the amount of resources \(142551.52 \times 10^{8} {\text{m}}^{3}\) in a certain year is defined as the accumulative proved degree of that year. Figure 10 shows the prediction results of annual proved reserves and accumulative proved reserves of multi-cycle Hubbert and Gauss model. The peak time of the two prediction results and the accumulated proved reserves at the peak value are close to each other. And the two annual proved reserves growth curves are similar in shape. However, the annual proved reserves and accumulated reserves of Hubbert model are closer to the original data than the Gauss model.

Fig. 10
figure 10

Comprehensive prediction results of new proved natural gas reserves in Sichuan Basin

Full size image

The reserves data of Hubbert and Gauss models are analyzed by correlation with the actual reserves data. The correlation coefficient between the actual reserves data and the Hubbert model reservres data is 0.954 and that of Gauss model is 0.912. A larger correlation coefficient represents a better prediction result. Therefore, the predicted reserves data of Hubbert multi-cycle model are selected as the reserves prediction results. At the same time, it is also proved that the improved multi-cycle Hubbert model combined with GM(1,3) model is suitable for the prediction of natural gas reserves in Sichuan Basin.

Comprehensive evaluation

Table 3 shows the natural gas reserves prediction results of Hubbert and Gauss model. The reserves data are recorded according to the four stages of natural gas exploration and development in the Sichuan Basin (2.1 Overview of natural gas resources in Sichuan Basin). The average annual average proved reserves of the four stages grow rapidly in time (reserves prediction results of Hubbert and Gauss model). The cumulative proved degree of natural gas will reach 52.34% in 2030 (predicted by Hubbert model), which is close to the peak of 50% in the normal distribution. As a result, the Sichuan Basin will reach its peak of proven lifetime reserves in the next few years.

Table 3 Table of comprehensive prediction results of natural gas reserves in Sichuan Basin
Full size table

Conclusions

During 1956–2018, the natural gas reserves curve grows exponentially with multi-cycle characteristics. According to the trend of reserves growth, three criteria for judging multi-cycle peak are set: 1. peak attribute judgment, 2. symmetry judgment of left and right peaks, 3. half peak elimination judgment. Based on these three judgment criteria, the conclusion is concluded that the cycle number of natural gas reserves curve from 1956 to 2018 is 13. It is concluded that natural gas reserves will keep the trend of rapid growth in the short term. It provides theoretical basis for parameter prediction of multi-cycle model.

Metabolism modified GM(1,3) gray prediction model is used to predict multi-cycle parameter model. The prediction accuracy of multi-cycle parameters is improved. The measured multi-cycle parameters are substituted into the improved multi-cycle Hubbert and Gauss models. It is found that the correlation coefficient between the predicted reserves and the actual reserves data of Hubbert model is 0.954, which is higher than the Gauss model’s correlation coefficient which is 0.912. Therefore, the Hubbert model prediction data are selected as the reserves prediction results. The reserves prediction results show that the reserves curve reach to one new peak of \(4718.86 \times 10^{8} {\text{m}}^{3}\) in 2022 and another new peak of \(5383.64 \times 10^{8} {\text{m}}^{3}\) in 2027. Predictions show that the cumulative proven level of natural gas will reach 52.34% in 2030.