Gas Production Prediction Model of Volcanic Reservoir Based on Data-Driven Method

Abstract

Based on on-site construction experience, considering the time-varying characteristics of gas well quantity, production time, effective reservoir thickness, controlled reserves, reserve abundance, formation pressure, and the energy storage coefficient, a data-driven method was used to establish a natural gas production prediction model based on differential simulation theory. The calculation results showed that the average error between the actual production and predicted production was 12.49%, and the model determination coefficient was 0.99, indicating that the model can effectively predict natural gas production. Additionally, we observed that the influence of factors such as reserve abundance, the number of wells in operation, controlled reserves, the previous year’s gas production, formation pressure, the energy storage coefficient, effective matrix thickness, and annual production time on the annual gas production increases progressively as the F-values decrease. These insights are pivotal to a more profound understanding of gas production dynamics in volcanic reservoirs and are instrumental in optimizing stimulation treatments and enhancing resource recovery in such reservoirs and other unconventional hydrocarbon formations.

1. Introduction

Hydrocarbon production from unconventional reservoirs requires the integration of different technologies, including long lateral horizontal drilling and multi-stage, multi-cluster hydraulic fracture systems that activate natural fracture networks in unconventional formations [1,2]. Dark box approaches, which involve using data analytics techniques, have recently gained significant attention in many areas [3,4,5]. A variety of data analytics methods, such as machine learning [6], linear regression [7], and neural networks [8], are used for predicting gas and oil production.

Many scholars [9,10,11,12,13] from all over the world have applied linear regression methods in predicting gas and oil production. Zhou et al. [14] conducted a multiple regression model for 173 wells in the Marcellus Formation, predicting the one-year cumulative gas production as a function of the proppant mass, fracture fluid volume, number of stages, treatment rate, vertical depth, and lateral length. Grujic et al. [15] developed a predictive model for 172 wells in North America, forecasting oil, gas, and water production as a function of the volume, petrophysics, temperature, pressure, and geographical and completion parameters. Zhong et al. [16] predicted oil production in 476 wells in the Wolfcamp Formation as a function of the well and completion designs, using a multiple regression method and comparing it with other methods. Lolon et al. [17] predicted the cumulative oil production as a function of the stage of cementing, percentage of ceramic proppant, stage spacing, proppant intensity, water cut, fracturing fluid, and maximum treatment rate, by comparing some machine learning methods. Khanal et al. [18] forecasted the gas rate, cumulative gas, and condensate-to-gas ratio (CGR) using linear regression analysis for 335 simulations and 46 wells in the Eagle Ford Formation. Xue et al. [19] predicted the dynamic shale gas production rate as a function of hydraulic fracturing and geological properties using a multi-objective random forest regression method. Their sensitivity analysis revealed that the most influential parameters were the geological properties, including initial pressure and formation thickness. Johan et al. [20] employed a genetic algorithm to optimize completion designs, using several predictor variables, including the depth, lateral length, azimuth, total fluid, fluid intensity, total proppant, proppant intensity, and additional engineering features to account for the influence of neighboring wells and depletion.

The overall goal of this study is to develop a time variation characteristic prediction model that can accurately predict gas production in volcanic gas reservoirs. Specifically, this study aims to address the limitations of existing prediction models in processing the dynamic production data of unconventional gas reservoirs, which typically require a large amount of geological and permeability data. Traditional prediction methods require extensive numerical simulations. To overcome these challenges, we derived a new prediction model based on differential simulation theory and fitted the model using dimensionless gas production data obtained from volcanic gas reservoirs. In addition, we also calculated the error, coefficient of determination (R²), and F-value between actual and predicted gas production to verify the effectiveness and accuracy of the model. Through this study, we hope to provide a new and more effective tool for predicting gas production in unconventional gas reservoirs, thereby providing a scientific basis for a development strategy and production parameter adjustments in oil and gas fields.

2. Methodology

2.1. Workflow

Data-driven prediction, a method grounded in using existing data for analyzing and modeling to forecast future trends, outcomes, or events, is utilized in this paper to develop a production forecast model with time-varying characteristics. The methodology is outlined in a flow chart, depicted in Figure 1, and involves four sequential steps: initial data preprocessing, which entails the collection and normalization of gas field data to neutralize dimensional influences, followed by two stages of accumulation to reduce historical data randomness and prepare the data for modeling; the subsequent establishment of a multiple linear regression model based on a differential equation; parameter estimation using the least squares method to ensure predictive accuracy; and finally, data dimension recovery following two reduction processes to facilitate the computation of projected gas production. This comprehensive approach ensures the creation of a precise and reliable prediction model, essential in the accurate forecasting of gas field production.

2.2. Data Preprocessing

2.2.1. Data Collection

The indices for gas production data in volcanic gas reservoirs are bifurcated into two primary categories: the prediction index and the main control factor index. The prediction index is Q(t) (gas production), and the main control factor index mainly contains U₁(t) (number of wells in operation), U₂(t) (production time), U₃(t) (effective thickness of reservoir), U₄(t) (control reserves), U₅(t) (reserve abundance), U₆(t) (formation pressure), and U₇(t) (energy storage factor). Then, all the data indexes are collected in an Excel sheet.

2.2.2. Data Nondimensionalization

To standardize the data and eliminate the discrepancies arising from different physical units among various parameters, each data index should be normalized using the following equation:

$$x^{\left(0\right)}={\displaystyle \frac{x}{{\sum }_1^{n}x}}$$

(1)

where $x^{\left(0\right)}$ is the dimensionless value; $x$ is the raw data of the gas field data; and $n$ is the number of data.

2.2.3. Data Accumulation

To mitigate the randomness inherent to historical data and enhance the stability and reliability of the analysis, the dimensionless data are subjected to a first accumulation process. This is typically accomplished using the following equation:

$$x_i^{\left(1\right)}={\sum }_1^i{x}^{\left(0\right)}$$

(2)

where $x_i^{\left(1\right)}$ is the dimensionless value after the first accumulation; x is the value before normalization; and $n$ is the number of data.

Then, the dimensionless data after the first accumulation can be accumulated for the second time to form a fitting sample library as the following equation:

$$x_i^{\left(2\right)}={\sum }_1^i{x}^{\left(1\right)}$$

(3)

2.3. Modeling

Differential simulation theory, a specialized data-driven approach, posits that the accumulated time series data exhibit exponential variation characteristics, a property that can be mathematically substantiated [21]. Leveraging this theory, a yield data-driven model can be formulated to capture the underlying trends and dynamics of the data. The model development typically involves the following steps:

$${\left[Q^{\left(2\right)}\left(t\right)\right]}^{\prime }=a{Q}^{\left(2\right)}\left(t\right)+B{U}^{\left(2\right)}\left(t\right)$$

(4)

where ${\left[Q^{\left(2\right)}\left(t\right)\right]}^{\prime }$ represents the derivative with respect to time $t$; $U^{\left(2\right)}\left(t\right)=\left(\begin{array}{c}{U}_1^{\left(2\right)}\left(t\right)\\ U_2^{\left(2\right)}\left(t\right)\\ \dots \\ U_3^{\left(2\right)}\left(t\right)\end{array}\right)$; a and B are the non-identified parameters, and B is shown as follows:

$$B=(B_1,B_2,B_3,\dots ,B_7)$$

The above output data-driven Equation (4) is discretized by first-order approximation over time as follows:

$$Q_{k+1}^{\left(2\right)}=a{Q}^{\left(2\right)}\left(t\right)+B{U}^{\left(2\right)}\left(t\right)$$

(5)

Then, this Equation can be converted to the following equation:

$$Q_{k+1}^{\left(2\right)}=a{Q}_k^{\left(2\right)}+B_1{U}_{1k+1}^{\left(2\right)}+B_2{U}_{2k+1}^{\left(2\right)}+\dots +B_7{U}_{7k+1}^{\left(2\right)}$$

(6)

The derived parameters a and B, along with the flowchart detailing the least squares fitting process, are depicted in Figure 2. The dimensionless data, following the second accumulation, were fed into Equation (6) for comprehensive multiple linear regression analysis. Employing the least squares method, we calculated the sum of squared residuals between the actual observed values and those predicted by the model. This approach was instrumental in determining the optimal parameter values, thereby enhancing the model’s predictive accuracy and reliability.

2.4. Dimension Recovery

The predicted gas production is derived after two subtractions, utilizing the fitted equation. Specifically, the model is applied with the output data from year $k+1$ and the main control factors for year $k+1$. Thereafter, the dimensionless production data for year $k+1$, after two accumulations, are computed. Ultimately, the predicted gas production is ascertained through a series of two subtractions and dimension recovery processes, thereby yielding a precise forecast.

The first subtraction can be calculated as follows:

$$Q_{k+1}^{\left(1\right)}=Q_{k+1}^{\left(2\right)}-Q_k^{\left(2\right)}$$

(7)

where $Q_{k+1}^{\left(1\right)}$ is the production dimensionless data in $k+1$ after the first accumulation; $Q_{k+1}^{\left(2\right)}$ is the production dimensionless data in $k$ after the second accumulation; and $Q_k^{\left(2\right)}$ is the production dimensionless data in $k$ after the second accumulation.

The second subtraction can be calculated as follows:

$$Q_{k+1}^{\left(0\right)}=Q_{k+1}^{\left(1\right)}-Q_k^{\left(1\right)}$$

(8)

where $Q_{k+1}^{\left(0\right)}$ is the production dimensionless data in $k+1$; $Q_k^{\left(1\right)}$ is the production dimensionless data in $k$ after the first accumulation.

The predicted gas production can be obtained as follows:

$$Q_{k+1}=Q_{k+1}^{\left(0\right)}{\sum }_1^n{Q}_n$$

(9)

where $Q_{k+1}$ is the production data in $k+1$; ${\sum }_1^n{Q}_n$ is the sum of the gas production data.

3. Field Example

3.1. Data Preprocessing

The lithology of the targeted volcanic gas reservoir predominantly comprises acid tuff breccias, rhyolite, and andesite. The internal interlayering within the gas-bearing strata of the main rock mass is underdeveloped, with an interlayer density of merely 0.03 m/m, and the individual gas layers exhibit substantial thickness. A tuffaceous breccia septum, approximately 13 m thick, is present in the upper section of the gas layer, while a more substantial septum, about 170 m thick, is observed in the lower part. The fractures within the lower septum are relatively well developed, constituting 50.5% of the septum’s total thickness. The porosity within the reservoir fluctuates between 7.1% and 22.2%, with an average value of 14.4%. The permeability varies widely, from 0.005 mD to 836.000 mD, averaging at 0.844 mD, indicating a heterogeneous subsurface environment.

The data indices for the targeted volcanic rock well area were categorized into two distinct groups: the prediction index, which represents gas production, and the main control factor index, encompassing the number of wells in operation, production time, effective thickness of the reservoir, controlled reserves, reserve abundance, formation pressure, and energy storage factor. Subsequently, all index data were meticulously organized in annual chronological order, as presented in

Table 1. Annual natural gas production and main control factor index table.

Time /Year	Q(t) /×108 m3	U1(t) /Well	U2(t) /×103 h	U3(t) /m	U4(t) /×108 m3	U5(t) /×108 m3/km2	U6(t) /MPa	U7(t) /-
2008	0.05	3	0.89	91.30	18.4	15.51	46.61	9.43
2009	1.15	8	26.65	94.64	22.6	15.41	45.05	9.46
2010	1.94	11	50.49	86.99	25.2	15.25	43.49	8.81
2011	1.96	13	57.51	80.92	32.3	15.09	41.92	8.31
2012	1.77	15	68.45	81.56	36.0	14.94	40.36	7.83
2013	2.27	17	85.78	86.77	43.3	14.75	38.80	9.01
2014	2.35	15	86.99	91.13	69.3	14.56	37.30	9.83
2015	2.71	17	96.56	89.47	77.1	14.33	35.90	9.61
2016	2.5	20	98.74	89.93	82.9	14.12	34.50	9.95
2017	2.66	21	107.23	86.47	90.1	13.90	33.10	9.38
2018	2.39	21	97.42	87.00	95.7	13.70	31.80	9.50

.

In accordance with Equation (1), the original data presented in Table 1 were subjected to dimensionless processing, and the outcomes are detailed in

Table 2. Data after dimensionless processing.

Time/Year	${\mathit{Q}}^0\left(\mathit{t}\right)$	${\mathit{U}}_1^0\left(\mathit{t}\right)$	${\mathit{U}}_2^0\left(\mathit{t}\right)$	${\mathit{U}}_3^0\left(\mathit{t}\right)$	${\mathit{U}}_4^0\left(\mathit{t}\right)$	${\mathit{U}}_5^0\left(\mathit{t}\right)$	${\mathit{U}}_6^0\left(\mathit{t}\right)$	${\mathit{U}}_7^0\left(\mathit{t}\right)$
2008	0.002	0.017	0.001	0.087	0.027	0.089	0.101	0.085
2009	0.047	0.044	0.030	0.090	0.033	0.088	0.098	0.086
2010	0.079	0.061	0.057	0.083	0.036	0.087	0.095	0.080
2011	0.080	0.072	0.064	0.077	0.047	0.086	0.091	0.075
2012	0.072	0.083	0.077	0.077	0.052	0.085	0.088	0.071
2013	0.093	0.094	0.096	0.082	0.062	0.084	0.084	0.081
2014	0.096	0.083	0.097	0.087	0.100	0.083	0.081	0.089
2015	0.111	0.094	0.108	0.085	0.111	0.082	0.078	0.087
2016	0.102	0.110	0.110	0.085	0.120	0.081	0.075	0.090
2017	0.109	0.116	0.120	0.082	0.130	0.079	0.072	0.085
2018	0.098	0.116	0.109	0.083	0.138	0.078	0.069	0.086
2019	0.110	0.110	0.131	0.083	0.144	0.077	0.066	0.086

. The distribution of these dimensionless data points is graphically represented in Figure 3. The maximum values of the dimensionless data range from 0.089 to 0.144, indicating a variation of 38.19%. Conversely, the minimum values span from 0.001 to 0.077, reflecting a substantial change of 98.70%. It is evident that there are considerable disparities among the various data indices, highlighting the importance of normalization in comparative analyses.

Utilizing Equation (2), the original data were subjected to the first accumulation, with the outcomes documented in

Table 3. Data after the first cumulative processing.

Time/Year	${\mathit{Q}}^1\left(\mathit{t}\right)$	${\mathit{U}}_1^1\left(\mathit{t}\right)$	${\mathit{U}}_2^1\left(\mathit{t}\right)$	${\mathit{U}}_3^1\left(\mathit{t}\right)$	${\mathit{U}}_4^1\left(\mathit{t}\right)$	${\mathit{U}}_5^1\left(\mathit{t}\right)$	${\mathit{U}}_6^1\left(\mathit{t}\right)$	${\mathit{U}}_7^1\left(\mathit{t}\right)$
2008	0.002	0.017	0.001	0.087	0.027	0.089	0.101	0.085
2009	0.049	0.061	0.031	0.177	0.059	0.177	0.200	0.171
2010	0.128	0.122	0.087	0.259	0.096	0.264	0.294	0.250
2011	0.209	0.193	0.152	0.336	0.142	0.350	0.385	0.326
2012	0.281	0.276	0.228	0.413	0.194	0.435	0.473	0.396
2013	0.374	0.370	0.324	0.496	0.257	0.519	0.558	0.478
2014	0.470	0.453	0.422	0.582	0.357	0.603	0.639	0.567
2015	0.581	0.547	0.530	0.667	0.468	0.684	0.717	0.654
2016	0.683	0.657	0.640	0.753	0.588	0.765	0.792	0.744
2017	0.792	0.773	0.760	0.835	0.718	0.844	0.864	0.829
2018	0.890	0.890	0.869	0.917	0.856	0.923	0.934	0.914
2019	1.000	1.000	1.000	1.000	1.000	1.000	1.000	1.000

. The distribution of the dimensionless data after the first accumulation is depicted in Figure 4. Notably, the maximum values of the dimensionless data after this initial accumulation uniformly reach 1, signifying no variation. The minimum values, on the other hand, range from 0.001 to 0.101, representing a reduction of 99.01%. This observation underscores that the variability in the dimensionless data following the first accumulation is markedly diminished compared to the initial dimensionless data set, thereby illustrating the effectiveness of the accumulation process in stabilizing data fluctuations.

In advance of the modeling phase, the dimensionless data following the first accumulation were further accumulated using Equation (3), with the results presented in

Table 4. Data after the second cumulative processing.

Time/Year	${\mathit{Q}}^2\left(\mathit{t}\right)$	${\mathit{U}}_1^2\left(\mathit{t}\right)$	${\mathit{U}}_2^2\left(\mathit{t}\right)$	${\mathit{U}}_3^2\left(\mathit{t}\right)$	${\mathit{U}}_4^2\left(\mathit{t}\right)$	${\mathit{U}}_5^2\left(\mathit{t}\right)$	${\mathit{U}}_6^2\left(\mathit{t}\right)$	${\mathit{U}}_7^2\left(\mathit{t}\right)$
2008	0.002	0.017	0.001	0.087	0.027	0.089	0.101	0.085
2009	0.051	0.077	0.032	0.263	0.086	0.265	0.301	0.256
2010	0.180	0.199	0.119	0.522	0.181	0.529	0.595	0.507
2011	0.388	0.392	0.271	0.858	0.323	0.879	0.981	0.832
2012	0.669	0.669	0.499	1.272	0.518	1.314	1.454	1.229
2013	1.043	1.039	0.823	1.768	0.774	1.833	2.012	1.707
2014	1.514	1.492	1.245	2.350	1.131	2.436	2.651	2.273
2015	2.095	2.039	1.775	3.017	1.599	3.120	3.368	2.927
2016	2.778	2.696	2.415	3.770	2.186	3.886	4.160	3.671
2017	3.570	3.470	3.175	4.605	2.904	4.730	5.025	4.499
2018	4.460	4.359	4.044	5.522	3.759	5.653	5.958	5.414

. The distribution of the dimensionless data subsequent to the second accumulation is illustrated in Figure 5. The maximum values of the dimensionless data after this second accumulation range from 0.089 to 0.144, demonstrating a variation of 38.19%. The minimum values span from 0.001 to 0.077, indicating a substantial change of 98.70%. These refined data values, exhibiting reduced variability, are now suitable for inclusion in the modeling process.

3.2. Modeling

The dimensionless data after the second accumulation are input into Equation (6), then the a and B parameters are obtained after multiple linear regression simulation training by using the least squares method, and the results are shown below:

$$\begin{gathered} Q_{k+1}^{\left(2\right)}=-0.514Q_k^{\left(2\right)}+{{0.147U}_1^{\left(2\right)}}_{k+1}+{{1.631U}_2^{\left(2\right)}}_{k+1}-{14.301U_3^{\left(2\right)}}_{k+1} \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{{0.266U}_4^{\left(2\right)}}_{k+1}+{{0.441U}_5^{\left(2\right)}}_{k+1}+{{4.995U}_6^{\left(2\right)}}_{k+1} \\ +{{8.314U}_7^{\left(2\right)}}_{k+1} \end{gathered}$$

(10)

The gas production data and predicted values for the years 2008 to 2019, along with the projected production for 2020 to 2024, are graphically represented in Figure 6. A clear observation is that the predicted production aligns closely with the actual production across different years, indicating that the gas field reached a stable production phase by 2019. Additionally, Figure 7 illustrates the discrepancies between the actual and predicted gas production. The prediction error ranges from 0.46% to 30.47%, with an average error of 12.49%. It is evident that the error tends to diminish progressively with the extension of the mining period and the accumulation of production data, suggesting an improvement in the predictive accuracy over time.

3.3. Modeling Evaluation

3.3.1. Coefficient of Determination

The coefficient of determination ($R^2$) is a commonly used statistical indicator to evaluate the fit between the regression model and observed data [14,15,16]. It obviously indicates that the fitting model can explain the proportion of observed data variance, with a range of values between 0 and 1. The higher the $R^2$, the better the fit, and it can be calculated as follows:

$$\tilde{y}={\displaystyle \frac{1}{n}}\sum _{i=1}^n{y}_i$$

(11)

$$SST=\sum _i{\left(y_i-\tilde{y}\right)}^2$$

(12)

$$SSR=\sum _i{\left(y_i-y_i^{\prime }\right)}^2$$

(13)

$$R^2=1-{\displaystyle \frac{SSR}{SST}}$$

(14)

where $\tilde{y}$ is the average of observed value; $y_i^{\prime }$ is the predicted value; $y_i$ is the true value; $SST$ is the sum of squares in the real data; and $SSR$ is the sum of residual squares.

Utilizing the gas production and predicted data from 2009 to 2019, the coefficient of determination was calculated to be 0.99, as per Equations (11)–(14). This high value for the coefficient of determination signifies that the gas production data are closely aligned with the predictions of the fitting model, thereby confirming the model’s efficacy in capturing the production trends over the specified period.

3.3.2. Significance Testing

According to the relevant data of the established multiple linear regression fitting model, the significance of the regression equation is tested [22,23,24,25].

Hypothesis: $H_{0i}:{\beta }_i=0$, $iϵ\{1,2,\cdots ,12\}$, the sum of squared deviations are calculated as follows:

$$S_T^2=\sum _i^n(y_i-\overline{y})$$

(15)

$$S_R^2=\sum _i^n(\widehat{y}-\overline{y})$$

(16)

$$S_E^2=S_T^2-S_R^2$$

(17)

$$F={\displaystyle \frac{S_R^2/12}{S_R^2/12}}$$

(18)

F-value is a statistical measure of analysis of variance used to test whether the regression equation is significant [26,27,28,29,30]. Taking a significance level of $F=0.05$, the distribution table shows $F_{1-\alpha }\left(k,n-k-1\right)=F_{0.95}\left(\mathrm{729,1974}\right)=1<19.249$; thus, the hypothesis should be rejected.

The F-values for Q(t) (gas production), U₁(t) (number of wells in operation), U₂(t) (production time), U₃(t) (effective thickness of reservoir), U₄(t) (controlled reserves), U₅(t) (reserve abundance), U₆(t) (formation pressure), and U₇(t) (energy storage factor) are graphically represented in Figure 8. A clear trend is evident, with the F-values arranged in descending order as follows: U₂(t), U₃(t), U₇(t), U₆(t), Q(t), U₄(t), U₁(t), and U₅(t). This ranking underscores the increasing influence of the annual production time, effective thickness of the reservoir, energy storage factor, formation pressure, gas production, controlled reserves, number of wells in operation, and reserve abundance on the annual gas production. The visualization of these F-values provides a quantitative assessment of the relative significance of each factor in the predictive model [31].

4. Conclusions

In this study, standardization ensures data consistency and comparability, enhancing data quality for analysis. Cumulative operations reduce random fluctuations, highlighting long-term trends for better model predictions. The application of two accumulations of dimensionless gas production data effectively mitigates the inherent randomness of historical data, thereby establishing a robust fitting sample library.

Based on on-site construction experience, considering the time-varying characteristics of the gas well quantity, production time, effective reservoir thickness, controlled reserves, reserve abundance, formation pressure, and energy storage coefficient, a data-driven method was used to establish a natural gas production prediction model based on differential simulation theory. The calculation results showed that the average error between the actual production and predicted production was 12.49%, and the model determination coefficient was 0.99, indicating that the model can effectively predict natural gas production.

The F-values in descending order are U₂(t), U₃(t), U₇(t), U₆(t), Q(t), U₄(t), U₁(t), and U₅(t), indicating that the influence of the reserve abundance, number of wells in operation, controlled reserves, previous year’s gas production, formation pressure, energy storage coefficient, effective thickness of matrix, and annual production time on the annual gas production gradually increases.